mirror of https://github.com/pfloos/quack
187 lines
5.7 KiB
Fortran
187 lines
5.7 KiB
Fortran
subroutine GMP2(dotest,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EGHF,eHF,EcMP2)
|
|
|
|
! Perform second-order Moller-Plesset calculation with and without regularizers
|
|
|
|
implicit none
|
|
|
|
! Input variables
|
|
|
|
logical,intent(in) :: dotest
|
|
|
|
logical,intent(in) :: regularize
|
|
integer,intent(in) :: nBas
|
|
integer,intent(in) :: nC
|
|
integer,intent(in) :: nO
|
|
integer,intent(in) :: nV
|
|
integer,intent(in) :: nR
|
|
double precision,intent(in) :: ENuc
|
|
double precision,intent(in) :: EGHF
|
|
double precision,intent(in) :: eHF(nBas)
|
|
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
|
|
|
! Local variables
|
|
|
|
integer :: i,j,a,b
|
|
|
|
double precision :: kappa,sigm1,sigm2
|
|
double precision :: num
|
|
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
|
|
|
|
double precision,intent(out) :: EcMP2
|
|
|
|
! Hello world
|
|
|
|
write(*,*)
|
|
write(*,*)'*******************************'
|
|
write(*,*)'* Generalized MP2 Calculation |'
|
|
write(*,*)'*******************************'
|
|
write(*,*)
|
|
|
|
!---------------------------------------------!
|
|
! 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
|
|
|
|
do i=nC+1,nO
|
|
do j=nC+1,nO
|
|
do a=nO+1,nBas-nR
|
|
do b=nO+1,nBas-nR
|
|
|
|
Dijab = eHF(a) + eHF(b) - eHF(i) - eHF(j)
|
|
|
|
! Second-order ring diagram
|
|
|
|
fs = (1d0 - exp(-sigm1*Dijab))/Dijab
|
|
fs2 = (1d0 - exp(-sigm2*Dijab*Dijab))/Dijab
|
|
fk = (1d0 - exp(-kappa*Dijab))**2/Dijab
|
|
|
|
num = 0.25d0*(ERI(i,j,a,b)**2 + ERI(i,j,b,a)**2)
|
|
E2d = E2d - num/Dijab
|
|
E2ds = E2ds - num*fs
|
|
E2ds2 = E2ds2 - num*fs2
|
|
E2dk = E2dk - num*fk
|
|
|
|
! Second-order exchange diagram
|
|
|
|
num = 0.5d0*ERI(i,j,a,b)*ERI(i,j,b,a)
|
|
E2x = E2x + num/Dijab
|
|
E2xs = E2xs + num*fs
|
|
E2xs2 = E2xs2 + num*fs2
|
|
E2xk = E2xk + num*fk
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
EcMP2 = E2d + E2x
|
|
EcsMP2 = E2ds + E2xs
|
|
Ecs2MP2 = E2ds2 + E2xs2
|
|
EckMP2 = E2dk + E2xk
|
|
|
|
!------------!
|
|
! MP2 energy !
|
|
!------------!
|
|
|
|
write(*,*)
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,'(A32)') ' GMP2 calculation '
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,'(A32,1X,F16.10)') ' GMP2 correlation energy = ',EcMP2
|
|
write(*,'(A32,1X,F16.10)') ' Direct part = ',E2d
|
|
write(*,'(A32,1X,F16.10)') ' Exchange part = ',E2x
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EGHF + EcMP2
|
|
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EGHF + EcMP2
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,*)
|
|
|
|
if(regularize) then
|
|
|
|
!-------------------!
|
|
! sigma1-MP2 energy !
|
|
!-------------------!
|
|
|
|
write(*,*)
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,'(A32)') ' sigma-GMP2 calculation '
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,'(A32,1X,F16.10)') ' GMP2 correlation energy = ',EcsMP2
|
|
write(*,'(A32,1X,F16.10)') ' Direct part = ',E2ds
|
|
write(*,'(A32,1X,F16.10)') ' Exchange part = ',E2xs
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EGHF + EcsMP2
|
|
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EGHF + EcsMP2
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,*)
|
|
|
|
!--------------------!
|
|
! sigma^2-MP2 energy !
|
|
!--------------------!
|
|
|
|
write(*,*)
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,'(A32)') ' sigma^2-GMP2 calculation '
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,'(A32,1X,F16.10)') ' GMP2 correlation energy = ',Ecs2MP2
|
|
write(*,'(A32,1X,F16.10)') ' Direct part = ',E2ds2
|
|
write(*,'(A32,1X,F16.10)') ' Exchange part = ',E2xs2
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EGHF + Ecs2MP2
|
|
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EGHF + Ecs2MP2
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,*)
|
|
|
|
!------------------!
|
|
! kappa-MP2 energy !
|
|
!------------------!
|
|
|
|
write(*,*)
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,'(A32)') ' kappa-GMP2 calculation '
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,'(A32,1X,F16.10)') ' GMP2 correlation energy = ',EckMP2
|
|
write(*,'(A32,1X,F16.10)') ' Direct part = ',E2dk
|
|
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xk
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EGHF + EckMP2
|
|
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EGHF + EckMP2
|
|
write(*,'(A32)') '---------------------------'
|
|
write(*,*)
|
|
|
|
end if
|
|
|
|
if(dotest) then
|
|
|
|
call dump_test_value('G','MP2 correlation energy',EcMP2)
|
|
|
|
end if
|
|
|
|
|
|
end subroutine
|