From 42f5439924860c2c96ce61fe72511bb4a89ac426 Mon Sep 17 00:00:00 2001 From: pfloos Date: Tue, 17 Sep 2024 19:12:45 +0200 Subject: [PATCH] ccG0W0 with minimal memory and DIIS --- src/GW/ccRG0W0.f90 | 103 ++++++++++------- src/GW/ccRG0W0.f90.bak | 252 ----------------------------------------- 2 files changed, 63 insertions(+), 292 deletions(-) delete mode 100644 src/GW/ccRG0W0.f90.bak diff --git a/src/GW/ccRG0W0.f90 b/src/GW/ccRG0W0.f90 index 0828ce5..d8a3a5a 100644 --- a/src/GW/ccRG0W0.f90 +++ b/src/GW/ccRG0W0.f90 @@ -31,10 +31,13 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH integer :: nSCF double precision :: Conv + double precision :: x + double precision,allocatable :: eGW(:) double precision,allocatable :: Z(:) double precision,allocatable :: del(:,:,:) + double precision,allocatable :: vec(:,:,:) double precision,allocatable :: res(:,:,:) double precision,allocatable :: amp(:,:,:) @@ -54,25 +57,48 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH ! Form energy denominator and guess amplitudes - allocate(del(nOrb,nOrb,nOrb)) - allocate(res(nOrb,nOrb,nOrb)) - allocate(amp(nOrb,nOrb,nOrb)) + allocate(del(nO,nV,nOrb)) + allocate(vec(nO,nV,nOrb)) + allocate(res(nO,nV,nOrb)) + allocate(amp(nO,nV,nOrb)) allocate(eGW(nOrb),Z(nOrb)) - allocate(r_diis(nOrb**3,max_diis)) - allocate(t_diis(nOrb**3,max_diis)) + allocate(r_diis(nO*nV*nOrb,max_diis)) + allocate(t_diis(nO*nV*nOrb,max_diis)) ! Initialization eGW(:) = eHF(:) - Z(:) = 1d0 + Z(:) = 1d0 + + ! Compute energy differences + + do i=nC+1,nO + do j=nC+1,nO + do a=1,nV-nR + + del(i,a,j) = eHF(i) + eHF(j) - eHF(nO+a) + + end do + end do + end do + + do i=nC+1,nO + do a=1,nV-nR + do b=1,nV-nR + + del(i,a,nO+b) = eHF(nO+a) + eHF(nO+b) - eHF(i) + + end do + end do + end do !-------------------------! ! Main loop over orbitals ! !-------------------------! - do p=nO,nO + do p=nC+1,nO+1 ! Initialization @@ -86,26 +112,25 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH amp(:,:,:) = 0d0 res(:,:,:) = 0d0 - del(:,:,:) = huge(1d0) - - ! Compute energy differences + ! Compute energy differences + do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nOrb-nR - - del(i,j,a) = eHF(i) + eHF(j) - eHF(a) - + do a=1,nV-nR + do j=nC+1,nO + + vec(i,a,j) = sqrt(2d0)*ERI(p,nO+a,i,j) + end do end do end do - + do i=nC+1,nO - do a=nO+1,nOrb-nR - do b=nO+1,nOrb-nR - - del(b,a,i) = eHF(a) + eHF(b) - eHF(i) - + do a=1,nV-nR + do b=1,nV-nR + + vec(i,a,nO+b) = sqrt(2d0)*ERI(p,i,nO+b,nO+a) + end do end do end do @@ -130,16 +155,16 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH ! Compute residual for 2h1p sector - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nOrb-nR + res(:,:,:) = vec(:,:,:) + (del(:,:,:) - eGW(p))*amp(:,:,:) - res(i,j,a) = sqrt(2d0)*ERI(p,a,i,j) + (del(i,j,a) - eGW(p))*amp(i,j,a) + do i=nC+1,nO + do a=1,nV-nR + do j=nC+1,nO do k=nC+1,nO - do c=nO+1,nOrb-nR + do c=1,nV-nR - res(i,j,a) = res(i,j,a) - 2d0*ERI(j,c,a,k)*amp(i,k,c) + res(i,a,j) = res(i,a,j) - 2d0*ERI(j,nO+c,nO+a,k)*amp(i,c,k) end do end do @@ -151,15 +176,13 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH ! Compute residual for 2p1h sector do i=nC+1,nO - do a=nO+1,nOrb-nR - do b=nO+1,nOrb-nR - - res(b,a,i) = sqrt(2d0)*ERI(p,i,b,a) + (del(b,a,i) - eGW(p))*amp(b,a,i) + do a=1,nV-nR + do b=1,nV-nR do k=nC+1,nO - do c=nO+1,nOrb-nR + do c=1,nV-nR - res(b,a,i) = res(b,a,i) + 2d0*ERI(a,k,i,c)*amp(b,c,k) + res(i,a,nO+b) = res(i,a,nO+b) + 2d0*ERI(nO+a,k,i,nO+c)*amp(k,c,nO+b) end do end do @@ -168,20 +191,20 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH end do end do - ! Check convergence + ! Check convergence Conv = maxval(abs(res)) ! Update amplitudes - amp(:,:,:) = amp(:,:,:) - res(:,:,:)/del(:,:,:) + amp(:,:,:) = amp(:,:,:) - res(:,:,:)/(del(:,:,:) - eHF(p)) ! DIIS extrapolation if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nOrb**3,nOrb**3,n_diis,r_diis,t_diis,res,amp) + call DIIS_extrapolation(rcond,nO*nV*nOrb,nO*nV*nOrb,n_diis,r_diis,t_diis,res,amp) end if @@ -189,11 +212,11 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH eGW(p) = eHF(p) - do q=nC+1,nOrb-nR - do r=nC+1,nOrb-nR - do s=nC+1,nOrb-nR + do i=nC+1,nO + do a=1,nV-nR + do q=nC+1,nOrb-nR - eGW(p) = eGW(p) + sqrt(2d0)*ERI(p,s,q,r)*amp(q,r,s) + eGW(p) = eGW(p) + vec(i,a,q)*amp(i,a,q) end do end do diff --git a/src/GW/ccRG0W0.f90.bak b/src/GW/ccRG0W0.f90.bak deleted file mode 100644 index 789c950..0000000 --- a/src/GW/ccRG0W0.f90.bak +++ /dev/null @@ -1,252 +0,0 @@ -subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) - -! CC-based GW module - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: maxSCF - double precision,intent(in) :: thresh - integer,intent(in) :: max_diis - - integer,intent(in) :: nBas - integer,intent(in) :: nOrb - integer,intent(in) :: nC - integer,intent(in) :: nO - integer,intent(in) :: nV - integer,intent(in) :: nR - double precision,intent(in) :: ENuc - double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nOrb) - double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) - -! Local variables - - integer :: p,q,r,s - integer :: i,j,k,l - integer :: a,b,c,d - - integer :: nSCF - double precision :: Conv - - double precision :: x - - double precision,allocatable :: eGW(:) - double precision,allocatable :: Z(:) - - double precision,allocatable :: del(:,:,:) - double precision,allocatable :: res(:,:,:) - double precision,allocatable :: amp(:,:,:) - - integer :: n_diis - double precision :: rcond - double precision,allocatable :: r_diis(:,:) - double precision,allocatable :: t_diis(:,:) - - -! Hello world - - write(*,*) - write(*,*)'*****************************' - write(*,*)'* CC-based G0W0 Calculation *' - write(*,*)'*****************************' - write(*,*) - -! Form energy denominator and guess amplitudes - - allocate(del(nO,nV,nOrb)) - allocate(res(nO,nV,nOrb)) - allocate(amp(nO,nV,nOrb)) - - allocate(eGW(nOrb),Z(nOrb)) - - allocate(r_diis(nO*nV*nOrb,max_diis)) - allocate(t_diis(nO*nV*nOrb,max_diis)) - -! Initialization - - eGW(:) = eHF(:) - Z(:) = 1d0 - - ! Compute energy differences - - do i=nC+1,nO - do j=nC+1,nO - do a=1,nV-nR - - del(i,a,j) = eHF(i) + eHF(j) - eHF(nO+a) - - end do - end do - end do - - do i=nC+1,nO - do a=1,nV-nR - do b=1,nV-nR - - del(i,a,nO+b) = eHF(nO+a) + eHF(nO+b) - eHF(i) - - end do - end do - end do - -!-------------------------! -! Main loop over orbitals ! -!-------------------------! - - do p=nO,nO - - ! Initialization - - Conv = 1d0 - nSCF = 0 - - n_diis = 0 - t_diis(:,:) = 0d0 - r_diis(:,:) = 0d0 - rcond = 0d0 - - amp(:,:,:) = 0d0 - res(:,:,:) = 0d0 - - !----------------------! - ! Loop over amplitudes ! - !----------------------! - - write(*,*) - write(*,*)'-------------------------------------------------------------' - write(*,*)'| CC-based G0W0 calculation |' - write(*,*)'-------------------------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & - '|','#','|','HF','|','G0W0','|','Conv','|' - write(*,*)'-------------------------------------------------------------' - - do while(Conv > thresh .and. nSCF < maxSCF) - - ! Increment - - nSCF = nSCF + 1 - - ! Compute residual for 2h1p sector - - do i=nC+1,nO - do a=1,nV-nR - do j=nC+1,nO - - res(i,a,j) = sqrt(2d0)*ERI(p,nO+a,i,j) + (del(i,a,j) - eGW(p))*amp(i,a,j) - - do k=nC+1,nO - do c=1,nV-nR - - res(i,a,j) = res(i,a,j) - 2d0*ERI(j,nO+c,nO+a,k)*amp(i,c,k) - - end do - end do - - end do - end do - end do - - ! Compute residual for 2p1h sector - - do i=nC+1,nO - do a=1,nV-nR - do b=1,nV-nR - - res(i,a,nO+b) = sqrt(2d0)*ERI(p,i,nO+b,nO+a) + (del(i,a,nO+b) - eGW(p))*amp(i,a,nO+b) - - do k=nC+1,nO - do c=1,nV-nR - - res(i,a,nO+b) = res(i,a,nO+b) + 2d0*ERI(nO+a,k,i,nO+c)*amp(k,c,nO+b) - - end do - end do - - end do - end do - end do - - ! Check convergence - - Conv = maxval(abs(res)) - - ! Update amplitudes - - amp(:,:,:) = amp(:,:,:) - res(:,:,:)/del(:,:,:) - - ! DIIS extrapolation - - if(max_diis > 1) then - - n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nO*nV*nOrb,nO*nV*nOrb,n_diis,r_diis,t_diis,res,amp) - - end if - - ! Compute quasiparticle energy - - eGW(p) = eHF(p) - - do i=nC+1,nOrb-nR - do a=1,nV-nR - do j=nC+1,nO - - eGW(p) = eGW(p) + sqrt(2d0)*ERI(p,nO+a,i,j)*amp(i,a,j) - - end do - end do - end do - - do i=nC+1,nOrb-nR - do a=1,nV-nR - do b=1,nV-nR - - eGW(p) = eGW(p) + sqrt(2d0)*ERI(p,i,nO+b,nO+a)*amp(i,a,nO+b) - - end do - end do - end do - - ! Dump results - - write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X)') & - '|',nSCF,'|',eHF(p)*HaToeV,'|',eGW(p)*HaToeV,'|',Conv,'|' - - end do - - write(*,*)'-------------------------------------------------------------' - !------------------------------------------------------------------------ - ! End of SCF loop - !------------------------------------------------------------------------ - - ! Did it actually converge? - - if(nSCF == maxSCF) then - - write(*,*) - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*)' Convergence failed ' - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*) - - stop - - end if - - write(*,*)'-------------------------------------------------------------------------------' - write(*,*)' CC-G0W0 calculation ' - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & - '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' - write(*,*)'-------------------------------------------------------------------------------' - - write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X)') & - '|',p,'|',eHF(p)*HaToeV,'|',(eGW(p)-eHF(p))*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' - write(*,*)'-------------------------------------------------------------------------------' - - end do - -end subroutine