PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00
Code Releases Activity

introduce SRG regularizer in GF, GW, and GT and add option for reg MP2

Browse Source
This commit is contained in:
Pierre-Francois Loos 2022-11-30 17:18:51 +01:00
parent 597c4a9885
commit 196e4c917a
13 changed files with 92 additions and 87 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -1,7 +1,7 @@
# RHF UHF KS MOM
T F F F
# MP2* MP3 MP2-F12
F F F
T F F
# CCD pCCD DCD CCSD CCSD(T)
F F F F F
# drCCD rCCD crCCD lCCD

4
input/options
View File

@ -1,7 +1,7 @@
# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability
512 0.0000001 T 5 1 1 F 0.0 F
# MP:
# MP: reg
F
# CC: maxSCF thresh DIIS n_diis
64 0.0000001 T 5
# spin: TDA singlet triplet spin_conserved spin_flip

12
src/GF/regularized_self_energy_GF2.f90
View File

@ -37,7 +37,7 @@ subroutine regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
kappa = 1d0
!----------------------------------------------------!
! Compute GF2 self-energy and renormalization factor !
@ -52,9 +52,8 @@ subroutine regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC
eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(q,a,i,j)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2)
SigC(p,q) = SigC(p,q) + num*fk
if(p == q) Z(p) = Z(p) - num*dfk
@ -74,9 +73,8 @@ subroutine regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC
eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(q,i,a,b)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2)
SigC(p,q) = SigC(p,q) + num*fk
if(p == q) Z(p) = Z(p) - num*dfk

12
src/GF/regularized_self_energy_GF2_diag.f90
View File

@ -37,7 +37,7 @@ subroutine regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
kappa = 1d0
!----------------------------------------------------!
! Compute GF2 self-energy and renormalization factor !
@ -51,9 +51,8 @@ subroutine regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI
eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2)
SigC(p) = SigC(p) + num*fk
Z(p) = Z(p) - num*dfk
@ -71,9 +70,8 @@ subroutine regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI
eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2)
SigC(p) = SigC(p) + num*fk
Z(p) = Z(p) - num*dfk

13
src/GT/regularized_renormalization_factor_Tmatrix.f90
View File

@ -33,7 +33,7 @@ subroutine regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,n
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
kappa = 1d0
! Occupied part of the T-matrix self-energy
@ -43,10 +43,8 @@ subroutine regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,n
eps = e(p) + e(i) - Omega1(cd)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2)
Z(p) = Z(p) - rho1(p,i,cd)**2*dfk
enddo
@ -61,9 +59,8 @@ subroutine regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,n
eps = e(p) + e(nO+a) - Omega2(kl)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2)
Z(p) = Z(p) - rho2(p,nO+a,kl)**2*dfk

10
src/GT/regularized_self_energy_Tmatrix.f90
View File

@ -37,7 +37,7 @@ subroutine regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
kappa = 1d0
!----------------------------------------------
! Occupied part of the T-matrix self-energy
@ -49,9 +49,9 @@ subroutine regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1
do cd=1,nVV
eps = e(p) + e(i) - Omega1(cd)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
SigT(p,q) = SigT(p,q) + rho1(p,i,cd)*rho1(q,i,cd)*eps/(eps**2 + eta**2)
SigT(p,q) = SigT(p,q) + rho1(p,i,cd)*rho1(q,i,cd)*fk
enddo
enddo
@ -68,9 +68,9 @@ subroutine regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1
do kl=1,nOO
eps = e(p) + e(a) - Omega2(kl)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
SigT(p,q) = SigT(p,q) + rho2(p,a,kl)*rho2(q,a,kl)*eps/(eps**2 + eta**2)
SigT(p,q) = SigT(p,q) + rho2(p,a,kl)*rho2(q,a,kl)*fk
enddo
enddo

6
src/GT/regularized_self_energy_Tmatrix_diag.f90
View File

@ -37,7 +37,7 @@ subroutine regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,O
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
kappa = 1d0
!----------------------------------------------
! Occupied part of the T-matrix self-energy
@ -48,7 +48,7 @@ subroutine regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,O
do cd=1,nVV
eps = e(p) + e(i) - Omega1(cd)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
SigT(p) = SigT(p) + rho1(p,i,cd)**2*fk
@ -65,7 +65,7 @@ subroutine regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,O
do kl=1,nOO
eps = e(p) + e(a) - Omega2(kl)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
SigT(p) = SigT(p) + rho2(p,a,kl)**2*fk

12
src/GW/regularized_renormalization_factor.f90
View File

@ -39,7 +39,7 @@ subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
kappa = 1d0
! static COHSEX approximation
@ -56,9 +56,8 @@ subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
do i=nC+1,nO
do jb=1,nS
eps = e(p) - e(i) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2)
Z(p) = Z(p) - 2d0*rho(p,i,jb)**2*dfk
end do
end do
@ -70,9 +69,8 @@ subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,O
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(p) - e(a) - Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps)))
dfk = dfk*fk
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2)
Z(p) = Z(p) - 2d0*rho(p,a,jb)**2*dfk
end do
end do

8
src/GW/regularized_self_energy_correlation.f90
View File

@ -42,7 +42,7 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,
! Parameters for regularized MP2 calculations !
!---------------------------------------------!
kappa = 1.1d0
kappa = 1d0
!-----------------------------!
! COHSEX static approximation !
@ -92,7 +92,7 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,
do i=nC+1,nO
do jb=1,nS
eps = e(p) - e(i) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,jb)*rho(q,i,jb)*fk
end do
end do
@ -106,7 +106,7 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(p) - e(a) - Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,jb)*rho(q,a,jb)*fk
end do
end do
@ -120,7 +120,7 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(a) - e(i) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
EcGM = EcGM - 4d0*rho(a,i,jb)**2*fk
end do
end do

8
src/GW/regularized_self_energy_correlation_diag.f90
View File

@ -40,7 +40,7 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,
! Parameters for regularized calculations !
!-----------------------------------------!
kappa = 1.1d0
kappa = 1d0
!-----------------------------
! COHSEX static self-energy
@ -87,7 +87,7 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,
do i=nC+1,nO
do jb=1,nS
eps = e(p) - e(i) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2*fk
end do
end do
@ -99,7 +99,7 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(p) - e(a) - Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
SigC(p) = SigC(p) + 2d0*rho(p,a,jb)**2*fk
end do
end do
@ -112,7 +112,7 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,
do a=nO+1,nBas-nR
do jb=1,nS
eps = e(a) - e(i) + Omega(jb)
fk = (1d0 - exp(-kappa*abs(eps)))**2/eps
fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps
EcGM = EcGM - 4d0*rho(a,i,jb)**2*fk
end do
end do

79
src/MP/MP2.f90
View File

@ -1,4 +1,4 @@
subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
subroutine MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
! Perform second-order Moller-Plesset calculation with and without regularizers
@ -6,6 +6,7 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
! Input variables
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -115,55 +116,59 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
write(*,'(A32)') '--------------------------'
write(*,*)
if(regularize) then
!-------------------!
! sigma1-MP2 energy !
!-------------------!
write(*,*)
write(*,'(A32)') '--------------------------'
write(*,'(A32)') ' sigma-MP2 calculation '
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EcsMP2
write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2ds
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xs
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EcsMP2
write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EcsMP2
write(*,'(A32)') '--------------------------'
write(*,*)
write(*,*)
write(*,'(A32)') '--------------------------'
write(*,'(A32)') ' sigma-MP2 calculation '
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EcsMP2
write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2ds
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xs
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EcsMP2
write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EcsMP2
write(*,'(A32)') '--------------------------'
write(*,*)
!--------------------!
! sigma^2-MP2 energy !
!--------------------!
write(*,*)
write(*,'(A32)') '--------------------------'
write(*,'(A32)') ' sigma^2-MP2 calculation '
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',Ecs2MP2
write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2ds2
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xs2
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + Ecs2MP2
write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + Ecs2MP2
write(*,'(A32)') '--------------------------'
write(*,*)
write(*,*)
write(*,'(A32)') '--------------------------'
write(*,'(A32)') ' sigma^2-MP2 calculation '
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',Ecs2MP2
write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2ds2
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xs2
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + Ecs2MP2
write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + Ecs2MP2
write(*,'(A32)') '--------------------------'
write(*,*)
!------------------!
! kappa-MP2 energy !
!------------------!
write(*,*)
write(*,'(A32)') '--------------------------'
write(*,'(A32)') ' kappa-MP2 calculation '
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EckMP2
write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2dk
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xk
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EckMP2
write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EckMP2
write(*,'(A32)') '--------------------------'
write(*,*)
write(*,*)
write(*,'(A32)') '--------------------------'
write(*,'(A32)') ' kappa-MP2 calculation '
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EckMP2
write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2dk
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xk
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EckMP2
write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EckMP2
write(*,'(A32)') '--------------------------'
write(*,*)
end if
end subroutine MP2

5
src/QuAcK/QuAcK.f90
View File

@ -108,6 +108,8 @@ program QuAcK
double precision :: thresh_HF,level_shift
logical :: DIIS_HF,guess_type,ortho_type,mix
logical :: regMP
integer :: maxSCF_CC,n_diis_CC
double precision :: thresh_CC
logical :: DIIS_CC
@ -174,6 +176,7 @@ program QuAcK
! Read options for methods
call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,level_shift,dostab, &
regMP, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
TDA,singlet,triplet,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, &
@ -479,7 +482,7 @@ program QuAcK
else
call MP2(nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,EcMP2)
call MP2(regMP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,EcMP2)
end if

8
src/QuAcK/read_options.f90
View File

@ -1,4 +1,5 @@
subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,level_shift,dostab, &
regMP, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
TDA,singlet,triplet,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, &
@ -24,6 +25,8 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
double precision,intent(out) :: level_shift
logical,intent(out) :: dostab
logical,intent(out) :: regMP
integer,intent(out) :: maxSCF_CC
double precision,intent(out) :: thresh_CC
logical,intent(out) :: DIIS_CC
@ -106,8 +109,11 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
! Read MPn options
read(1,*)
regMP = .false.
read(1,*)
read(1,*) answer1
if(answer1 == 'T') regMP = .true.
! Read CC options