From e006fb5a01e4fb3850997af242e1268df06f6bb9 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 16 Dec 2021 14:24:02 +0100 Subject: [PATCH] regularized MP2 --- input/methods | 4 +- src/MP/MP2.f90 | 130 ++++++++++++++++++++++++++++++++++++++++++------- 2 files changed, 114 insertions(+), 20 deletions(-) diff --git a/input/methods b/input/methods index 77c21ce..a1580bd 100644 --- a/input/methods +++ b/input/methods @@ -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 diff --git a/src/MP/MP2.f90 b/src/MP/MP2.f90 index 7dd921f..cfca705 100644 --- a/src/MP/MP2.f90 +++ b/src/MP/MP2.f90 @@ -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