10
1
mirror of https://github.com/pfloos/quack synced 2024-10-20 06:48:15 +02:00

remove RMP3

This commit is contained in:
Pierre-Francois Loos 2024-09-10 17:08:16 +02:00
parent 66c59f85a1
commit 34fba9ed9c
2 changed files with 4 additions and 197 deletions

View File

@ -31,7 +31,7 @@ subroutine RMP(dotest,doMP2,doMP3,regularize,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
! Output variables ! Output variables
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Compute MP3 energy ! Compute MP2 energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------
if(doMP2) then if(doMP2) then
@ -53,11 +53,12 @@ subroutine RMP(dotest,doMP2,doMP3,regularize,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
if(doMP3) then if(doMP3) then
call wall_time(start_MP) call wall_time(start_MP)
call RMP3(nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) write(*,*) 'Restricted MP3 NYI... Sorry'
write(*,*)
call wall_time(end_MP) call wall_time(end_MP)
t_MP = end_MP - start_MP t_MP = end_MP - start_MP
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for MP2 = ',t_MP,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for MP3 = ',t_MP,' seconds'
write(*,*) write(*,*)
end if end if

View File

@ -1,194 +0,0 @@
subroutine RMP3(nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,EHF,e)
! Perform third-order Moller-Plesset calculation
implicit none
! Input variables
integer,intent(in) :: nBasin
integer,intent(in) :: nCin
integer,intent(in) :: nOin
integer,intent(in) :: nVin
integer,intent(in) :: nRin
double precision,intent(in) :: ENuc,EHF
double precision,intent(in) :: e(nBasin)
double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin)
! Local variables
double precision :: eps,E2,EcMP2
double precision :: eps1,eps2,E3a,E3b,E3c
double precision :: EcMP3
integer :: nBas
integer :: nC
integer :: nO
integer :: nV
integer :: nR
integer :: i,j,k,l,a,b,c,d
double precision,allocatable :: se(:)
double precision,allocatable :: eO(:)
double precision,allocatable :: eV(:)
double precision,allocatable :: sERI(:,:,:,:)
double precision,allocatable :: dbERI(:,:,:,:)
double precision,allocatable :: OOOO(:,:,:,:)
double precision,allocatable :: OOVV(:,:,:,:)
double precision,allocatable :: OVVO(:,:,:,:)
double precision,allocatable :: VVOO(:,:,:,:)
double precision,allocatable :: VVVV(:,:,:,:)
! Hello world
write(*,*)
write(*,*)'******************************'
write(*,*)'* Restricted MP3 Calculation *'
write(*,*)'******************************'
write(*,*)
! Spatial to spin orbitals
nBas = 2*nBasin
nC = 2*nCin
nO = 2*nOin
nV = 2*nVin
nR = 2*nRin
allocate(se(nBas),sERI(nBas,nBas,nBas,nBas))
call spatial_to_spin_MO_energy(nBasin,e,nBas,se)
call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI)
! Antysymmetrize ERIs
allocate(dbERI(nBas,nBas,nBas,nBas))
call antisymmetrize_ERI(2,nBas,sERI,dbERI)
deallocate(sERI)
! Form energy denominator
allocate(eO(nO),eV(nV))
eO(:) = se(1:nO)
eV(:) = se(nO+1:nBas)
deallocate(se)
! Create integral batches
allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VVOO(nV,nV,nO,nO),VVVV(nV,nV,nV,nV))
OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO )
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
VVOO(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas, 1:nO , 1:nO )
VVVV(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas,nO+1:nBas,nO+1:nBas)
deallocate(dbERI)
! Compute MP2 energy
E2 = 0d0
do i=nC+1,nO
do j=nC+1,nO
do a=1,nV-nR
do b=1,nV-nR
eps = eO(i) + eO(j) - eV(a) - eV(b)
E2 = E2 + OOVV(i,j,a,b)*OOVV(i,j,a,b)/eps
end do
end do
end do
end do
EcMP2 = 0.25d0*E2
! Compute MP3 energy
E3a = 0d0
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do l=nC+1,nO
do a=1,nV-nR
do b=1,nV-nR
eps1 = eO(i) + eO(j) - eV(a) - eV(b)
eps2 = eO(k) + eO(l) - eV(a) - eV(b)
E3a = E3a + OOVV(i,j,a,b)*OOOO(k,l,i,j)*VVOO(a,b,k,l)/(eps1*eps2)
end do
end do
end do
end do
end do
end do
E3b = 0d0
do i=nC+1,nO
do j=nC+1,nO
do a=1,nV-nR
do b=1,nV-nR
do c=1,nV-nR
do d=1,nV-nR
eps1 = eO(i) + eO(j) - eV(a) - eV(b)
eps2 = eO(i) + eO(j) - eV(c) - eV(d)
E3b = E3b + OOVV(i,j,a,b)*VVVV(a,b,c,d)*VVOO(c,d,i,j)/(eps1*eps2)
end do
end do
end do
end do
end do
end do
E3c = 0d0
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=1,nV-nR
do b=1,nV-nR
do c=1,nV-nR
eps1 = eO(i) + eO(j) - eV(a) - eV(b)
eps2 = eO(i) + eO(k) - eV(a) - eV(c)
E3c = E3c + OOVV(i,j,a,b)*OVVO(k,b,c,j)*VVOO(a,c,i,k)/(eps1*eps2)
end do
end do
end do
end do
end do
end do
EcMP3 = 0.125d0*E3a + 0.125d0*E3b + E3c
write(*,*)
write(*,'(A32)') '-----------------------'
write(*,'(A32)') ' MP3 calculation '
write(*,'(A32)') '-----------------------'
write(*,'(A32,1X,F16.10)') ' MP2 contribution ',EcMP2
write(*,'(A32,1X,F16.10)') ' MP3 contribution ',EcMP3
write(*,'(A32)') '-----------------------'
write(*,'(A32,1X,F16.10)') ' MP3 correlation energy', EcMP2 + EcMP3
write(*,'(A32,1X,F16.10)') ' MP3 total energy',ENuc + EHF + EcMP2 + EcMP3
write(*,'(A32)') '-----------------------'
write(*,*)
end subroutine