diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90 index b0965fe..ee696cb 100644 --- a/src/GW/SRG_qsGW.f90 +++ b/src/GW/SRG_qsGW.f90 @@ -100,9 +100,9 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! Hello world write(*,*) - write(*,*)'************************************************' - write(*,*)'| Self-consistent SRG-qsGW calculation |' - write(*,*)'************************************************' + write(*,*)'***********************************' + write(*,*)'* Restricted SRG-qsGW Calculation *' + write(*,*)'***********************************' write(*,*) ! Warning diff --git a/src/GW/SRG_qsUGW.f90 b/src/GW/SRG_qsUGW.f90 new file mode 100644 index 0000000..0705847 --- /dev/null +++ b/src/GW/SRG_qsUGW.f90 @@ -0,0 +1,425 @@ +subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & + eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, & + dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) + +! Perform a quasiparticle self-consistent GW calculation + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dotest + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: BSE + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + double precision,intent(in) :: eta + logical,intent(in) :: regularize + + integer,intent(in) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + double precision,intent(in) :: ENuc + + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS(nspin) + + double precision,intent(in) :: EUHF + double precision,intent(in) :: eHF(nBas,nspin) + double precision,intent(in) :: cHF(nBas,nBas,nspin) + double precision,intent(in) :: PHF(nBas,nBas,nspin) + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: T(nBas,nBas) + double precision,intent(in) :: V(nBas,nBas) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_bbbb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + double precision,intent(inout):: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(inout):: dipole_int_bb(nBas,nBas,ncart) + +! Local variables + + logical :: dRPA + integer :: nSCF + integer :: nBasSq + integer :: ispin + integer :: ixyz + integer :: is + integer :: n_diis + integer :: nSa,nSb,nSt + double precision :: dipole(ncart) + + double precision :: ET(nspin) + double precision :: EV(nspin) + double precision :: EJ(nsp) + double precision :: EK(nspin) + double precision :: EcRPA(nspin) + double precision :: EcGM(nspin) + double precision :: EqsGW + double precision :: EcBSE(nspin) + double precision :: Conv + double precision :: rcond(nspin) + double precision,external :: trace_matrix + double precision,allocatable :: err_diis(:,:,:) + double precision,allocatable :: F_diis(:,:,:) + + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) + double precision,allocatable :: rho(:,:,:,:) + double precision,allocatable :: c(:,:,:) + double precision,allocatable :: cp(:,:,:) + double precision,allocatable :: eGW(:,:) + double precision,allocatable :: P(:,:,:) + double precision,allocatable :: F(:,:,:) + double precision,allocatable :: Fp(:,:,:) + double precision,allocatable :: J(:,:,:) + double precision,allocatable :: K(:,:,:) + double precision,allocatable :: SigC(:,:,:) + double precision,allocatable :: SigCp(:,:,:) + double precision,allocatable :: Z(:,:) + double precision,allocatable :: err(:,:,:) + +! Hello world + + write(*,*) + write(*,*)'*************************************' + write(*,*)'* Unrestricted SRG-qsGW Calculation *' + write(*,*)'*************************************' + write(*,*) + +! Warning + + write(*,*) '!! ERIs in MO basis will be overwritten in qsUGW !!' + write(*,*) + +! Stuff + + nBasSq = nBas*nBas + dRPA = .true. + +! TDA for W + + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if + +! TDA + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! Memory allocation + + nSa = nS(1) + nSb = nS(2) + nSt = nSa + nSb + + allocate(Aph(nSt,nSt),Bph(nSt,nSt),eGW(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin), & + F(nBas,nBas,nspin),Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin), & + SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin),Z(nBas,nspin),Om(nSt),XpY(nSt,nSt),XmY(nSt,nSt), & + rho(nBas,nBas,nSt,nspin),err(nBas,nBas,nspin),err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin)) + +! Initialization + + nSCF = -1 + n_diis = 0 + ispin = 1 + Conv = 1d0 + P(:,:,:) = PHF(:,:,:) + eGW(:,:) = eHF(:,:) + c(:,:,:) = cHF(:,:,:) + F_diis(:,:,:) = 0d0 + err_diis(:,:,:) = 0d0 + rcond(:) = 0d0 + +!------------------------------------------------------------------------ +! Main loop +!------------------------------------------------------------------------ + + do while(Conv > thresh .and. nSCF < maxSCF) + + ! Increment + + nSCF = nSCF + 1 + + ! Buid Hartree matrix + + do is=1,nspin + call Hartree_matrix_AO_basis(nBas,P(:,:,is),ERI_AO(:,:,:,:),J(:,:,is)) + end do + + ! Compute exchange part of the self-energy + + do is=1,nspin + call exchange_matrix_AO_basis(nBas,P(:,:,is),ERI_AO(:,:,:,:),K(:,:,is)) + end do + + !-------------------------------------------------- + ! AO to MO transformation of two-electron integrals + !-------------------------------------------------- + + do ixyz=1,ncart + call AOtoMO(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) + call AOtoMO(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) + end do + + ! 4-index transform for (aa|aa) block + + call AOtoMO_ERI_UHF(1,1,nBas,c,ERI_AO,ERI_aaaa) + + ! 4-index transform for (aa|bb) block + + call AOtoMO_ERI_UHF(1,2,nBas,c,ERI_AO,ERI_aabb) + + ! 4-index transform for (bb|bb) block + + call AOtoMO_ERI_UHF(2,2,nBas,c,ERI_AO,ERI_bbbb) + + ! Compute linear response + + call phULR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,Aph) + if(.not.TDA) call phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph) + + call phULR(TDA_W,nSa,nSb,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY) + + !----------------------! + ! Excitation densities ! + !----------------------! + + call UGW_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) + + !------------------------------------------------! + ! Compute self-energy and renormalization factor ! + !------------------------------------------------! + + if(regularize) then + do is=1,nspin + call GW_regularization(nBas,nC(is),nO(is),nV(is),nR(is),nSt,eGW(:,is),Om,rho(:,:,:,is)) + end do + end if + + call USRG_self_energy(eta,nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,EcGM,SigC,Z) + + ! Make correlation self-energy Hermitian and transform it back to AO basis + + do is=1,nspin + call MOtoAO(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) + end do + + ! Solve the quasi-particle equation + + do is=1,nspin + F(:,:,is) = Hc(:,:) + J(:,:,is) + J(:,:,mod(is,2)+1) + K(:,:,is) + SigCp(:,:,is) + end do + + ! Check convergence + + do is=1,nspin + err(:,:,is) = matmul(F(:,:,is),matmul(P(:,:,is),S(:,:))) - matmul(matmul(S(:,:),P(:,:,is)),F(:,:,is)) + end do + + if(nSCF > 1) Conv = maxval(abs(err)) + + ! DIIS extrapolation + + if(max_diis > 1) then + + n_diis = min(n_diis+1,max_diis) + do is=1,nspin + if(nO(is) > 1) call DIIS_extrapolation(rcond(is),nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis,is), & + F_diis(:,1:n_diis,is),err(:,:,is),F(:,:,is)) + end do + + end if + + ! Transform Fock matrix in orthogonal basis + + do is=1,nspin + Fp(:,:,is) = matmul(transpose(X(:,:)),matmul(F(:,:,is),X(:,:))) + end do + + ! Diagonalize Fock matrix to get eigenvectors and eigenvalues + + cp(:,:,:) = Fp(:,:,:) + do is=1,nspin + call diagonalize_matrix(nBas,cp(:,:,is),eGW(:,is)) + end do + + ! Back-transform eigenvectors in non-orthogonal basis + + do is=1,nspin + c(:,:,is) = matmul(X(:,:),cp(:,:,is)) + end do + + ! Back-transform self-energy + + do is=1,nspin + call AOtoMO(nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is)) + end do + + ! Compute density matrix + + do is=1,nspin + P(:,:,is) = matmul(c(:,1:nO(is),is),transpose(c(:,1:nO(is),is))) + end do + + !------------------------------------------------------------------------ + ! Compute total energy + !------------------------------------------------------------------------ + + ! Kinetic energy + + do is=1,nspin + ET(is) = trace_matrix(nBas,matmul(P(:,:,is),T(:,:))) + end do + + ! Potential energy + + do is=1,nspin + EV(is) = trace_matrix(nBas,matmul(P(:,:,is),V(:,:))) + end do + + ! Hartree energy + + EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) + EJ(2) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) & + + 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,1))) + EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2))) + + ! Exchange energy + + do is=1,nspin + EK(is) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,is),K(:,:,is))) + end do + + ! Total energy + + EqsGW = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(EK(:)) + + !------------------------------------------------------------------------ + ! Print results + !------------------------------------------------------------------------ + + call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int_AO,dipole) + call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,S,ENuc,ET,EV,EJ,EK,EcGM,EcRPA(ispin),EqsGW,SigCp,Z,dipole) + + end do +!------------------------------------------------------------------------ +! End main loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + end if + +! Deallocate memory + + deallocate(cp,P,F,Fp,J,K,SigC,SigCp,Z,Om,XpY,XmY,rho,err,err_diis,F_diis) + +! Perform BSE calculation + + if(BSE) then + + call UGW_phBSE(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & + S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,c,eGW,eGW,EcBSE) + + if(exchange_kernel) then + + EcBSE(1) = 0.5d0*EcBSE(1) + EcBSE(2) = 0.5d0*EcBSE(2) + + else + + EcBSE(2) = 0.0d0 + + end if + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@UHF correlation energy (spin-conserved) = ',EcBSE(1),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@UHF correlation energy (spin-flip) = ',EcBSE(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@UHF correlation energy = ',sum(EcBSE),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@UHF total energy = ',ENuc + EqsGW + sum(EcBSE),' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Compute the BSE correlation energy via the adiabatic connection + + if(doACFDT) then + + write(*,*) '--------------------------------------------------------------' + write(*,*) ' Adiabatic connection version of BSE@qsUGW correlation energy ' + write(*,*) '--------------------------------------------------------------' + write(*,*) + + if(doXBS) then + + write(*,*) '*** scaled screening version (XBS) ***' + write(*,*) + + end if + + call UGW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,spin_conserved,spin_flip, & + eta,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eGW,eGW,EcRPA) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@UHF correlation energy (spin-conserved) = ',EcRPA(1),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@UHF correlation energy (spin-flip) = ',EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@UHF correlation energy = ',sum(EcRPA),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@UHF total energy = ',ENuc + EqsGW + sum(EcRPA),' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if + + end if + +! Testing zone + + if(dotest) then + + call dump_test_value('U','qsGW correlation energy',EcRPA) + call dump_test_value('U','qsGW HOMOa energy',eGW(nO(1),1)) + call dump_test_value('U','qsGW LUMOa energy',eGW(nO(1)+1,1)) + call dump_test_value('U','qsGW HOMOa energy',eGW(nO(2),2)) + call dump_test_value('U','qsGW LUMOa energy',eGW(nO(2)+1,2)) + + end if + +end subroutine diff --git a/src/GW/UGW.f90 b/src/GW/UGW.f90 index ffcda9a..8cdb3e9 100644 --- a/src/GW/UGW.f90 +++ b/src/GW/UGW.f90 @@ -114,9 +114,9 @@ subroutine UGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre if(doqsGW) then call wall_time(start_GW) - call qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & - eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, & - dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) + call qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA, & + spin_conserved,spin_flip,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF, & + S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -132,7 +132,9 @@ subroutine UGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre if(doSRGqsGW) then call wall_time(start_GW) - print*,'Unrestricted version of SRG-qsGW not yet implemented! Sorry.' + call SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA, & + spin_conserved,spin_flip,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF, & + S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW diff --git a/src/GW/USRG_self_energy.f90 b/src/GW/USRG_self_energy.f90 new file mode 100644 index 0000000..4dbd0fe --- /dev/null +++ b/src/GW/USRG_self_energy.f90 @@ -0,0 +1,150 @@ +subroutine USRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) + +! Compute correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS,nspin) + +! Local variables + + integer :: ispin + integer :: i,j,a,b + integer :: p,q,r + integer :: m + double precision :: Dpim,Dqim,Dpam,Dqam,Diam + double precision :: t1,t2 + +! Output variables + + double precision,intent(out) :: EcGM + double precision,intent(out) :: SigC(nBas,nBas,nspin) + double precision,intent(out) :: Z(nBas,nspin) + +! Initialize + + SigC(:,:,:) = 0d0 + +!--------------------! +! SRG-GW self-energy ! +!--------------------! + + ! Occupied part of the correlation self-energy + + call wall_time(t1) + + !$OMP PARALLEL & + !$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Om) & + !$OMP PRIVATE(ispin,m,i,q,p,Dpim,Dqim) & + !$OMP DEFAULT(NONE) + !$OMP DO + do ispin=1,nspin + do q=nC(ispin)+1,nBas-nR(ispin) + do p=nC(ispin)+1,nBas-nR(ispin) + do m=1,nS + do i=nC(ispin)+1,nO(ispin) + Dpim = e(p,ispin) - e(i,ispin) + Om(m) + Dqim = e(q,ispin) - e(i,ispin) + Om(m) + SigC(p,q,ispin) = SigC(p,q,ispin) & + + rho(p,i,m,ispin)*rho(q,i,m,ispin)*(1d0-dexp(-eta*Dpim*Dpim)*dexp(-eta*Dqim*Dqim)) & + *(Dpim + Dqim)/(Dpim*Dpim + Dqim*Dqim) + end do + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(t2) + print *, "first loop", (t2-t1) + +! Virtual part of the correlation self-energy + + call wall_time(t1) + !$OMP PARALLEL & + !$OMP SHARED(SigC,rho,eta,nS,nC,nO,nR,nBas,e,Om) & + !$OMP PRIVATE(ispin,m,a,q,p,Dpam,Dqam) & + !$OMP DEFAULT(NONE) + !$OMP DO + do ispin=1,nspin + do q=nC(ispin)+1,nBas-nR(ispin) + do p=nC(ispin)+1,nBas-nR(ispin) + do m=1,nS + do a=nO(ispin)+1,nBas-nR(ispin) + Dpam = e(p,ispin) - e(a,ispin) - Om(m) + Dqam = e(q,ispin) - e(a,ispin) - Om(m) + SigC(p,q,ispin) = SigC(p,q,ispin) & + + rho(p,a,m,ispin)*rho(q,a,m,ispin)*(1d0-exp(-eta*Dpam*Dpam)*exp(-eta*Dqam*Dqam)) & + *(Dpam + Dqam)/(Dpam*Dpam + Dqam*Dqam) + end do + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(t2) + print *, "second loop", (t2-t1) + + +! Initialize + + Z(:,:) = 0d0 + + do ispin=1,nspin + do p=nC(ispin)+1,nBas-nR(ispin) + do i=nC(ispin)+1,nO(ispin) + do m=1,nS + Dpim = e(p,ispin) - e(i,ispin) + Om(m) + Z(p,ispin) = Z(p,ispin) - rho(p,i,m,ispin)**2*(1d0-dexp(-2d0*eta*Dpim*Dpim))/Dpim**2 + end do + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do ispin=1,nspin + do p=nC(ispin)+1,nBas-nR(ispin) + do a=nO(ispin)+1,nBas-nR(ispin) + do m=1,nS + Dpam = e(p,ispin) - e(a,ispin) - Om(m) + Z(p,ispin) = Z(p,ispin) - rho(p,a,m,ispin)**2*(1d0-dexp(-2d0*eta*Dpam*Dpam))/Dpam**2 + end do + end do + end do + end do + +! Compute renormalization factor from derivative of SigC + + Z(:,:) = 1d0/(1d0 - Z(:,:)) + +! Galitskii-Migdal correlation energy + + EcGM = 0d0 + do ispin=1,nspin + do i=nC(ispin)+1,nO(ispin) + do a=nO(ispin)+1,nBas-nR(ispin) + do m=1,nS + Diam = e(a,ispin) - e(i,ispin) + Om(m) + EcGM = EcGM - rho(a,i,m,ispin)*rho(a,i,m,ispin)*(1d0-exp(-2d0*eta*Diam*Diam))/Diam + end do + end do + end do + end do + +end subroutine diff --git a/src/GW/print_qsUGW.f90 b/src/GW/print_qsUGW.f90 index 3638c59..161e534 100644 --- a/src/GW/print_qsUGW.f90 +++ b/src/GW/print_qsUGW.f90 @@ -33,26 +33,20 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,Ov,ENuc,ET,EV,EJ,Ex,Ec logical :: dump_orb = .false. integer :: p integer :: ispin,ixyz - double precision :: HOMO(nspin) - double precision :: LUMO(nspin) - double precision :: Gap(nspin) + double precision :: eHOMO(nspin) + double precision :: eLUMO(nspin) + double precision :: Gap double precision :: Sz double precision :: Sx2,Sy2,Sz2 double precision,external :: trace_matrix ! HOMO and LUMO - do ispin=1,nspin - if(nO(ispin) > 0) then - HOMO(ispin) = eGW(nO(ispin),ispin) - LUMO(ispin) = eGW(nO(ispin)+1,ispin) - Gap(ispin) = LUMO(ispin) - HOMO(ispin) - else - HOMO(ispin) = 0d0 - LUMO(ispin) = eGW(1,ispin) - Gap(ispin) = 0d0 - end if + do ispin=1,nspin + eHOMO(ispin) = maxval(eGW(1:nO(ispin),ispin)) + eLUMO(ispin) = minval(eGW(nO(ispin)+1:nBas,ispin)) end do + Gap = minval(eLUMO) -maxval(eHOMO) Sz = 0.5d0*dble(nO(1) - nO(2)) Sx2 = 0.25d0*dble(nO(1) - nO(2)) + 0.5d0*nO(2) - 0.5d0*sum(matmul(transpose(c(:,1:nO(1),1)),matmul(Ov,c(:,1:nO(2),2)))**2) @@ -91,9 +85,9 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,Ov,ENuc,ET,EV,EJ,Ex,Ec write(*,'(2X,A14,F15.5)')'Convergence = ',Conv write(*,*)'----------------------------------------------------------------'// & '----------------------------------------------------------------' - write(*,'(2X,A60,F15.6,A3)') 'qsGW@UHF HOMO energy = ',maxval(HOMO)*HaToeV,' eV' - write(*,'(2X,A60,F15.6,A3)') 'qsGW@UHF LUMO energy = ',minval(LUMO)*HaToeV,' eV' - write(*,'(2X,A60,F15.6,A3)') 'qsGW@UHF HOMO-LUMO gap = ',(minval(LUMO)-maxval(HOMO))*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'qsGW@UHF HOMO energy = ',maxval(eHOMO)*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'qsGW@UHF LUMO energy = ',minval(eLUMO)*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'qsGW@UHF HOMO-LUMO gap = ',(minval(eLUMO)-maxval(eHOMO))*HaToeV,' eV' write(*,*)'----------------------------------------------------------------'// & '----------------------------------------------------------------' write(*,'(2X,A60,F15.6,A3)') ' qsGW@UHF total energy = ',ENuc + EqsGW,' au'