4
1
mirror of https://github.com/pfloos/quack synced 2024-11-08 07:03:59 +01:00
This commit is contained in:
Pierre-Francois Loos 2023-07-04 09:27:35 +02:00
parent b701886e5c
commit 2fc348ebe4
7 changed files with 713 additions and 180 deletions

View File

@ -1,14 +0,0 @@
# Debuggin mode?
F
# Chemist notation for two-electron integral?
T
# Exposant of the Slater geminal
1.0
# One-electron integrals: Ov Kin Nuc
T T T
# Two-electron integrals: ERI F12 Yuk Erf
T F F F
# Three-electron integrals: Type1 Type2 Type3
F F F
# Four-electron integrals: Type1 Type2 Type3
F F F

View File

@ -15,5 +15,5 @@
# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
T F F F F F T F F F F F
# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
F F F T T F F F F F F T
# * unrestricted version available # * unrestricted version available

View File

@ -11,7 +11,7 @@
# GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W reg # GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W reg
256 0.00001 T 5 T 0.0 F F F F 256 0.00001 T 5 T 0.0 F F F F
# GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg
10 0.00001 T 5 T 0.1 T F 256 0.00001 T 5 T 0.0 T F
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS
F T T F T T
# BSE: BSE dBSE dTDA evDyn ppBSE BSE2 # BSE: BSE dBSE dTDA evDyn ppBSE BSE2

View File

@ -0,0 +1,94 @@
subroutine GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC)
! Compute correlation part of the self-energy for GTeh
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Om(nS)
double precision,intent(in) :: rhoL(nBas,nBas,nS)
double precision,intent(in) :: rhoR(nBas,nBas,nS)
! Local variables
integer :: i,j,a,b
integer :: p,q,r
integer :: m
double precision :: eps
! Output variables
double precision,intent(out) :: EcGM
double precision,intent(out) :: SigC(nBas,nBas)
! Initialize
SigC(:,:) = 0d0
!----------------!
! GW self-energy !
!----------------!
! Occupied part of the correlation self-energy
!$OMP PARALLEL &
!$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Om) &
!$OMP PRIVATE(m,i,q,p,eps) &
!$OMP DEFAULT(NONE)
!$OMP DO
do q=nC+1,nBas-nR
do p=nC+1,nBas-nR
do m=1,nS
do i=nC+1,nO
eps = e(p) - e(i) + Om(m)
SigC(p,q) = SigC(p,q) + rhoL(i,p,m)*rhoR(i,q,m)*eps/(eps**2 + eta**2)
end do
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! Virtual part of the correlation self-energy
!$OMP PARALLEL &
!$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Om) &
!$OMP PRIVATE(m,a,q,p,eps) &
!$OMP DEFAULT(NONE)
!$OMP DO
do q=nC+1,nBas-nR
do p=nC+1,nBas-nR
do m=1,nS
do a=nO+1,nBas-nR
eps = e(p) - e(a) - Om(m)
SigC(p,q) = SigC(p,q) + rhoL(p,a,m)*rhoR(q,a,m)*eps/(eps**2 + eta**2)
end do
end do
end do
end do
!$OMP END DO
!$OMP END PARALLEL
! Galitskii-Migdal correlation energy
EcGM = 0d0
do m=1,nS
do a=nO+1,nBas-nR
do i=nC+1,nO
eps = e(a) - e(i) + Om(m)
EcGM = EcGM - rhoL(i,a,m)*rhoR(i,a,m)*eps/(eps**2 + eta**2)
end do
end do
end do
end subroutine

120
src/GT/print_qsGTeh.f90 Normal file
View File

@ -0,0 +1,120 @@
subroutine print_qsGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
! Print one-electron energies and other stuff for qsGTeh
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nO
integer,intent(in) :: nSCF
double precision,intent(in) :: ENuc
double precision,intent(in) :: ET
double precision,intent(in) :: EV
double precision,intent(in) :: EJ
double precision,intent(in) :: Ex
double precision,intent(in) :: EcGM
double precision,intent(in) :: EcRPA(nspin)
double precision,intent(in) :: Conv
double precision,intent(in) :: thresh
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGT(nBas)
double precision,intent(in) :: c(nBas)
double precision,intent(in) :: SigC(nBas,nBas)
double precision,intent(in) :: Z(nBas)
double precision,intent(in) :: dipole(ncart)
! Local variables
integer :: p,ixyz,HOMO,LUMO
double precision :: Gap
double precision,external :: trace_matrix
! Output variables
double precision,intent(out) :: EqsGT
! HOMO and LUMO
HOMO = nO
LUMO = HOMO + 1
Gap = eGT(LUMO)-eGT(HOMO)
! Compute energies
! Dump results
write(*,*)'-------------------------------------------------------------------------------'
if(nSCF < 10) then
write(*,'(1X,A21,I1,A3,I1,A12)')' Self-consistent qsG',nSCF,'Teh',nSCF,' calculation'
else
write(*,'(1X,A21,I2,A3,I2,A12)')' Self-consistent qsG',nSCF,'Teh',nSCF,' calculation'
endif
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_T (eV)','|','Z','|','e_QP (eV)','|'
write(*,*)'-------------------------------------------------------------------------------'
do p=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A10,I3)') 'Iteration ',nSCF
write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
write(*,*)'-------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') 'qsGTeh HOMO energy:',eGT(HOMO)*HaToeV,' eV'
write(*,'(2X,A30,F15.6,A3)') 'qsGTeh LUMO energy:',eGT(LUMO)*HaToeV,' eV'
write(*,'(2X,A30,F15.6,A3)') 'qsGTeh HOMO-LUMO gap :',Gap*HaToeV,' eV'
write(*,*)'-------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') ' qsGTeh total energy:',ENuc + EqsGT,' au'
write(*,'(2X,A30,F15.6,A3)') ' qsGTeh exchange energy:',Ex,' au'
write(*,'(2X,A30,F15.6,A3)') ' GM@qsGTeh correlation energy:',EcGM,' au'
write(*,'(2X,A30,F15.6,A3)') 'ppRPA@qsGTeh correlation energy:',sum(EcRPA(:)),' au'
write(*,*)'-------------------------------------------'
write(*,*)
! Dump results for final iteration
if(Conv < thresh) then
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' Summary '
write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10,A3)') ' One-electron energy: ',ET + EV,' au'
write(*,'(A32,1X,F16.10,A3)') ' Kinetic energy: ',ET,' au'
write(*,'(A32,1X,F16.10,A3)') ' Potential energy: ',EV,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10,A3)') ' Two-electron energy: ',EJ + Ex,' au'
write(*,'(A32,1X,F16.10,A3)') ' Hartree energy: ',EJ,' au'
write(*,'(A32,1X,F16.10,A3)') ' Exchange energy: ',Ex,' au'
write(*,'(A32,1X,F16.10,A3)') ' Correlation energy: ',EcGM,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10,A3)') ' Electronic energy: ',EqsGT,' au'
write(*,'(A32,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au'
write(*,'(A32,1X,F16.10,A3)') ' qsGTeh energy: ',ENuc + EqsGT,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A35)') ' Dipole moment (Debye) '
write(*,'(10X,4A10)') 'X','Y','Z','Tot.'
write(*,'(10X,4F10.6)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD
write(*,'(A50)') '-----------------------------------------'
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGTeh MO coefficients'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,c)
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGTeh MO energies'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,1,eGT)
write(*,*)
endif
end subroutine

341
src/GT/qsGTeh.f90 Normal file
View File

@ -0,0 +1,341 @@
subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA, &
dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! Perform a quasiparticle self-consistent GTeh 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) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS
logical,intent(in) :: BSE
logical,intent(in) :: BSE2
logical,intent(in) :: TDA_T
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: singlet
logical,intent(in) :: triplet
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
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
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_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_MO(nBas,nBas,ncart)
! Local variables
integer :: nSCF
integer :: nBasSq
integer :: ispin
integer :: n_diis
double precision :: ET
double precision :: EV
double precision :: EJ
double precision :: Ex
double precision :: EqsGT
double precision :: EcRPA
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: EcGM
double precision :: Conv
double precision :: rcond
double precision,external :: trace_matrix
double precision :: dipole(ncart)
logical :: print_T = .true.
double precision,allocatable :: error_diis(:,:)
double precision,allocatable :: F_diis(:,:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rhoL_RPA(:,:,:)
double precision,allocatable :: rhoR_RPA(:,:,:)
double precision,allocatable :: c(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: eGT(:)
double precision,allocatable :: eOld(:)
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 qsGTeh calculation |'
write(*,*)'************************************************'
write(*,*)
! Warning
write(*,*) '!! ERIs in MO basis will be overwritten in qsGTeh !!'
write(*,*)
! Stuff
nBasSq = nBas*nBas
! TDA for T
if(TDA_T) 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
allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), &
OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS), &
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Initialization
nSCF = -1
n_diis = 0
ispin = 2
Conv = 1d0
P(:,:) = PHF(:,:)
eGT(:) = eHF(:)
eOld(:) = eHF(:)
c(:,:) = cHF(:,:)
F_diis(:,:) = 0d0
error_diis(:,:) = 0d0
rcond = 0d0
!------------------------------------------------------------------------
! Main loop
!------------------------------------------------------------------------
do while(Conv > thresh .and. nSCF <= maxSCF)
! Increment
nSCF = nSCF + 1
! Buid Coulomb matrix
call Coulomb_matrix_AO_basis(nBas,P,ERI_AO,J)
! Compute exchange part of the self-energy
call exchange_matrix_AO_basis(nBas,P,ERI_AO,K)
! AO to MO transformation of two-electron integrals
call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_MO)
! Compute linear response
call linear_response(ispin,.false.,TDA_T,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO, &
EcRPA,OmRPA,XpY_RPA,XmY_RPA)
if(print_T) call print_excitation('RPA@qsGTeh ',ispin,nS,OmRPA)
! Compute correlation part of the self-energy
call GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA)
if(regularize) then
! call regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC)
! call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,Z)
else
call GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC)
call GTeh_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,Z)
endif
! Make correlation self-energy Hermitian and transform it back to AO basis
SigCp = 0.5d0*(SigC + transpose(SigC))
SigCm = 0.5d0*(SigC - transpose(SigC))
call MOtoAO_transform(nBas,S,c,SigCp)
! Solve the quasi-particle equation
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:)
! Compute commutator and convergence criteria
error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
! DIIS extrapolation
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
end if
! Diagonalize Hamiltonian in AO basis
Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGT)
c = matmul(X,cp)
SigCp = matmul(transpose(c),matmul(SigCp,c))
! Compute new density matrix in the AO basis
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Save quasiparticles energy for next cycle
Conv = maxval(abs(error))
eOld(:) = eGT(:)
!------------------------------------------------------------------------
! Compute total energy
!------------------------------------------------------------------------
! Kinetic energy
ET = trace_matrix(nBas,matmul(P,T))
! Potential energy
EV = trace_matrix(nBas,matmul(P,V))
! Coulomb energy
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
! Exchange energy
Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
! Total energy
EqsGT = ET + EV + EJ + Ex
! Print results
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigCp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
enddo
!------------------------------------------------------------------------
! End main loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF+1) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
stop
endif
! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,OmRPA,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA,error,error_diis,F_diis)
! Perform BSE calculation
! if(BSE) then
! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, &
! eGW,eGW,EcBSE)
! if(exchange_kernel) then
! EcBSE(1) = 0.5d0*EcBSE(1)
! EcBSE(2) = 1.5d0*EcBSE(2)
! end if
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1)
! write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy (triplet) =',EcBSE(2)
! write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2)
! write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
! Compute the BSE correlation energy via the adiabatic connection
! if(doACFDT) then
! write(*,*) '------------------------------------------------------'
! write(*,*) 'Adiabatic connection version of BSE correlation energy'
! write(*,*) '------------------------------------------------------'
! write(*,*)
! if(doXBS) then
! write(*,*) '*** scaled screening version (XBS) ***'
! write(*,*)
! end if
! call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC)
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (singlet) =',EcAC(1)
! write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (triplet) =',EcAC(2)
! write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy =',EcAC(1) + EcAC(2)
! write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW total energy =',ENuc + EqsGW + EcAC(1) + EcAC(2)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
! end if
! end if
end subroutine

View File

@ -40,19 +40,6 @@ program QuAcK
logical :: exchange_kernel logical :: exchange_kernel
logical :: doXBS logical :: doXBS
! integer :: nShell
! integer,allocatable :: TotAngMomShell(:)
! integer,allocatable :: KShell(:)
! double precision,allocatable :: CenterShell(:,:)
! double precision,allocatable :: DShell(:,:)
! double precision,allocatable :: ExpShell(:,:)
! integer,allocatable :: max_ang_mom(:)
! double precision,allocatable :: min_exponent(:,:)
! double precision,allocatable :: max_exponent(:)
integer :: TrialType
double precision,allocatable :: cTrial(:),gradient(:),hessian(:,:)
double precision,allocatable :: S(:,:) double precision,allocatable :: S(:,:)
double precision,allocatable :: T(:,:) double precision,allocatable :: T(:,:)
double precision,allocatable :: V(:,:) double precision,allocatable :: V(:,:)
@ -80,28 +67,14 @@ program QuAcK
double precision :: start_HF ,end_HF ,t_HF double precision :: start_HF ,end_HF ,t_HF
double precision :: start_stab ,end_stab ,t_stab double precision :: start_stab ,end_stab ,t_stab
double precision :: start_KS ,end_KS ,t_KS double precision :: start_KS ,end_KS ,t_KS
double precision :: start_MOM ,end_MOM ,t_MOM
double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO
double precision :: start_CCD ,end_CCD ,t_CCD double precision :: start_CC ,end_CC ,t_CC
double precision :: start_DCD ,end_DCD ,t_DCD double precision :: start_CI ,end_CI ,t_CI
double precision :: start_CCSD ,end_CCSD ,t_CCSD
double precision :: start_CIS ,end_CIS ,t_CIS
double precision :: start_CID ,end_CID ,t_CID
double precision :: start_CISD ,end_CISD ,t_CISD
double precision :: start_FCI ,end_FCI ,t_FCI
double precision :: start_RPA ,end_RPA ,t_RPA double precision :: start_RPA ,end_RPA ,t_RPA
double precision :: start_ADC ,end_ADC ,t_ADC double precision :: start_GF ,end_GF ,t_GF
double precision :: start_GF2 ,end_GF2 ,t_GF2 double precision :: start_GW ,end_GW ,t_GW
double precision :: start_GF3 ,end_GF3 ,t_GF3 double precision :: start_GT ,end_GT ,t_GT
double precision :: start_G0W0 ,end_G0W0 ,t_G0W0 double precision :: start_MP ,end_MP ,t_MP
double precision :: start_evGW ,end_evGW ,t_evGW
double precision :: start_qsGW ,end_qsGW ,t_qsGW
double precision :: start_ufGW ,end_ufGW ,t_ufGW
double precision :: start_G0T0 ,end_G0T0 ,t_G0T0
double precision :: start_evGT ,end_evGT ,t_evGT
double precision :: start_qsGT ,end_qsGT ,t_qsGT
double precision :: start_MP2 ,end_MP2 ,t_MP2
double precision :: start_MP3 ,end_MP3 ,t_MP3
integer :: maxSCF_HF,n_diis_HF integer :: maxSCF_HF,n_diis_HF
double precision :: thresh_HF,level_shift double precision :: thresh_HF,level_shift
@ -136,10 +109,6 @@ program QuAcK
logical :: BSE,dBSE,dTDA,evDyn,ppBSE,BSE2 logical :: BSE,dBSE,dTDA,evDyn,ppBSE,BSE2
integer :: nMC,nEq,nWalk,nPrint,iSeed
double precision :: dt
logical :: doDrift
! Hello World ! Hello World
write(*,*) write(*,*)
@ -328,7 +297,7 @@ program QuAcK
if(doMOM) then if(doMOM) then
call cpu_time(start_MOM) call cpu_time(start_HF)
if(unrestricted) then if(unrestricted) then
@ -341,10 +310,10 @@ program QuAcK
end if end if
call cpu_time(end_MOM) call cpu_time(end_HF)
t_MOM = end_MOM - start_MOM t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM = ',t_MOM,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM = ',t_HF,' seconds'
write(*,*) write(*,*)
end if end if
@ -475,7 +444,7 @@ program QuAcK
if(doMP2) then if(doMP2) then
call cpu_time(start_MP2) call cpu_time(start_MP)
if(unrestricted) then if(unrestricted) then
@ -487,10 +456,10 @@ program QuAcK
end if end if
call cpu_time(end_MP2) call cpu_time(end_MP)
t_MP2 = end_MP2 - start_MP2 t_MP = end_MP - start_MP
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP2,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds'
write(*,*) write(*,*)
end if end if
@ -501,7 +470,7 @@ program QuAcK
if(doMP3) then if(doMP3) then
call cpu_time(start_MP3) call cpu_time(start_MP)
if(unrestricted) then if(unrestricted) then
@ -514,10 +483,10 @@ program QuAcK
end if end if
call cpu_time(end_MP3) call cpu_time(end_MP)
t_MP3 = end_MP3 - start_MP3 t_MP = end_MP - start_MP
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP3 = ',t_MP3,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP3 = ',t_MP,' seconds'
write(*,*) write(*,*)
end if end if
@ -528,12 +497,12 @@ program QuAcK
if(doCCD) then if(doCCD) then
call cpu_time(start_CCD) call cpu_time(start_CC)
call CCD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call CCD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD) call cpu_time(end_CC)
t_CCD = end_CCD - start_CCD t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CCD,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CC,' seconds'
write(*,*) write(*,*)
end if end if
@ -544,13 +513,13 @@ program QuAcK
if(doDCD) then if(doDCD) then
call cpu_time(start_DCD) call cpu_time(start_CC)
call DCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & call DCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, &
ERI_MO,ENuc,ERHF,eHF) ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_DCD) call cpu_time(end_CC)
t_DCD = end_DCD - start_DCD t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for DCD = ',t_DCD,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for DCD = ',t_CC,' seconds'
write(*,*) write(*,*)
end if end if
@ -563,12 +532,12 @@ program QuAcK
if(doCCSD) then if(doCCSD) then
call cpu_time(start_CCSD) call cpu_time(start_CC)
call CCSD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call CCSD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCSD) call cpu_time(end_CC)
t_CCSD = end_CCSD - start_CCSD t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CCSD,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CC,' seconds'
write(*,*) write(*,*)
end if end if
@ -579,12 +548,12 @@ program QuAcK
if(do_drCCD) then if(do_drCCD) then
call cpu_time(start_CCD) call cpu_time(start_CC)
call drCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call drCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD) call cpu_time(end_CC)
t_CCD = end_CCD - start_CCD t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for direct ring CCD = ',t_CCD,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for direct ring CCD = ',t_CC,' seconds'
write(*,*) write(*,*)
end if end if
@ -595,12 +564,12 @@ program QuAcK
if(do_rCCD) then if(do_rCCD) then
call cpu_time(start_CCD) call cpu_time(start_CC)
call rCCD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,eHF) call rCCD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,eHF)
call cpu_time(end_CCD) call cpu_time(end_CC)
t_CCD = end_CCD - start_CCD t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD = ',t_CCD,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD = ',t_CC,' seconds'
write(*,*) write(*,*)
end if end if
@ -611,13 +580,12 @@ program QuAcK
if(do_crCCD) then if(do_crCCD) then
call cpu_time(start_CCD) call cpu_time(start_CC)
call crCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call crCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CC)
call cpu_time(end_CCD) t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CC,' seconds'
t_CCD = end_CCD - start_CCD
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CCD,' seconds'
write(*,*) write(*,*)
end if end if
@ -628,13 +596,13 @@ program QuAcK
if(do_lCCD) then if(do_lCCD) then
call cpu_time(start_CCD) call cpu_time(start_CC)
call lCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & call lCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, &
ERI_MO,ENuc,ERHF,eHF) ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD) call cpu_time(end_CC)
t_CCD = end_CCD - start_CCD t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ladder CCD = ',t_CCD,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ladder CCD = ',t_CC,' seconds'
write(*,*) write(*,*)
end if end if
@ -645,12 +613,12 @@ program QuAcK
if(dopCCD) then if(dopCCD) then
call cpu_time(start_CCD) call cpu_time(start_CC)
call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD) call cpu_time(end_CC)
t_CCD = end_CCD - start_CCD t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pair CCD = ',t_CCD,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pair CCD = ',t_CC,' seconds'
write(*,*) write(*,*)
end if end if
@ -661,7 +629,7 @@ program QuAcK
if(doCIS) then if(doCIS) then
call cpu_time(start_CIS) call cpu_time(start_CI)
if(unrestricted) then if(unrestricted) then
@ -674,10 +642,10 @@ program QuAcK
end if end if
call cpu_time(end_CIS) call cpu_time(end_CI)
t_CIS = end_CIS - start_CIS t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CIS = ',t_CIS,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CIS = ',t_CI,' seconds'
write(*,*) write(*,*)
end if end if
@ -688,12 +656,12 @@ program QuAcK
if(doCID) then if(doCID) then
call cpu_time(start_CID) call cpu_time(start_CI)
call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,F_MO,ERHF) call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,F_MO,ERHF)
call cpu_time(end_CID) call cpu_time(end_CI)
t_CID = end_CID - start_CID t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CID = ',t_CID,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CID = ',t_CI,' seconds'
write(*,*) write(*,*)
end if end if
@ -704,12 +672,12 @@ program QuAcK
if(doCISD) then if(doCISD) then
call cpu_time(start_CISD) call cpu_time(start_CI)
call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,F_MO,ERHF) call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,F_MO,ERHF)
call cpu_time(end_CISD) call cpu_time(end_CI)
t_CISD = end_CISD - start_CISD t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CISD = ',t_CISD,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CISD = ',t_CI,' seconds'
write(*,*) write(*,*)
end if end if
@ -829,7 +797,7 @@ program QuAcK
if(doG0F2) then if(doG0F2) then
call cpu_time(start_GF2) call cpu_time(start_GF)
if(unrestricted) then if(unrestricted) then
@ -844,10 +812,10 @@ program QuAcK
end if end if
call cpu_time(end_GF2) call cpu_time(end_GF)
t_GF2 = end_GF2 - start_GF2 t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF2,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds'
write(*,*) write(*,*)
end if end if
@ -858,7 +826,7 @@ program QuAcK
if(doevGF2) then if(doevGF2) then
call cpu_time(start_GF2) call cpu_time(start_GF)
if(unrestricted) then if(unrestricted) then
@ -874,10 +842,10 @@ program QuAcK
end if end if
call cpu_time(end_GF2) call cpu_time(end_GF)
t_GF2 = end_GF2 - start_GF2 t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF2,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds'
write(*,*) write(*,*)
end if end if
@ -888,7 +856,7 @@ program QuAcK
if(doqsGF2) then if(doqsGF2) then
call cpu_time(start_GF2) call cpu_time(start_GF)
if(unrestricted) then if(unrestricted) then
@ -903,10 +871,10 @@ program QuAcK
end if end if
call cpu_time(end_GF2) call cpu_time(end_GF)
t_GF2 = end_GF2 - start_GF2 t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGF2 = ',t_GF2,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGF2 = ',t_GF,' seconds'
write(*,*) write(*,*)
end if end if
@ -917,7 +885,7 @@ program QuAcK
if(doG0F3) then if(doG0F3) then
call cpu_time(start_GF3) call cpu_time(start_GF)
if(unrestricted) then if(unrestricted) then
@ -929,10 +897,10 @@ program QuAcK
end if end if
call cpu_time(end_GF3) call cpu_time(end_GF)
t_GF3 = end_GF3 - start_GF3 t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF3,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds'
write(*,*) write(*,*)
end if end if
@ -943,7 +911,7 @@ program QuAcK
if(doevGF3) then if(doevGF3) then
call cpu_time(start_GF3) call cpu_time(start_GF)
if(unrestricted) then if(unrestricted) then
@ -955,10 +923,10 @@ program QuAcK
end if end if
call cpu_time(end_GF3) call cpu_time(end_GF)
t_GF3 = end_GF3 - start_GF3 t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF3,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds'
write(*,*) write(*,*)
end if end if
@ -971,7 +939,7 @@ program QuAcK
if(doG0W0) then if(doG0W0) then
call cpu_time(start_G0W0) call cpu_time(start_GW)
if(unrestricted) then if(unrestricted) then
call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, &
@ -987,16 +955,14 @@ program QuAcK
else else
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, & call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, &
linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
! call ehTM(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, &
! linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
end if end if
end if end if
call cpu_time(end_G0W0) call cpu_time(end_GW)
t_G0W0 = end_G0W0 - start_G0W0 t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_G0W0,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_GW,' seconds'
write(*,*) write(*,*)
end if end if
@ -1007,7 +973,7 @@ program QuAcK
if(doevGW) then if(doevGW) then
call cpu_time(start_evGW) call cpu_time(start_GW)
if(unrestricted) then if(unrestricted) then
call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
@ -1021,10 +987,10 @@ program QuAcK
BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet,linGW,eta_GW,regGW, & BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet,linGW,eta_GW,regGW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
end if end if
call cpu_time(end_evGW) call cpu_time(end_GW)
t_evGW = end_evGW - start_evGW t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGW = ',t_evGW,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGW = ',t_GW,' seconds'
write(*,*) write(*,*)
end if end if
@ -1035,7 +1001,7 @@ program QuAcK
if(doqsGW) then if(doqsGW) then
call wall_time(start_qsGW) call wall_time(start_GW)
if(unrestricted) then if(unrestricted) then
@ -1052,10 +1018,10 @@ program QuAcK
end if end if
call wall_time(end_qsGW) call wall_time(end_GW)
t_qsGW = end_qsGW - start_qsGW t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_qsGW,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds'
write(*,*) write(*,*)
end if end if
@ -1066,7 +1032,7 @@ program QuAcK
if(doSRGqsGW) then if(doSRGqsGW) then
call wall_time(start_qsGW) call wall_time(start_GW)
if(unrestricted) then if(unrestricted) then
@ -1080,10 +1046,10 @@ program QuAcK
end if end if
call wall_time(end_qsGW) call wall_time(end_GW)
t_qsGW = end_qsGW - start_qsGW t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_qsGW,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds'
write(*,*) write(*,*)
end if end if
@ -1094,12 +1060,12 @@ program QuAcK
if(doufG0W0) then if(doufG0W0) then
call cpu_time(start_ufGW) call cpu_time(start_GW)
call ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,TDA_W) call ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,TDA_W)
call cpu_time(end_ufGW) call cpu_time(end_GW)
t_ufGW = end_ufGW - start_ufGW t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0W0 = ',t_ufGW,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0W0 = ',t_GW,' seconds'
write(*,*) write(*,*)
if(BSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0) if(BSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0)
@ -1112,13 +1078,13 @@ program QuAcK
if(doufGW) then if(doufGW) then
call cpu_time(start_ufGW) call cpu_time(start_GW)
call ufGW(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call ufGW(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
! call CCGW(maxSCF_CC,thresh_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) ! call CCGW(maxSCF_CC,thresh_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_ufGW) call cpu_time(end_GW)
t_ufGW = end_ufGW - start_ufGW t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufGW = ',t_ufGW,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufGW = ',t_GW,' seconds'
write(*,*) write(*,*)
if(BSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0) if(BSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0)
@ -1133,7 +1099,7 @@ program QuAcK
if(doG0T0pp) then if(doG0T0pp) then
call cpu_time(start_G0T0) call cpu_time(start_GT)
if(unrestricted) then if(unrestricted) then
@ -1152,10 +1118,10 @@ program QuAcK
end if end if
call cpu_time(end_G0T0) call cpu_time(end_GT)
t_G0T0 = end_G0T0 - start_G0T0 t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_GT,' seconds'
write(*,*) write(*,*)
end if end if
@ -1166,7 +1132,7 @@ program QuAcK
if(doevGTpp) then if(doevGTpp) then
call cpu_time(start_evGT) call cpu_time(start_GT)
if(unrestricted) then if(unrestricted) then
@ -1185,10 +1151,10 @@ program QuAcK
end if end if
call cpu_time(end_evGT) call cpu_time(end_GT)
t_evGT = end_evGT - start_evGT t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_evGT,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_GT,' seconds'
write(*,*) write(*,*)
end if end if
@ -1199,7 +1165,7 @@ program QuAcK
if(doqsGTpp) then if(doqsGTpp) then
call cpu_time(start_qsGT) call cpu_time(start_GT)
if(unrestricted) then if(unrestricted) then
@ -1216,10 +1182,10 @@ program QuAcK
end if end if
call cpu_time(end_qsGT) call cpu_time(end_GT)
t_qsGT = end_qsGT - start_qsGT t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGT = ',t_qsGT,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGT = ',t_GT,' seconds'
write(*,*) write(*,*)
end if end if
@ -1232,7 +1198,7 @@ program QuAcK
if(doG0T0eh) then if(doG0T0eh) then
call cpu_time(start_G0T0) call cpu_time(start_GT)
if(unrestricted) then if(unrestricted) then
@ -1245,10 +1211,10 @@ program QuAcK
end if end if
call cpu_time(end_G0T0) call cpu_time(end_GT)
t_G0T0 = end_G0T0 - start_G0T0 t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_GT,' seconds'
write(*,*) write(*,*)
end if end if
@ -1259,7 +1225,7 @@ program QuAcK
if(doevGTeh) then if(doevGTeh) then
call cpu_time(start_evGT) call cpu_time(start_GT)
if(unrestricted) then if(unrestricted) then
else else
@ -1268,10 +1234,36 @@ program QuAcK
BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet,linGT,eta_GT,regGT, & BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet,linGT,eta_GT,regGT, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0T0) nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0T0)
end if end if
call cpu_time(end_evGT) call cpu_time(end_GT)
t_evGT = end_evGT - start_evGT t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_evGT,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_GT,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform qsGTeh calculation
!------------------------------------------------------------------------
if(doqsGTeh) then
call wall_time(start_GT)
if(unrestricted) then
else
call qsGTeh(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, &
BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GT,regGT,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
end if
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GT,' seconds'
write(*,*) write(*,*)
end if end if
@ -1282,13 +1274,13 @@ program QuAcK
if(doFCI) then if(doFCI) then
call cpu_time(start_FCI) call cpu_time(start_CI)
write(*,*) ' FCI is not yet implemented! Sorry.' write(*,*) ' FCI is not yet implemented! Sorry.'
! call FCI(nBas,nC,nO,nV,nR,ERI_MO,eHF) ! call FCI(nBas,nC,nO,nV,nR,ERI_MO,eHF)
call cpu_time(end_FCI) call cpu_time(end_CI)
t_FCI = end_FCI - start_FCI t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for FCI = ',t_FCI,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for FCI = ',t_CI,' seconds'
write(*,*) write(*,*)
end if end if