From 154b395af2c564b63a357edd9b2ba7ecb9782cc6 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sun, 17 Mar 2019 10:27:06 +0100 Subject: [PATCH] MP3 is working --- input/basis | 49 +++++------- input/methods | 4 +- input/molecule | 4 +- input/weight | 49 +++++------- src/MCQC/MCQC.f90 | 4 +- src/MCQC/MP2F12.f90 | 45 +++++++++-- src/MCQC/MP3.f90 | 182 +++++++++++++++++++++++++++++++------------- 7 files changed, 213 insertions(+), 124 deletions(-) diff --git a/input/basis b/input/basis index eb58543..c12e4d6 100644 --- a/input/basis +++ b/input/basis @@ -1,29 +1,20 @@ -1 6 -S 8 1.00 - 17880.0000000 0.0007380 - 2683.0000000 0.0056770 - 611.5000000 0.0288830 - 173.5000000 0.1085400 - 56.6400000 0.2909070 - 20.4200000 0.4483240 - 7.8100000 0.2580260 - 1.6530000 0.0150630 -S 8 1.00 - 17880.0000000 -0.0001720 - 2683.0000000 -0.0013570 - 611.5000000 -0.0067370 - 173.5000000 -0.0276630 - 56.6400000 -0.0762080 - 20.4200000 -0.1752270 - 7.8100000 -0.1070380 - 1.6530000 0.5670500 -S 1 1.00 - 0.4869000 1.0000000 -P 3 1.00 - 28.3900000 0.0460870 - 6.2700000 0.2401810 - 1.6950000 0.5087440 -P 1 1.00 - 0.4317000 1.0000000 -D 1 1.00 - 2.2020000 1.0000000 +1 5 +S 6 1.00 + 1264.5857000 0.0019448 + 189.9368100 0.0148351 + 43.1590890 0.0720906 + 12.0986630 0.2371542 + 3.8063232 0.4691987 + 1.2728903 0.3565202 +S 3 1.00 + 3.1964631 -0.1126487 + 0.7478133 -0.2295064 + 0.2199663 1.1869167 +S 1 1.00 + 0.0823099 1.0000000 +P 3 1.00 + 3.1964631 0.0559802 + 0.7478133 0.2615506 + 0.2199663 0.7939723 +P 1 1.00 + 0.0823099 1.0000000 diff --git a/input/methods b/input/methods index ff4ad8e..a3433d1 100644 --- a/input/methods +++ b/input/methods @@ -1,9 +1,9 @@ # HF MOM T F # MP2 MP3 MP2-F12 - T F F + T T F # CCD CCSD CCSD(T) - F T T + F F F # CIS TDHF ADC F F F # GF2 GF3 diff --git a/input/molecule b/input/molecule index a0a7bc2..67ee9ac 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEl nCore nRyd - 1 10 0 0 + 1 4 0 0 # Znuc x y z - Ne 0.0 0.0 0.0 + Be 0.0 0.0 0.0 diff --git a/input/weight b/input/weight index eb58543..c12e4d6 100644 --- a/input/weight +++ b/input/weight @@ -1,29 +1,20 @@ -1 6 -S 8 1.00 - 17880.0000000 0.0007380 - 2683.0000000 0.0056770 - 611.5000000 0.0288830 - 173.5000000 0.1085400 - 56.6400000 0.2909070 - 20.4200000 0.4483240 - 7.8100000 0.2580260 - 1.6530000 0.0150630 -S 8 1.00 - 17880.0000000 -0.0001720 - 2683.0000000 -0.0013570 - 611.5000000 -0.0067370 - 173.5000000 -0.0276630 - 56.6400000 -0.0762080 - 20.4200000 -0.1752270 - 7.8100000 -0.1070380 - 1.6530000 0.5670500 -S 1 1.00 - 0.4869000 1.0000000 -P 3 1.00 - 28.3900000 0.0460870 - 6.2700000 0.2401810 - 1.6950000 0.5087440 -P 1 1.00 - 0.4317000 1.0000000 -D 1 1.00 - 2.2020000 1.0000000 +1 5 +S 6 1.00 + 1264.5857000 0.0019448 + 189.9368100 0.0148351 + 43.1590890 0.0720906 + 12.0986630 0.2371542 + 3.8063232 0.4691987 + 1.2728903 0.3565202 +S 3 1.00 + 3.1964631 -0.1126487 + 0.7478133 -0.2295064 + 0.2199663 1.1869167 +S 1 1.00 + 0.0823099 1.0000000 +P 3 1.00 + 3.1964631 0.0559802 + 0.7478133 0.2615506 + 0.2199663 0.7939723 +P 1 1.00 + 0.0823099 1.0000000 diff --git a/src/MCQC/MCQC.f90 b/src/MCQC/MCQC.f90 index 9a2386a..35f2538 100644 --- a/src/MCQC/MCQC.f90 +++ b/src/MCQC/MCQC.f90 @@ -231,7 +231,7 @@ program MCQC if(doMP3) then call cpu_time(start_MP3) - call MP3(nBas,nC,nO,nV,nR,ERI_MO_basis,ERHF,EcMP2(1),eHF,EcMP3) + call MP3(nBas,nEl,ERI_MO_basis,eHF,ENuc,ERHF) call cpu_time(end_MP3) t_MP3 = end_MP3 - start_MP3 @@ -251,7 +251,7 @@ program MCQC allocate(F12(nBas,nBas,nBas,nBas),Yuk(nBas,nBas,nBas,nBas),FC(nBas,nBas,nBas,nBas,nBas,nBas)) ! Read integrals call read_F12_integrals(nBas,S,ERI_AO_basis,F12,Yuk,FC) - call MP2F12(nBas,nC,nO,nV,ERI_AO_basis,F12,Yuk,FC,ERHF,EcMP2(1),cHF,EcMP2F12) + call MP2F12(nBas,nC,nO,nV,ERI_AO_basis,F12,Yuk,FC,ERHF,eHF,cHF) call cpu_time(end_MP2F12) t_MP2F12 = end_MP2F12 - start_MP2F12 diff --git a/src/MCQC/MP2F12.f90 b/src/MCQC/MP2F12.f90 index f4f9efa..a549e6a 100644 --- a/src/MCQC/MP2F12.f90 +++ b/src/MCQC/MP2F12.f90 @@ -1,4 +1,4 @@ -subroutine MP2F12(nBas,nC,nO,nV,ERI,F12,Yuk,FC,EHF,EcMP2,c,EcMP2F12) +subroutine MP2F12(nBas,nC,nO,nV,ERI,F12,Yuk,FC,EHF,e,c) ! Perform MP2-F12 calculation @@ -7,7 +7,8 @@ subroutine MP2F12(nBas,nC,nO,nV,ERI,F12,Yuk,FC,EHF,EcMP2,c,EcMP2F12) ! Input variables integer,intent(in) :: nBas,nC,nO,nV - double precision,intent(in) :: EHF,EcMP2 + double precision,intent(in) :: EHF + double precision,intent(in) :: e(nBas) double precision,intent(in) :: c(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: F12(nBas,nBas,nBas,nBas) @@ -22,20 +23,52 @@ subroutine MP2F12(nBas,nC,nO,nV,ERI,F12,Yuk,FC,EHF,EcMP2,c,EcMP2F12) double precision,allocatable :: ooCvv(:,:,:,:) double precision,allocatable :: ooFvv(:,:,:,:) double precision,allocatable :: oooFCooo(:,:,:,:,:,:) + double precision,allocatable :: eO(:),eV(:) double precision,allocatable :: cO(:,:),cV(:,:) double precision :: E2a,E2b,E3a,E3b,E4a,E4b,E4c,E4d integer :: i,j,k,l,a,b + double precision :: EcMP2F12(4) -! Output variables - - double precision,intent(out) :: EcMP2F12(3) + double precision :: EcMP2a + double precision :: EcMP2b + double precision :: EcMP2 + double precision :: eps + ! Split MOs into occupied and virtual sets + allocate(eO(nO),eV(nV)) + + eO(1:nO) = e(nC+1:nC+nO) + eV(1:nV) = e(nC+nO+1:nBas) + allocate(cO(nBas,nO),cV(nBas,nV)) + cO(1:nBas,1:nO) = c(1:nBas,nC+1:nC+nO) cV(1:nBas,1:nV) = c(1:nBas,nC+nO+1:nBas) +! Compute conventional MP2 energy + + allocate(ooCvv(nO,nO,nV,nV)) + call AOtoMO_oovv(nBas,nO,nV,cO,cV,ERI,ooCvv) + + EcMP2a = 0d0 + EcMP2b = 0d0 + + do i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + eps = eO(i) + eO(j) - eV(a) - eV(b) + EcMP2a = EcMP2a + ooCoo(i,j,a,b)/eps + EcMP2b = EcMP2b + ooCoo(i,j,b,a)/eps + enddo + enddo + enddo + enddo + + EcMP2 = EcMP2a + EcMP2b + ! Compute the two-electron part of the MP2-F12 energy allocate(ooYoo(nO,nO,nO,nO)) @@ -120,7 +153,7 @@ subroutine MP2F12(nBas,nC,nO,nV,ERI,F12,Yuk,FC,EHF,EcMP2,c,EcMP2F12) write(*,'(A32)') '-----------------------' write(*,'(A32)') ' MP2-F12 calculation ' write(*,'(A32)') '-----------------------' - write(*,'(A32,1X,F16.10)') ' MP2 ',EcMP2 + write(*,'(A32,1X,F16.10)') ' MP2 ',+EcMP2 write(*,'(A32,1X,F16.10)') ' MP2-F12 E(2) ',-EcMP2F12(1) write(*,'(A32,1X,F16.10)') ' MP2-F12 E(3) ',-EcMP2F12(2) write(*,'(A32,1X,F16.10)') ' MP2-F12 E(4) ',-EcMP2F12(3) diff --git a/src/MCQC/MP3.f90 b/src/MCQC/MP3.f90 index c84d7bc..0539d31 100644 --- a/src/MCQC/MP3.f90 +++ b/src/MCQC/MP3.f90 @@ -1,4 +1,4 @@ -subroutine MP3(nBas,nC,nO,nV,nR,V,EHF,EcMP2,e,EcMP3) +subroutine MP3(nBas,nEl,ERI,e,ENuc,EHF) ! Perform third-order Moller-Plesset calculation @@ -6,19 +6,33 @@ subroutine MP3(nBas,nC,nO,nV,nR,V,EHF,EcMP2,e,EcMP3) ! Input variables - integer,intent(in) :: nBas,nC,nO,nV,nR - double precision,intent(in) :: EHF,EcMP2 - double precision,intent(in) :: V(nBas,nBas,nBas,nBas),e(nBas) + integer,intent(in) :: nBas,nEl + double precision,intent(in) :: ENuc,EHF + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables + double precision :: eps,E2,EcMP2 double precision :: eps1,eps2,E3a,E3b,E3c + double precision :: EcMP3 + integer :: nBas2,nO,nV integer :: i,j,k,l,a,b,c,d -! Output variables + 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(:,:,:,:) - double precision,intent(out) :: EcMP3 ! Hello world @@ -28,75 +42,135 @@ subroutine MP3(nBas,nC,nO,nV,nR,V,EHF,EcMP2,e,EcMP3) write(*,*)'************************************************' write(*,*) +! Spatial to spin orbitals + + nBas2 = 2*nBas + + allocate(se(nBas2),sERI(nBas2,nBas2,nBas2,nBas2)) + + call spatial_to_spin_MO_energy(nBas,e,nBas2,se) + call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI) + +! Antysymmetrize ERIs + + allocate(dbERI(nBas2,nBas2,nBas2,nBas2)) + + call antisymmetrize_ERI(2,nBas2,sERI,dbERI) + + deallocate(sERI) + +! Define occupied and virtual spaces + + nO = nEl + nV = nBas2 - nO + +! Form energy denominator + + allocate(eO(nO),eV(nV)) + + eO(:) = se(1:nO) + eV(:) = se(nO+1:nBas2) + + 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:nBas2,nO+1:nBas2) + OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2,nO+1:nBas2, 1:nO ) + VVOO(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2, 1:nO , 1:nO ) + VVVV(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2) + + deallocate(dbERI) + +! Compute MP2 energy + + E2 = 0d0 + + do i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + + eps = eO(i) + eO(j) - eV(a) - eV(b) + + E2 = E2 + OOVV(i,j,a,b)*OOVV(i,j,a,b)/eps + + enddo + enddo + enddo + enddo + + 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=nO+1,nBas-nR - do b=nO+1,nBas-nR - eps1 = e(i) + e(j) - e(a) - e(b) - eps2 = e(k) + e(l) - e(a) - e(b) + do i=1,nO + do j=1,nO + do k=1,nO + do l=1,nO + do a=1,nV + do b=1,nV - E3a = E3a + (V(i,j,a,b) - V(i,j,b,a))* & - (V(k,l,i,j) - V(k,l,j,i))* & - (V(a,b,k,l) - V(a,b,l,k))/(eps1*eps2) + 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) + enddo + enddo + enddo + enddo enddo - enddo - enddo - enddo - enddo enddo E3b = 0d0 - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR - do c=nO+1,nBas-nR - do d=nO+1,nBas-nR - eps1 = e(i) + e(j) - e(a) - e(b) - eps2 = e(i) + e(j) - e(c) - e(d) + do i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + do c=1,nV + do d=1,nV - E3b = E3b + (V(i,j,a,b) - V(i,j,b,a))* & - (V(a,b,c,d) - V(a,b,d,c))* & - (V(c,d,i,j) - V(c,d,j,i))/(eps1*eps2) + 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) + enddo + enddo + enddo + enddo enddo - enddo - enddo - enddo - enddo enddo E3c = 0d0 - do i=nC+1,nO - do j=nC+1,nO - do k=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR - do c=nO+1,nBas-nR - eps1 = e(i) + e(j) - e(a) - e(b) - eps2 = e(i) + e(k) - e(a) - e(c) + do i=1,nO + do j=1,nO + do k=1,nO + do a=1,nV + do b=1,nV + do c=1,nV - E3c = E3c + (V(i,j,a,b) - V(i,j,b,a))* & - (V(k,b,c,j) - V(k,b,j,c))* & - (V(a,c,i,k) - V(a,c,k,i))/(eps1*eps2) + 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) - enddo - enddo + enddo + enddo + enddo + enddo enddo enddo - enddo - enddo - EcMP3 = 0.25d0*E3a + 0.25d0*E3b + E3c + EcMP3 = 0.125d0*E3a + 0.125d0*E3b + E3c write(*,*) write(*,'(A32)') '-----------------------' @@ -105,8 +179,8 @@ subroutine MP3(nBas,nC,nO,nV,nR,V,EHF,EcMP2,e,EcMP3) 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',EHF + EcMP2 + EcMP3 + 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(*,*)