mirror of
https://github.com/pfloos/quack
synced 2024-07-25 04:07:35 +02:00
qsUGF2
This commit is contained in:
parent
533e88e8fa
commit
29ed74cba0
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
178
src/GF/print_qsUGF2.f90
Normal file
178
src/GF/print_qsUGF2.f90
Normal file
@ -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)') ' <S**2> (exact) :',S2_exact
|
||||
write(*,'(A40,F13.6)') ' <S**2> :',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
|
||||
@ -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)
|
||||
|
||||
|
||||
309
src/GF/qsUGF2.f90
Normal file
309
src/GF/qsUGF2.f90
Normal file
@ -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
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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(*,*)
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user