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

regularized MP2

Browse Source
This commit is contained in:
Pierre-Francois Loos 2021-12-16 14:24:02 +01:00
parent 5d5824ce14
commit e006fb5a01
2 changed files with 114 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
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
@ -13,7 +13,7 @@
# G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F
# G0W0* evGW* qsGW* ufG0W0 ufGW
F F T F F
F F F F F
# G0T0 evGT qsGT
F F F
# MCMP2

130
src/MP/MP2.f90
View File

@ -1,6 +1,6 @@
subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
! Perform second-order Moller-Plesset calculation
! Perform second-order Moller-Plesset calculation with and without regularizers
implicit none
@ -19,7 +19,15 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
! Local variables
integer :: i,j,a,b
double precision :: eps,E2a,E2b
double precision :: kappa,sigm1,sigm2
double precision :: Dijab
double precision :: fs,fs2,fk
double precision :: E2d,E2ds,E2ds2,E2dk
double precision :: E2x,E2xs,E2xs2,E2xk
double precision :: EcsMP2,Ecs2MP2,EckMP2
! Output variables
@ -33,43 +41,129 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
write(*,*)'************************************************'
write(*,*)
! Compute MP2 energy
!---------------------------------------------!
! Parameters for regularized MP2 calculations !
!---------------------------------------------!
kappa = 1.1d0
sigm1 = 0.7d0
sigm2 = 0.4d0
!--------------------------------------------------!
! Compute conventinal and regularized MP2 energies !
!--------------------------------------------------!
E2d = 0d0
E2ds = 0d0
E2ds2 = 0d0
E2dk = 0d0
E2x = 0d0
E2xs = 0d0
E2xs2 = 0d0
E2xk = 0d0
E2a = 0d0
E2b = 0d0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
do j=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = e(i) + e(j) - e(a) - e(b)
Dijab = e(a) + e(b) - e(i) - e(j)
! Second-order ring diagram
! Second-order ring diagram
E2a = E2a + ERI(i,j,a,b)*ERI(i,j,a,b)/eps
fs = (1d0 - exp(-kappa*Dijab))/Dijab
fs2 = (1d0 - exp(-sigm1*Dijab*Dijab))/Dijab
fk = (1d0 - exp(-sigm2*Dijab))**2/Dijab
! Second-order exchange diagram
E2d = E2d - ERI(i,j,a,b)*ERI(i,j,a,b)/Dijab
E2ds = E2ds - ERI(i,j,a,b)*ERI(i,j,a,b)*fs
E2ds2 = E2ds2 - ERI(i,j,a,b)*ERI(i,j,a,b)*fs2
E2dk = E2dk - ERI(i,j,a,b)*ERI(i,j,a,b)*fk
E2b = E2b + ERI(i,j,a,b)*ERI(i,j,b,a)/eps
! Second-order exchange diagram
enddo
E2x = E2x - ERI(i,j,a,b)*ERI(i,j,b,a)/Dijab
E2xs = E2xs - ERI(i,j,a,b)*ERI(i,j,b,a)*fs
E2xs2 = E2xs2 - ERI(i,j,a,b)*ERI(i,j,b,a)*fs2
E2xk = E2xk - ERI(i,j,a,b)*ERI(i,j,b,a)*fk
enddo
enddo
enddo
enddo
enddo
EcMP2 = 2d0*E2a - E2b
EcMP2 = 2d0*E2d - E2x
EcsMP2 = 2d0*E2ds - E2xs
Ecs2MP2 = 2d0*E2ds2 - E2xs2
EckMP2 = 2d0*E2dk - E2xk
!------------!
! MP2 energy !
!------------!
write(*,*)
write(*,'(A32)') '--------------------------'
write(*,'(A32)') ' MP2 calculation '
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EcMP2
write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2a
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2b
write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2d
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2x
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EcMP2
write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EcMP2
write(*,'(A32)') '--------------------------'
write(*,*)
!-------------------!
! 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(*,*)
!--------------------!
! 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(*,*)
!------------------!
! 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(*,*)
end subroutine MP2