diff --git a/src/GF/G0F2.f90 b/src/GF/G0F2.f90 index bc1c7ab..cfc7750 100644 --- a/src/GF/G0F2.f90 +++ b/src/GF/G0F2.f90 @@ -57,7 +57,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO ! Frequency-dependent second-order contribution - call self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,ERI,SigC,Z,Ec) + call self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,ERI,SigC,Z) if(linearize) then @@ -71,6 +71,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO ! Print results + call MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF2,Ec) call print_G0F2(nBas,nO,eHF,SigC,eGF2,Z,ENuc,ERHF,Ec) ! Perform BSE2 calculation diff --git a/src/GF/evGF2.f90 b/src/GF/evGF2.f90 index eb5d775..2203533 100644 --- a/src/GF/evGF2.f90 +++ b/src/GF/evGF2.f90 @@ -77,7 +77,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet, ! Frequency-dependent second-order contribution - call self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z,Ec) + call self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) if(linearize) then @@ -93,6 +93,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet, ! Print results + call MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF2,Ec) call print_evGF2(nBas,nO,nSCF,Conv,eHF,SigC,Z,eGF2,ENuc,ERHF,Ec) ! DIIS extrapolation diff --git a/src/GF/print_qsGF2.f90 b/src/GF/print_qsGF2.f90 index 558e8b8..02c694d 100644 --- a/src/GF/print_qsGF2.f90 +++ b/src/GF/print_qsGF2.f90 @@ -20,12 +20,13 @@ subroutine print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigC double precision,intent(in) :: T(nBas,nBas),V(nBas,nBas) double precision,intent(in) :: J(nBas,nBas),K(nBas,nBas),F(nBas,nBas) double precision,intent(in) :: Z(nBas),SigC(nBas,nBas) + double precision,intent(in) :: Ec double precision,intent(in) :: dipole(ncart) ! Local variables integer :: q,ixyz,HOMO,LUMO - double precision :: Gap,ET,EV,EJ,Ex,Ec + double precision :: Gap,ET,EV,EJ,Ex double precision,external :: trace_matrix ! Output variables diff --git a/src/GF/print_qsUGF2.f90 b/src/GF/print_qsUGF2.f90 new file mode 100644 index 0000000..5b825e9 --- /dev/null +++ b/src/GF/print_qsUGF2.f90 @@ -0,0 +1,178 @@ +subroutine print_qsUGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,cGF2,PGF2,Ov,T,V,J,K, & + ENuc,ET,EV,EJ,Ex,Ec,EqsGF2,SigC,Z,dipole) + +! Print one-electron energies and other stuff for qsUGF2 + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nSCF + double precision,intent(in) :: ENuc + double precision,intent(in) :: ET(nspin) + double precision,intent(in) :: EV(nspin) + double precision,intent(in) :: EJ(nsp) + double precision,intent(in) :: Ex(nspin) + double precision,intent(in) :: Ec(nsp) + double precision,intent(in) :: EqsGF2 + double precision,intent(in) :: Conv + double precision,intent(in) :: thresh + double precision,intent(in) :: eHF(nBas,nspin) + double precision,intent(in) :: eGF2(nBas,nspin) + double precision,intent(in) :: cGF2(nBas,nBas,nspin) + double precision,intent(in) :: PGF2(nBas,nBas,nspin) + double precision,intent(in) :: Ov(nBas,nBas) + double precision,intent(in) :: T(nBas,nBas) + double precision,intent(in) :: V(nBas,nBas) + double precision,intent(in) :: J(nBas,nBas,nspin) + double precision,intent(in) :: K(nBas,nBas,nspin) + double precision,intent(in) :: SigC(nBas,nBas,nspin) + double precision,intent(in) :: Z(nBas,nspin) + double precision,intent(in) :: dipole(ncart) + +! Local variables + + integer :: p + integer :: ispin,ixyz + double precision :: HOMO(nspin) + double precision :: LUMO(nspin) + double precision :: Gap(nspin) + double precision :: S_exact,S2_exact + double precision :: S,S2 + double precision,external :: trace_matrix + +! HOMO and LUMO + + do ispin=1,nspin + if(nO(ispin) > 0) then + HOMO(ispin) = eGF2(nO(ispin),ispin) + LUMO(ispin) = eGF2(nO(ispin)+1,ispin) + Gap(ispin) = LUMO(ispin) - HOMO(ispin) + else + HOMO(ispin) = 0d0 + LUMO(ispin) = eGF2(1,ispin) + Gap(ispin) = 0d0 + end if + end do + + S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0) + S2 = S2_exact + nO(2) - sum(matmul(transpose(cGF2(:,1:nO(1),1)),matmul(Ov,cGF2(:,1:nO(2),2)))**2) + + S_exact = 0.5d0*dble(nO(1) - nO(2)) + S = -0.5d0 + 0.5d0*sqrt(1d0 + 4d0*S2) + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + if(nSCF < 10) then + write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation' + else + write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation' + endif + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(A1,A3,A1,A30,A1,A30,A1,A30,A1,A30,A1)') & + '|',' ','|','e_HF ','|','Sig_c ','|','Z ','|','e_QP ','|' + write(*,'(A1,A3,A1,2A15,A1,2A15,A1,2A15,A1,2A15,A1)') & + '|','#','|','up ','dw ','|','up ','dw ','|','up ','dw ','|','up ','dw ','|' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + + do p=1,nBas + write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & + '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigC(p,p,1)*HaToeV,SigC(p,p,2)*HaToeV,'|', & + Z(p,1),Z(p,2),'|',eGF2(p,1)*HaToeV,eGF2(p,2)*HaToeV,'|' + enddo + + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A10,I3)') 'Iteration ',nSCF + write(*,'(2X,A19,F15.5)')'max(|FPS - SPF|) = ',Conv + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A30,F15.6,A3)') 'qsUGF2 HOMO energy:',maxval(HOMO(:))*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'qsUGF2 LUMO energy:',minval(LUMO(:))*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'qsUGF2 HOMO-LUMO gap :',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV,' eV' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A30,F15.6,A3)') ' qsUGF2 total energy:',ENuc + EqsGF2 + sum(Ec(:)),' au' + write(*,'(2X,A30,F15.6,A3)') ' qsUGF2 exchange energy:',sum(Ex(:)),' au' + write(*,'(2X,A30,F15.6,A3)') ' qsUGF2 correlation energy:',sum(Ec(:)),' au' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,*) + +! Dump results for final iteration + + if(Conv < thresh) then + + write(*,*) + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40)') ' Summary ' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.10,A3)') ' One-electron energy: ',sum(ET(:)) + sum(EV(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' One-electron a energy: ',ET(1) + EV(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' One-electron b energy: ',ET(2) + EV(2),' au' + write(*,*) + write(*,'(A40,1X,F16.10,A3)') ' Kinetic energy: ',sum(ET(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Kinetic a energy: ',ET(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Kinetic b energy: ',ET(2),' au' + write(*,*) + write(*,'(A40,1X,F16.10,A3)') ' Potential energy: ',sum(EV(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Potential a energy: ',EV(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Potential b energy: ',EV(2),' au' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron energy: ',sum(EJ(:)) + sum(Ex(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron aa energy: ',EJ(1) + Ex(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy: ',EJ(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy: ',EJ(3) + Ex(2),' au' + write(*,*) + write(*,'(A40,1X,F16.10,A3)') ' Hartree energy: ',sum(EJ(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Hartree aa energy: ',EJ(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Hartree ab energy: ',EJ(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Hartree bb energy: ',EJ(3),' au' + write(*,*) + write(*,'(A40,1X,F16.10,A3)') ' Exchange energy: ',sum(Ex(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Exchange a energy: ',Ex(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Exchange b energy: ',Ex(2),' au' + write(*,*) + write(*,'(A40,1X,F16.10,A3)') ' Correlation energy: ',sum(Ec(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Correlation aa energy: ',Ec(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Correlation ab energy: ',Ec(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Correlation bb energy: ',Ec(3),' au' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.10,A3)') ' Electronic energy: ',EqsGF2 + sum(Ec(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au' + write(*,'(A40,1X,F16.10,A3)') ' qsUGF2 energy: ',ENuc + EqsGF2 + sum(Ec(:)),' au' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,F13.6)') ' S (exact) :',2d0*S_exact + 1d0 + write(*,'(A40,F13.6)') ' S :',2d0*S + 1d0 + write(*,'(A40,F13.6)') ' (exact) :',S2_exact + write(*,'(A40,F13.6)') ' :',S2 + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A45)') ' Dipole moment (Debye) ' + write(*,'(19X,4A10)') 'X','Y','Z','Tot.' + write(*,'(19X,4F10.6)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + + ! Print orbitals + + write(*,'(A50)') '-----------------------------------------' + write(*,'(A50)') 'qsUGF2 spin-up orbital coefficients ' + write(*,'(A50)') '-----------------------------------------' + call matout(nBas,nBas,cGF2(:,:,1)) + write(*,*) + write(*,'(A50)') '-----------------------------------------' + write(*,'(A50)') 'qsUGF2 spin-down orbital coefficients ' + write(*,'(A50)') '-----------------------------------------' + call matout(nBas,nBas,cGF2(:,:,2)) + write(*,*) + + endif + +end subroutine print_qsUGF2 diff --git a/src/GF/qsGF2.f90 b/src/GF/qsGF2.f90 index 7ae192b..5f810cf 100644 --- a/src/GF/qsGF2.f90 +++ b/src/GF/qsGF2.f90 @@ -137,7 +137,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, ! Compute self-energy and renormalization factor - call self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI_MO,SigC,Z,Ec) + call self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI_MO,SigC,Z) ! Make correlation self-energy Hermitian and transform it back to AO basis @@ -177,6 +177,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, ! Print results + call MP2(nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,eGF2,Ec) call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) call print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigCp,Z,EqsGF2,Ec,dipole) diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 new file mode 100644 index 0000000..93d2833 --- /dev/null +++ b/src/GF/qsUGF2.f90 @@ -0,0 +1,309 @@ +subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, & + 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 an unrestricted quasiparticle self-consistent GF2 calculation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + logical,intent(in) :: BSE + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: evDyn + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + double precision,intent(in) :: eta + + 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(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) + +! Local variables + + integer :: nSCF + integer :: nBasSq + integer :: ispin + integer :: is + integer :: n_diis + integer :: nS_aa,nS_bb,nS_sc + double precision :: dipole(ncart) + + double precision :: ET(nspin) + double precision :: EV(nspin) + double precision :: EJ(nsp) + double precision :: Ex(nspin) + double precision :: Ec(nspin) + double precision :: EqsGF2 + double precision :: EcBSE(nspin) + double precision :: EcAC(nspin) + double precision :: Conv + double precision :: rcond(nspin) + double precision,external :: trace_matrix + double precision,allocatable :: error_diis(:,:,:) + double precision,allocatable :: F_diis(:,:,:) + double precision,allocatable :: c(:,:,:) + double precision,allocatable :: cp(:,:,:) + double precision,allocatable :: eGF2(:,:) + 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 :: SigCm(:,:,:) + double precision,allocatable :: Z(:,:) + double precision,allocatable :: error(:,:,:) + +! Hello world + + write(*,*) + write(*,*)'**************************************************' + write(*,*)'| Self-consistent unrestricted qsGF2 calculation |' + write(*,*)'**************************************************' + write(*,*) + +! Warning + + write(*,*) '!! ERIs in MO basis will be overwritten in qsUGF2 !!' + write(*,*) + +! Stuff + + nBasSq = nBas*nBas + +! TDA + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(eGF2(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),SigCm(nBas,nBas,nspin), & + Z(nBas,nspin),error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin)) + +! Initialization + + nSCF = -1 + n_diis = 0 + ispin = 1 + Conv = 1d0 + P(:,:,:) = PHF(:,:,:) + eGF2(:,:) = eHF(:,:) + c(:,:,:) = cHF(:,:,:) + F_diis(:,:,:) = 0d0 + error_diis(:,:,:) = 0d0 + +!------------------------------------------------------------------------ +! Main loop +!------------------------------------------------------------------------ + + do while(Conv > thresh .and. nSCF < maxSCF) + + ! Increment + + nSCF = nSCF + 1 + + ! Buid Coulomb matrix + + do is=1,nspin + call Coulomb_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 + !-------------------------------------------------- + + ! 4-index transform for (aa|aa) block + + call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_aaaa) + + ! 4-index transform for (aa|bb) block + + call AOtoMO_integral_transform(1,1,2,2,nBas,c,ERI_AO,ERI_aabb) + + ! 4-index transform for (bb|bb) block + + call AOtoMO_integral_transform(2,2,2,2,nBas,c,ERI_AO,ERI_bbbb) + + !------------------------------------------------! + ! Compute self-energy and renormalization factor ! + !------------------------------------------------! + + call unrestricted_self_energy_GF2(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z) + + ! Make correlation self-energy Hermitian and transform it back to AO basis + + do is=1,nspin + SigCp(:,:,is) = 0.5d0*(SigC(:,:,is) + transpose(SigC(:,:,is))) + SigCm(:,:,is) = 0.5d0*(SigC(:,:,is) - transpose(SigC(:,:,is))) + end do + + do is=1,nspin + call MOtoAO_transform(nBas,S,c(:,:,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 + error(:,:,is) = matmul(F(:,:,is),matmul(P(:,:,is),S(:,:))) - matmul(matmul(S(:,:),P(:,:,is)),F(:,:,is)) + end do + + if(nSCF > 1) conv = maxval(abs(error(:,:,:))) + + ! DIIS extrapolation + + 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,error_diis(:,1:n_diis,is), & + F_diis(:,1:n_diis,is),error(:,:,is),F(:,:,is)) + end do + + ! Reset DIIS if required + + if(minval(rcond(:)) < 1d-15) n_diis = 0 + + ! 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),eGF2(:,is)) + end do + + ! Back-transform eigenvectors in non-orthogonal basis + + do is=1,nspin + c(:,:,is) = matmul(X(:,:),cp(:,:,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 + + ! Coulomb 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 + Ex(is) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,is),K(:,:,is))) + end do + + ! Correlation energy + + call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,eGF2,Ec) + + ! Total energy + + EqsGF2 = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + + !------------------------------------------------------------------------ + ! Print results + !------------------------------------------------------------------------ + + call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int_AO,dipole) + call print_qsUGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,P,S,T,V,J,K,ENuc,ET,EV,EJ,Ex,Ec,EqsGF2,SigCp,Z,dipole) + + enddo +!------------------------------------------------------------------------ +! End main loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + endif + +! Deallocate memory + + deallocate(cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,error,error_diis,F_diis) + +end subroutine qsUGF2 diff --git a/src/GF/self_energy_GF2.f90 b/src/GF/self_energy_GF2.f90 index 8ab7bad..173ecaf 100644 --- a/src/GF/self_energy_GF2.f90 +++ b/src/GF/self_energy_GF2.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z,Ec) +subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) ! Compute GF2 self-energy and its renormalization factor @@ -24,7 +24,6 @@ subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z,Ec) double precision,intent(out) :: SigC(nBas,nBas) double precision,intent(out) :: Z(nBas) - double precision,intent(out) :: Ec ! Initialize @@ -71,23 +70,4 @@ subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z,Ec) Z(:) = 1d0/(1d0 - Z(:)) -! Compute correlaiton energy - - Ec = 0d0 - - do j=nC+1,nO - do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR - - eps = eGF2(j) + eHF(i) - eHF(a) - eHF(b) - num = (2d0*ERI(j,i,a,b) - ERI(j,i,b,a))*ERI(j,i,a,b) - - Ec = Ec + num*eps/(eps**2 + eta**2) - - end do - end do - end do - end do - end subroutine self_energy_GF2 diff --git a/src/GF/self_energy_GF2_diag.f90 b/src/GF/self_energy_GF2_diag.f90 index e0611a0..c64d507 100644 --- a/src/GF/self_energy_GF2_diag.f90 +++ b/src/GF/self_energy_GF2_diag.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z,Ec) +subroutine self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) ! Compute diagonal part of the GF2 self-energy and its renormalization factor @@ -24,7 +24,6 @@ subroutine self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z,Ec) double precision,intent(out) :: SigC(nBas) double precision,intent(out) :: Z(nBas) - double precision,intent(out) :: Ec ! Initialize @@ -67,23 +66,4 @@ subroutine self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z,Ec) Z(:) = 1d0/(1d0 - Z(:)) -! Compute correlaiton energy - - Ec = 0d0 - - do j=nC+1,nO - do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR - - eps = eGF2(j) + eHF(i) - eHF(a) - eHF(b) - num = (2d0*ERI(j,i,a,b) - ERI(j,i,b,a))*ERI(j,i,a,b) - - Ec = Ec + num*eps/(eps**2 + eta**2) - - end do - end do - end do - end do - end subroutine self_energy_GF2_diag diff --git a/src/MP/UMP2.f90 b/src/MP/UMP2.f90 index 674eae9..acdfbd2 100644 --- a/src/MP/UMP2.f90 +++ b/src/MP/UMP2.f90 @@ -32,7 +32,7 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec) ! Output variables - double precision,intent(out) :: Ec + double precision,intent(out) :: Ec(nsp) ! Hello world @@ -72,6 +72,7 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec) enddo Ecaa = Edaa + Exaa + Ec(1) = Ecaa ! aabb block @@ -97,6 +98,7 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec) enddo Ecab = Edab + Exab + Ec(2) = Ecab ! bbbb block @@ -124,18 +126,18 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec) enddo Ecbb = Edbb + Exbb + Ec(3) = Ecbb ! Final flush Ed = Edaa + Edab + Edbb Ex = Exaa + Exab + Exbb - Ec = Ed + Ex write(*,*) write(*,'(A32)') '--------------------------' write(*,'(A32)') ' MP2 calculation ' write(*,'(A32)') '--------------------------' - write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',Ec + write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',sum(Ec(:)) write(*,'(A32,1X,F16.10)') ' alpha-alpha = ',Ecaa write(*,'(A32,1X,F16.10)') ' alpha-beta = ',Ecab write(*,'(A32,1X,F16.10)') ' beta-beta = ',Ecbb @@ -150,8 +152,8 @@ subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec) write(*,'(A32,1X,F16.10)') ' alpha-beta = ',Exab write(*,'(A32,1X,F16.10)') ' beta-beta = ',Exbb write(*,'(A32)') '--------------------------' - write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ', EHF + Ec - write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + Ec + write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ', EHF + sum(Ec(:)) + write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + sum(Ec(:)) write(*,'(A32)') '--------------------------' write(*,*)