diff --git a/src/GW/ccRG0W0.f90 b/src/GW/ccRG0W0.f90 index 1090288..0828ce5 100644 --- a/src/GW/ccRG0W0.f90 +++ b/src/GW/ccRG0W0.f90 @@ -31,8 +31,6 @@ 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(:) @@ -96,7 +94,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH do j=nC+1,nO do a=nO+1,nOrb-nR - del(i,j,a) = eHF(i) + eHF(j) - eHF(a) - eHF(p) + del(i,j,a) = eHF(i) + eHF(j) - eHF(a) end do end do @@ -106,7 +104,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH do a=nO+1,nOrb-nR do b=nO+1,nOrb-nR - del(b,a,i) = eHF(a) + eHF(b) - eHF(i) - eHF(p) + del(b,a,i) = eHF(a) + eHF(b) - eHF(i) end do end do @@ -130,27 +128,13 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH nSCF = nSCF + 1 - ! Compute intermediates x_2h1p and x_2p1h - - x = 0d0 - - do q=nC+1,nOrb-nR - do r=nC+1,nOrb-nR - do s=nC+1,nOrb-nR - - x = x + sqrt(2d0)*ERI(p,s,q,r)*amp(q,r,s) - - end do - end do - end do - ! Compute residual for 2h1p sector do i=nC+1,nO do j=nC+1,nO do a=nO+1,nOrb-nR - res(i,j,a) = sqrt(2d0)*ERI(p,a,i,j) + (del(i,j,a) - x)*amp(i,j,a) + res(i,j,a) = sqrt(2d0)*ERI(p,a,i,j) + (del(i,j,a) - eGW(p))*amp(i,j,a) do k=nC+1,nO do c=nO+1,nOrb-nR @@ -170,7 +154,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH 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) - x)*amp(b,a,i) + res(b,a,i) = sqrt(2d0)*ERI(p,i,b,a) + (del(b,a,i) - eGW(p))*amp(b,a,i) do k=nC+1,nO do c=nO+1,nOrb-nR diff --git a/src/GW/ccRG0W0.f90.bak b/src/GW/ccRG0W0.f90.bak new file mode 100644 index 0000000..789c950 --- /dev/null +++ b/src/GW/ccRG0W0.f90.bak @@ -0,0 +1,252 @@ +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