From 66dfcd6ae7a99d9ff1a80b21cf721e5217be043e Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 10 Sep 2024 16:43:58 +0200 Subject: [PATCH 1/3] move orthogonalization matrix --- src/QuAcK/QuAcK.f90 | 64 +------------ src/utils/orthogonalization_matrix.f90 | 121 ++++++------------------- 2 files changed, 31 insertions(+), 154 deletions(-) diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index b01002f..df2b8bc 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -30,7 +30,7 @@ program QuAcK double precision,allocatable :: T(:,:) double precision,allocatable :: V(:,:) double precision,allocatable :: Hc(:,:) - double precision,allocatable :: X(:,:) + double precision,allocatable :: X(:,:),X_tmp(:,:) double precision,allocatable :: dipole_int_AO(:,:,:) double precision,allocatable :: ERI_AO(:,:,:,:) double precision,allocatable :: Uvec(:,:), Uval(:) @@ -71,10 +71,6 @@ program QuAcK logical :: dotest,doRtest,doUtest,doGtest - integer :: i, j, j0 - double precision :: acc_d, acc_nd - double precision, allocatable :: tmp1(:,:), tmp2(:,:) - !-------------! ! Hello World ! !-------------! @@ -180,61 +176,11 @@ program QuAcK ! Compute orthogonalization matrix - !call orthogonalization_matrix(nBas, S, X) - - allocate(Uvec(nBas,nBas), Uval(nBas)) - - Uvec(1:nBas,1:nBas) = S(1:nBas,1:nBas) - call diagonalize_matrix(nBas, Uvec, Uval) - - nOrb = 0 - do i = 1, nBas - if(Uval(i) > 1d-6) then - Uval(i) = 1d0 / dsqrt(Uval(i)) - nOrb = nOrb + 1 - else - write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization' - end if - end do - - write(*,'(A38)') '--------------------------------------' - write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas - write(*,'(A38,1X,I16)') 'Number of basis functions (MOs)', nOrb - write(*,'(A38,1X,F9.3)') ' % of discarded orbitals = ', 100.d0 * (1.d0 - dble(nOrb)/dble(nBas)) - write(*,'(A38)') '--------------------------------------' - write(*,*) - - j0 = nBas - nOrb + allocate(X_tmp(nBas,nBas)) + call orthogonalization_matrix(nBas,nOrb,S,X_tmp) allocate(X(nBas,nOrb)) - do j = j0+1, nBas - do i = 1, nBas - X(i,j-j0) = Uvec(i,j) * Uval(j) - enddo - enddo - - deallocate(Uvec, Uval) - - !! check if X.T S X = 1_(nOrb,nOrb) - !allocate(tmp1(nOrb,nBas), tmp2(nOrb,nOrb)) - !call dgemm("T", "N", nOrb, nBas, nBas, 1.d0, & - ! X(1,1), nBas, S(1,1), nBas, & - ! 0.d0, tmp1(1,1), nOrb) - !call dgemm("N", "N", nOrb, nOrb, nBas, 1.d0, & - ! tmp1(1,1), nOrb, X(1,1), nBas, & - ! 0.d0, tmp2(1,1), nOrb) - !acc_d = 0.d0 - !acc_nd = 0.d0 - !do i = 1, nOrb - ! !write(*,'(1000(F15.7,2X))') (tmp2(i,j), j = 1, nOrb) - ! acc_d = acc_d + tmp2(i,i) - ! do j = 1, nOrb - ! if(j == i) cycle - ! acc_nd = acc_nd + dabs(tmp2(j,i)) - ! enddo - !enddo - !print*, ' diag part: ', dabs(acc_d - dble(nOrb)) / dble(nOrb) - !print*, ' non-diag part: ', acc_nd - !deallocate(tmp1, tmp2) + X(1:nBas,1:nOrb) = X_tmp(1:nBas,1:nOrb) + deallocate(X_tmp) !---------------------! ! Choose QuAcK branch ! diff --git a/src/utils/orthogonalization_matrix.f90 b/src/utils/orthogonalization_matrix.f90 index af781e1..fd9e59c 100644 --- a/src/utils/orthogonalization_matrix.f90 +++ b/src/utils/orthogonalization_matrix.f90 @@ -1,7 +1,4 @@ - -! --- - -subroutine orthogonalization_matrix(nBas, S, X) +subroutine orthogonalization_matrix(nBas,nOrb,S,X) ! Compute the orthogonalization matrix X @@ -17,113 +14,47 @@ subroutine orthogonalization_matrix(nBas, S, X) logical :: debug double precision,allocatable :: UVec(:,:),Uval(:) double precision,parameter :: thresh = 1d-6 - integer,parameter :: ortho_type = 1 - integer :: i + integer :: i,j,j0 ! Output variables + integer :: nOrb double precision,intent(out) :: X(nBas,nBas) debug = .false. -! Type of orthogonalization ortho_type -! -! 1 = Lowdin -! 2 = Canonical -! 3 = SVD -! - allocate(Uvec(nBas,nBas),Uval(nBas)) - if(ortho_type == 1) then - - ! - ! S V = V s where - ! - ! V.T V = 1 and s > 0 (S is positive def) - ! - ! S = V s V.T - ! = V s^0.5 s^0.5 V.T - ! = V s^0.5 V.T V s^0.5 V.T - ! = S^0.5 S^0.5 - ! - ! where - ! - ! S^0.5 = V s^0.5 V.T - ! - ! X = S^(-0.5) - ! = V s^(-0.5) V.T - ! - -! write(*,*) -! write(*,*) ' Lowdin orthogonalization' -! write(*,*) - - Uvec = S - call diagonalize_matrix(nBas, Uvec, Uval) - - do i = 1, nBas - - if(Uval(i) < thresh) then - - write(*,*) 'Eigenvalue',i,' is very small in Lowdin orthogonalization = ',Uval(i) - - end if - - Uval(i) = 1d0 / dsqrt(Uval(i)) - - end do - - call ADAt(nBas, Uvec(1,1), Uval(1), X(1,1)) - - elseif(ortho_type == 2) then - -! write(*,*) -! write(*,*) 'Canonical orthogonalization' -! write(*,*) - - Uvec = S - call diagonalize_matrix(nBas, Uvec, Uval) - - do i = 1, nBas - - if(Uval(i) > thresh) then + Uvec(1:nBas,1:nBas) = S(1:nBas,1:nBas) + call diagonalize_matrix(nBas,Uvec,Uval) + nOrb = 0 + do i = 1,nBas + if(Uval(i) > thresh) then Uval(i) = 1d0 / dsqrt(Uval(i)) + nOrb = nOrb + 1 + else + write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization' + end if + end do - else + write(*,'(A50)') '------------------------------------------------' + write(*,'(A40,1X,I5)') 'Number of basis functions = ',nBas + write(*,'(A40,1X,I5)') 'Number of spatial orbitals = ',nOrb + write(*,'(A40,1X,I5)') 'Number of discarded functions = ',nBas - nOrb + write(*,'(A50)') '------------------------------------------------' + write(*,*) - write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization' + j0 = nBas - nOrb - end if + do j = j0+1,nBas + do i = 1,nBas + X(i,j-j0) = Uvec(i,j) * Uval(j) + enddo + enddo - end do - - call AD(nBas, Uvec, Uval) - X = Uvec - - elseif(ortho_type == 3) then - -! write(*,*) -! write(*,*) ' SVD-based orthogonalization NYI' -! write(*,*) - -! Uvec = S -! call diagonalize_matrix(nBas,Uvec,Uval) - -! do i=1,nBas -! if(Uval(i) > thresh) then -! Uval(i) = 1d0/sqrt(Uval(i)) -! else -! write(*,*) 'Eigenvalue',i,'too small for canonical orthogonalization' -! end if -! end do -! -! call AD(nBas,Uvec,Uval) -! X = Uvec - - end if + deallocate(Uvec,Uval) ! Print results From 66c59f85a1c2b6cff6acd5533e88da6aea649027 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 10 Sep 2024 17:05:03 +0200 Subject: [PATCH 2/3] remove unnecessary code in AOtoMO --- src/AOtoMO/AOtoMO_oooa.f90 | 85 ------------------------- src/AOtoMO/AOtoMO_oooo.f90 | 85 ------------------------- src/AOtoMO/AOtoMO_oovv.f90 | 77 ---------------------- src/AOtoMO/Hartree_matrix_MO_basis.f90 | 26 -------- src/AOtoMO/exchange_matrix_MO_basis.f90 | 26 -------- 5 files changed, 299 deletions(-) delete mode 100644 src/AOtoMO/AOtoMO_oooa.f90 delete mode 100644 src/AOtoMO/AOtoMO_oooo.f90 delete mode 100644 src/AOtoMO/AOtoMO_oovv.f90 delete mode 100644 src/AOtoMO/Hartree_matrix_MO_basis.f90 delete mode 100644 src/AOtoMO/exchange_matrix_MO_basis.f90 diff --git a/src/AOtoMO/AOtoMO_oooa.f90 b/src/AOtoMO/AOtoMO_oooa.f90 deleted file mode 100644 index c8b60ea..0000000 --- a/src/AOtoMO/AOtoMO_oooa.f90 +++ /dev/null @@ -1,85 +0,0 @@ -subroutine AOtoMO_oooa(nBas,nO,nA,cO,cA,O,ooOoa) - -! AO to MO transformation of two-electron integrals for the block oooa -! Semi-direct O(N^5) algorithm - - implicit none - -! Input variables - - integer,intent(in) :: nBas,nO,nA - double precision,intent(in) :: cO(nBas,nO),cA(nBas,nA),O(nBas,nBas,nBas,nBas) - -! Local variables - - double precision,allocatable :: scr1(:,:,:,:),scr2(:,:,:,:) - integer :: mu,nu,la,si,i,j,k,x - -! Output variables - - double precision,intent(out) :: ooOoa(nO,nO,nO,nA) - -! Memory allocation - allocate(scr1(nBas,nBas,nBas,nBas),scr2(nBas,nBas,nBas,nBas)) - - write(*,*) - write(*,'(A42)') '----------------------------------------' - write(*,'(A42)') ' AO to MO transformation for oooa block ' - write(*,'(A42)') '----------------------------------------' - write(*,*) - - scr1 = 0d0 - do mu=1,nBas - do nu=1,nBas - do la=1,nBas - do si=1,nBas - do x=1,nA - scr1(mu,nu,la,x) = scr1(mu,nu,la,x) + O(mu,nu,la,si)*cA(si,x) - end do - end do - end do - end do - end do - - scr2 = 0d0 - do mu=1,nBas - do nu=1,nBas - do la=1,nBas - do i=1,nO - do x=1,nA - scr2(i,nu,la,x) = scr2(i,nu,la,x) + cO(mu,i)*scr1(mu,nu,la,x) - end do - end do - end do - end do - end do - - scr1 = 0d0 - do nu=1,nBas - do la=1,nBas - do i=1,nO - do k=1,nO - do x=1,nA - scr1(i,nu,k,x) = scr1(i,nu,k,x) + scr2(i,nu,la,x)*cO(la,k) - end do - end do - end do - end do - end do - - ooOoa = 0d0 - do nu=1,nBas - do i=1,nO - do j=1,nO - do k=1,nO - do x=1,nA - ooOoa(i,j,k,x) = ooOoa(i,j,k,x) + cO(nu,j)*scr1(i,nu,k,x) - end do - end do - end do - end do - end do - - deallocate(scr1,scr2) - -end subroutine diff --git a/src/AOtoMO/AOtoMO_oooo.f90 b/src/AOtoMO/AOtoMO_oooo.f90 deleted file mode 100644 index 75286bd..0000000 --- a/src/AOtoMO/AOtoMO_oooo.f90 +++ /dev/null @@ -1,85 +0,0 @@ -subroutine AOtoMO_oooo(nBas,nO,cO,O,ooOoo) - -! AO to MO transformation of two-electron integrals for the block oooo -! Semi-direct O(N^5) algorithm - - implicit none - -! Input variables - - integer,intent(in) :: nBas,nO - double precision,intent(in) :: cO(nBas,nO),O(nBas,nBas,nBas,nBas) - -! Local variables - - double precision,allocatable :: scr1(:,:,:,:),scr2(:,:,:,:) - integer :: mu,nu,la,si,i,j,k,l - -! Output variables - - double precision,intent(out) :: ooOoo(nO,nO,nO,nO) - -! Memory allocation - allocate(scr1(nBas,nBas,nBas,nBas),scr2(nBas,nBas,nBas,nBas)) - - write(*,*) - write(*,'(A42)') '----------------------------------------' - write(*,'(A42)') ' AO to MO transformation for oooo block ' - write(*,'(A42)') '----------------------------------------' - write(*,*) - - scr1 = 0d0 - do mu=1,nBas - do nu=1,nBas - do la=1,nBas - do si=1,nBas - do l=1,nO - scr1(mu,nu,la,l) = scr1(mu,nu,la,l) + O(mu,nu,la,si)*cO(si,l) - end do - end do - end do - end do - end do - - scr2 = 0d0 - do mu=1,nBas - do nu=1,nBas - do la=1,nBas - do i=1,nO - do l=1,nO - scr2(i,nu,la,l) = scr2(i,nu,la,l) + cO(mu,i)*scr1(mu,nu,la,l) - end do - end do - end do - end do - end do - - scr1 = 0d0 - do nu=1,nBas - do la=1,nBas - do i=1,nO - do k=1,nO - do l=1,nO - scr1(i,nu,k,l) = scr1(i,nu,k,l) + scr2(i,nu,la,l)*cO(la,k) - end do - end do - end do - end do - end do - - ooOoo = 0d0 - do nu=1,nBas - do i=1,nO - do j=1,nO - do k=1,nO - do l=1,nO - ooOoo(i,j,k,l) = ooOoo(i,j,k,l) + cO(nu,j)*scr1(i,nu,k,l) - end do - end do - end do - end do - end do - - deallocate(scr1,scr2) - -end subroutine diff --git a/src/AOtoMO/AOtoMO_oovv.f90 b/src/AOtoMO/AOtoMO_oovv.f90 deleted file mode 100644 index ad5e5d3..0000000 --- a/src/AOtoMO/AOtoMO_oovv.f90 +++ /dev/null @@ -1,77 +0,0 @@ -subroutine AOtoMO_oovv(nBas,nO,nV,cO,cV,O,ooOvv) - -! AO to MO transformation of two-electron integrals for the block oovv -! Semi-direct O(N^5) algorithm - - implicit none - -! Input variables - - integer,intent(in) :: nBas,nO,nV - double precision,intent(in) :: cO(nBas,nO),cV(nBas,nV),O(nBas,nBas,nBas,nBas) - -! Local variables - - double precision,allocatable :: scr1(:,:,:,:),scr2(:,:,:,:) - integer :: mu,nu,la,si,i,j,a,b - -! Output variables - - double precision,intent(out) :: ooOvv(nO,nO,nV,nV) - -! Memory allocation - allocate(scr1(nBas,nBas,nBas,nBas),scr2(nBas,nBas,nBas,nBas)) - - scr1 = 0d0 - do mu=1,nBas - do nu=1,nBas - do la=1,nBas - do si=1,nBas - do b=1,nV - scr1(mu,nu,la,b) = scr1(mu,nu,la,b) + O(mu,nu,la,si)*cV(si,b) - end do - end do - end do - end do - end do - - scr2 = 0d0 - do mu=1,nBas - do nu=1,nBas - do la=1,nBas - do i=1,nO - do b=1,nV - scr2(i,nu,la,b) = scr2(i,nu,la,b) + cO(mu,i)*scr1(mu,nu,la,b) - end do - end do - end do - end do - end do - - scr1 = 0d0 - do nu=1,nBas - do la=1,nBas - do i=1,nO - do a=1,nV - do b=1,nV - scr1(i,nu,a,b) = scr1(i,nu,a,b) + scr2(i,nu,la,b)*cV(la,a) - end do - end do - end do - end do - end do - - ooOvv = 0d0 - do nu=1,nBas - do i=1,nO - do j=1,nO - do a=1,nV - do b=1,nV - ooOvv(i,j,a,b) = ooOvv(i,j,a,b) + cO(nu,j)*scr1(i,nu,a,b) - end do - end do - end do - end do - end do - -end subroutine diff --git a/src/AOtoMO/Hartree_matrix_MO_basis.f90 b/src/AOtoMO/Hartree_matrix_MO_basis.f90 deleted file mode 100644 index c1496b6..0000000 --- a/src/AOtoMO/Hartree_matrix_MO_basis.f90 +++ /dev/null @@ -1,26 +0,0 @@ -subroutine Hartree_matrix_MO_basis(nBas,c,P,G,H) - -! Compute Hartree matrix in the MO basis - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nBas - double precision,intent(in) :: c(nBas,nBas),P(nBas,nBas) - double precision,intent(in) :: G(nBas,nBas,nBas,nBas) - -! Output variables - - double precision,intent(out) :: H(nBas,nBas) - -! Compute Hartree matrix in the AO basis - - call Hartree_matrix_AO_basis(nBas,P,G,H) - -! Transform Hartree matrix in the MO basis - - H = matmul(transpose(c),matmul(H,c)) - -end subroutine diff --git a/src/AOtoMO/exchange_matrix_MO_basis.f90 b/src/AOtoMO/exchange_matrix_MO_basis.f90 deleted file mode 100644 index 5cb13b1..0000000 --- a/src/AOtoMO/exchange_matrix_MO_basis.f90 +++ /dev/null @@ -1,26 +0,0 @@ -subroutine exchange_matrix_MO_basis(nBas,c,P,G,K) - -! Compute exchange matrix in the MO basis - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nBas - double precision,intent(in) :: c(nBas,nBas),P(nBas,nBas) - double precision,intent(in) :: G(nBas,nBas,nBas,nBas) - -! Output variables - - double precision,intent(out) :: K(nBas,nBas) - -! Compute Hartree Hamiltonian in the AO basis - - call exchange_matrix_AO_basis(nBas,P,G,K) - -! Transform Coulomb matrix in the MO basis - - K = matmul(transpose(c),matmul(K,c)) - -end subroutine exchange_matrix_MO_basis From 34fba9ed9c3d3364abfd92113cb6a5a90f91e1ae Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 10 Sep 2024 17:08:16 +0200 Subject: [PATCH 3/3] remove RMP3 --- src/MP/RMP.f90 | 7 +- src/MP/RMP3.f90 | 194 ------------------------------------------------ 2 files changed, 4 insertions(+), 197 deletions(-) delete mode 100644 src/MP/RMP3.f90 diff --git a/src/MP/RMP.f90 b/src/MP/RMP.f90 index 01e2288..8b08858 100644 --- a/src/MP/RMP.f90 +++ b/src/MP/RMP.f90 @@ -31,7 +31,7 @@ subroutine RMP(dotest,doMP2,doMP3,regularize,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) ! Output variables !------------------------------------------------------------------------ -! Compute MP3 energy +! Compute MP2 energy !------------------------------------------------------------------------ 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 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) 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(*,*) end if diff --git a/src/MP/RMP3.f90 b/src/MP/RMP3.f90 deleted file mode 100644 index 788737b..0000000 --- a/src/MP/RMP3.f90 +++ /dev/null @@ -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