From 5aaa2addcfc51908068645de4c413af7bdbc81e6 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sun, 17 Mar 2019 15:37:55 +0100 Subject: [PATCH] adding missing files --- src/MCQC/CCD.f90 | 193 ++++++++++++++++++ src/MCQC/CCSD.f90 | 259 +++++++++++++++++++++++++ src/MCQC/CCSDT.f90 | 45 +++++ src/MCQC/EOM_CCSD.f90 | 78 ++++++++ src/MCQC/chem_to_phys_ERI.f90 | 33 ++++ src/MCQC/form_T.f90 | 46 +++++ src/MCQC/form_X.f90 | 92 +++++++++ src/MCQC/form_abh.f90 | 105 ++++++++++ src/MCQC/form_delta_OOOVVV.f90 | 37 ++++ src/MCQC/form_delta_OOVV.f90 | 33 ++++ src/MCQC/form_delta_OV.f90 | 27 +++ src/MCQC/form_g.f90 | 53 +++++ src/MCQC/form_h.f90 | 79 ++++++++ src/MCQC/form_r1.f90 | 77 ++++++++ src/MCQC/form_r2.f90 | 139 +++++++++++++ src/MCQC/form_tau.f90 | 34 ++++ src/MCQC/form_ub.f90 | 48 +++++ src/MCQC/form_ubb.f90 | 67 +++++++ src/MCQC/form_v.f90 | 79 ++++++++ src/MCQC/spatial_to_spin_ERI.f90 | 35 ++++ src/MCQC/spatial_to_spin_MO_energy.f90 | 29 +++ 21 files changed, 1588 insertions(+) create mode 100644 src/MCQC/CCD.f90 create mode 100644 src/MCQC/CCSD.f90 create mode 100644 src/MCQC/CCSDT.f90 create mode 100644 src/MCQC/EOM_CCSD.f90 create mode 100644 src/MCQC/chem_to_phys_ERI.f90 create mode 100644 src/MCQC/form_T.f90 create mode 100644 src/MCQC/form_X.f90 create mode 100644 src/MCQC/form_abh.f90 create mode 100644 src/MCQC/form_delta_OOOVVV.f90 create mode 100644 src/MCQC/form_delta_OOVV.f90 create mode 100644 src/MCQC/form_delta_OV.f90 create mode 100644 src/MCQC/form_g.f90 create mode 100644 src/MCQC/form_h.f90 create mode 100644 src/MCQC/form_r1.f90 create mode 100644 src/MCQC/form_r2.f90 create mode 100644 src/MCQC/form_tau.f90 create mode 100644 src/MCQC/form_ub.f90 create mode 100644 src/MCQC/form_ubb.f90 create mode 100644 src/MCQC/form_v.f90 create mode 100644 src/MCQC/spatial_to_spin_ERI.f90 create mode 100644 src/MCQC/spatial_to_spin_MO_energy.f90 diff --git a/src/MCQC/CCD.f90 b/src/MCQC/CCD.f90 new file mode 100644 index 0000000..de3ff2c --- /dev/null +++ b/src/MCQC/CCD.f90 @@ -0,0 +1,193 @@ +subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) + +! CCD module + + implicit none + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas,nEl + double precision,intent(in) :: ENuc,ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: nBas2 + integer :: nO + integer :: nV + integer :: nSCF + double precision :: Conv + double precision :: EcMP2 + double precision :: ECCD,EcCCD + double precision,allocatable :: seHF(:) + double precision,allocatable :: sERI(:,:,:,:) + double precision,allocatable :: dbERI(:,:,:,:) + + double precision,allocatable :: eO(:) + double precision,allocatable :: eV(:) + double precision,allocatable :: delta_OOVV(:,:,:,:) + + double precision,allocatable :: OOOO(:,:,:,:) + double precision,allocatable :: OOVV(:,:,:,:) + double precision,allocatable :: OVOV(:,:,:,:) + double precision,allocatable :: VVVV(:,:,:,:) + + double precision,allocatable :: X1(:,:,:,:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: X3(:,:) + double precision,allocatable :: X4(:,:,:,:) + + double precision,allocatable :: u(:,:,:,:) + double precision,allocatable :: v(:,:,:,:) + + double precision,allocatable :: r2(:,:,:,:) + double precision,allocatable :: t2(:,:,:,:) + +! Hello world + + write(*,*) + write(*,*)'**************************************' + write(*,*)'| CCD calculation |' + write(*,*)'**************************************' + write(*,*) + +! Spatial to spin orbitals + + nBas2 = 2*nBas + + allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2)) + + call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF) + 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)) + allocate(delta_OOVV(nO,nO,nV,nV)) + + eO(:) = seHF(1:nO) + eV(:) = seHF(nO+1:nBas2) + + call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV) + + deallocate(seHF) + +! Create integral batches + + allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV),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) + OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2, 1:nO ,nO+1:nBas2) + VVVV(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2) + + deallocate(dbERI) + +! MP2 guess amplitudes + + allocate(t2(nO,nO,nV,nV)) + + t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + + EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) + write(*,'(1X,A10,1X,F10.6)') 'Ec(MP2) = ',EcMP2 + +! Initialization + + allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV)) + allocate(X1(nO,nO,nO,nO),X2(nV,nV),X3(nO,nO),X4(nO,nO,nV,nV)) + + Conv = 1d0 + nSCF = 0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)'| CCD calculation |' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(CCD)','|','Ec(CCD)','|','Conv','|' + write(*,*)'----------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + +! Increment + + nSCF = nSCF + 1 + +! Form linear array + + call form_u(nO,nV,OOOO,VVVV,OVOV,t2,u) + +! Form interemediate arrays + + call form_X(nO,nV,OOVV,t2,X1,X2,X3,X4) + +! Form quadratic array + + call form_v(nO,nV,X1,X2,X3,X4,t2,v) + +! Compute residual + + r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + u(:,:,:,:) + v(:,:,:,:) + +! Check convergence + + Conv = maxval(abs(r2(:,:,:,:))) + +! Update amplitudes + + t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:) + +! Compute correlation energy + + EcCCD = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) + +! Dump results + + ECCD = ERHF + EcCCD + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|' + + enddo + write(*,*)'----------------------------------------------------' +!------------------------------------------------------------------------ +! End of SCF loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + endif + +end subroutine CCD diff --git a/src/MCQC/CCSD.f90 b/src/MCQC/CCSD.f90 new file mode 100644 index 0000000..124b2d9 --- /dev/null +++ b/src/MCQC/CCSD.f90 @@ -0,0 +1,259 @@ +subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF) + +! CCSD module + + implicit none + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + logical,intent(in) :: doCCSDT + integer,intent(in) :: nBas,nEl + double precision,intent(in) :: ENuc,ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: start_CCSDT,end_CCSDT,t_CCSDT + integer :: nBas2 + integer :: nO + integer :: nV + integer :: nSCF + double precision :: Conv + double precision :: EcMP2 + double precision :: ECCSD,EcCCSD + double precision :: EcCCT + + double precision,allocatable :: seHF(:) + double precision,allocatable :: sERI(:,:,:,:) + double precision,allocatable :: dbERI(:,:,:,:) + double precision,allocatable :: delta_OV(:,:) + double precision,allocatable :: delta_OOVV(:,:,:,:) + + double precision,allocatable :: OOOO(:,:,:,:) + double precision,allocatable :: OOOV(:,:,:,:) + double precision,allocatable :: OVOO(:,:,:,:) + double precision,allocatable :: VOOO(:,:,:,:) + double precision,allocatable :: OOVV(:,:,:,:) + double precision,allocatable :: OVVO(:,:,:,:) + double precision,allocatable :: OVVV(:,:,:,:) + double precision,allocatable :: VOVV(:,:,:,:) + double precision,allocatable :: VVVO(:,:,:,:) + double precision,allocatable :: VVVV(:,:,:,:) + + double precision,allocatable :: eO(:) + double precision,allocatable :: eV(:) + double precision,allocatable :: hvv(:,:) + double precision,allocatable :: hoo(:,:) + double precision,allocatable :: hvo(:,:) + double precision,allocatable :: gvv(:,:) + double precision,allocatable :: goo(:,:) + double precision,allocatable :: aoooo(:,:,:,:) + double precision,allocatable :: bvvvv(:,:,:,:) + double precision,allocatable :: hovvo(:,:,:,:) + + double precision,allocatable :: r1(:,:) + double precision,allocatable :: r2(:,:,:,:) + + double precision,allocatable :: t1(:,:) + double precision,allocatable :: t2(:,:,:,:) + double precision,allocatable :: tau(:,:,:,:) + +! Hello world + + write(*,*) + write(*,*)'**************************************' + write(*,*)'| CCSD calculation |' + write(*,*)'**************************************' + write(*,*) + +! Spatial to spin orbitals + + nBas2 = 2*nBas + + allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2)) + + call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF) + 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)) + allocate(delta_OV(nO,nV),delta_OOVV(nO,nO,nV,nV)) + + eO(:) = seHF(1:nO) + eV(:) = seHF(nO+1:nBas2) + + call form_delta_OV(nO,nV,eO,eV,delta_OV) + call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV) + + deallocate(seHF) + +! Create integral batches + + allocate(OOOO(nO,nO,nO,nO), & + OOOV(nO,nO,nO,nV),OVOO(nO,nV,nO,nO),VOOO(nV,nO,nO,nO), & + OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO), & + OVVV(nO,nV,nV,nV),VOVV(nV,nO,nV,nV),VVVO(nV,nV,nV,nO), & + VVVV(nV,nV,nV,nV)) + + OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO ) + OOOV(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO ,nO+1:nBas2) + OVOO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2, 1:nO , 1:nO ) + VOOO(:,:,:,:) = dbERI(nO+1:nBas2, 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 ) + OVVV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2) + VOVV(:,:,:,:) = dbERI(nO+1:nBas2, 1:nO ,nO+1:nBas2,nO+1:nBas2) + VVVO(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2, 1:nO ) + VVVV(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2) + + deallocate(dbERI) + +! MP2 guess amplitudes + + allocate(t1(nO,nV),t2(nO,nO,nV,nV),tau(nO,nO,nV,nV)) + + t1(:,:) = 0d0 + t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + call form_tau(nO,nV,t1,t2,tau) + + EcMP2 = 0.5d0*dot_product(pack(OOVV,.true.),pack(tau,.true.)) + write(*,'(1X,A10,1X,F10.6)') 'Ec(MP2) = ',EcMP2 + +! Initialization + + allocate(hvv(nV,nV),hoo(nO,nO),hvo(nV,nO), & + gvv(nV,nV),goo(nO,nO), & + aoooo(nO,nO,nO,nO),bvvvv(nV,nV,nV,nV),hovvo(nO,nV,nV,nO), & + r1(nO,nV),r2(nO,nO,nV,nV)) + + Conv = 1d0 + nSCF = 0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)'| CCSD calculation |' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(CCSD)','|','Ec(CCSD)','|','Conv','|' + write(*,*)'----------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + +! Increment + + nSCF = nSCF + 1 + +! Scuseria Eqs. (5), (6) and (7) + + call form_h(nO,nV,eO,eV,OOVV,t1,tau,hvv,hoo,hvo) + +! Scuseria Eqs. (9), (10), (11), (12) and (13) + + call form_g(nO,nV,hvv,hoo,VOVV,OOOV,t1,gvv,goo) + + call form_abh(nO,nV,OOOO,OVOO,OOVV,VVVV,VOVV,OVVO,OVVV,t1,tau,aoooo,bvvvv,hovvo) + +! Compute residuals + + call form_r1(nO,nV,OVVO,OVVV,OOOV,hvv,hoo,hvo,t1,t2,tau,r1) + + call form_r2(nO,nV,OOVV,OVOO,OVVV,OVVO,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau,r2) + +! Check convergence + + Conv = max(maxval(abs(r1(:,:))),maxval(abs(r2(:,:,:,:)))) + +! Update + + t1(:,:) = t1(:,:) - r1(:,:) /delta_OV (:,:) + t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:) + + call form_tau(nO,nV,t1,t2,tau) + +! Compute correlation energy + + EcCCSD = 0.5d0*dot_product(pack(OOVV,.true.),pack(tau,.true.)) + +! Dump results + + ECCSD = ERHF + EcCCSD + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nSCF,'|',ECCSD+ENuc,'|',EcCCSD,'|',Conv,'|' + + end do + write(*,*)'----------------------------------------------------' +!------------------------------------------------------------------------ +! End of SCF loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + end if + +! Deallocate memory + + deallocate(hvv,hoo,hvo, & + delta_OV,delta_OOVV, & + gvv,goo, & + aoooo,bvvvv,hovvo, & + tau, & + r1,r2) + +!------------------------------------------------------------------------ +! (T) correction +!------------------------------------------------------------------------ + if(doCCSDT) then + + call cpu_time(start_CCSDT) + call CCSDT(nO,nV,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT) + call cpu_time(end_CCSDT) + + t_CCSDT = end_CCSDT - start_CCSDT + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for (T) = ',t_CCSDT,' seconds' + write(*,*) + + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)' CCSDT(T) energy ' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A20,1X,F15.10)')' E(CCSD(T)) = ',ECCSD + EcCCT + write(*,'(1X,A20,1X,F10.6)') ' Ec(CCSD(T)) = ',EcCCSD + EcCCT + write(*,*)'----------------------------------------------------' + write(*,*) + + end if + +end subroutine CCSD diff --git a/src/MCQC/CCSDT.f90 b/src/MCQC/CCSDT.f90 new file mode 100644 index 0000000..2af9b3b --- /dev/null +++ b/src/MCQC/CCSDT.f90 @@ -0,0 +1,45 @@ +subroutine CCSDT(nO,nV,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT) + +! Compute the (T) correction of the CCSD(T) energy + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + + double precision,intent(in) :: eO(nO) + double precision,intent(in) :: eV(nV) + + double precision,intent(in) :: OOVV(nO,nO,nV,nV) + double precision,intent(in) :: VVVO(nV,nV,nV,nO) + double precision,intent(in) :: VOOO(nV,nO,nO,nO) + + double precision,intent(in) :: t1(nO,nV) + double precision,intent(in) :: t2(nO,nO,nV,nV) + +! Local variables + + double precision,allocatable :: delta_OOOVVV(:,:,:,:,:,:) + double precision,allocatable :: ub(:,:,:,:,:,:) + double precision,allocatable :: ubb(:,:,:,:,:,:) + +! Output variables + + double precision,intent(out) :: EcCCT + +! Memory allocation + + allocate(delta_OOOVVV(nO,nO,nO,nV,nV,nV),ub(nO,nO,nO,nV,nV,nV),ubb(nO,nO,nO,nV,nV,nV)) + +! Form CCSD(T) quantities + + call form_delta_OOOVVV(nO,nV,eO,eV,delta_OOOVVV) + + call form_ub(nO,nV,OOVV,t1,ub) + + call form_ubb(nO,nV,VVVO,VOOO,t2,ubb) + + call form_T(nO,nV,delta_OOOVVV,ub,ubb,EcCCT) + +end subroutine CCSDT diff --git a/src/MCQC/EOM_CCSD.f90 b/src/MCQC/EOM_CCSD.f90 new file mode 100644 index 0000000..2be73aa --- /dev/null +++ b/src/MCQC/EOM_CCSD.f90 @@ -0,0 +1,78 @@ +subroutine EOM_CCSD(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) + +! Compute EOM-CCSD excitation energies: see Stanton & Bartlett JCP 98 7029 (1993) + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: maxSCF + double precision,intent(in) :: thresh + integer,intent(in) :: max_diis + integer,intent(in) :: nBas,nC,nO,nV,nR + double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: nH,nP,nHH,nPP,nSCF,n_diis + double precision :: Conv + double precision,external :: Kronecker_delta + double precision,allocatable :: B_ADC(:,:),X_ADC(:,:),e_ADC(:),SigInf(:,:),G_ADC(:,:) + double precision,allocatable :: db_ERI(:,:,:,:),eOld(:),error_diis(:,:),e_diis(:,:) + + integer :: i,j,k,l + integer :: a,b,c,d + integer :: p,q,r,s + integer :: nADC,iADC,jADC + + +! Hello world + + write(*,*) + write(*,*)'***********************************' + write(*,*)'| EOM-CCSD calculation |' + write(*,*)'***********************************' + write(*,*) + +! Number of holes + + nH = nO + nHH = nH*nH + +! Number of particles + + nP = nV + nPP = nP*nP + + write(*,*) 'Total states: ',nH + nP + write(*,*) 'Hole states: ',nH + write(*,*) 'Particle states: ',nP + +! Size of EOM-CCSD matrices + + nEOM = nH + nP + nH*nPP + nHH*nP + nHH*nPP + write(*,'(1X,A25,I3,A6,I6)') 'Size of EOM-CCSD matrix: ',nEOM,' x ',nEOM + +! Memory allocation + + allocate() + +! Construct EOM-CCSD matrix + + H(:,:) = 0d0 + + iEOM = 1 + jEOM = 1 + + H(iEOM,jEOM) = ECCSD + + do p=1,nO + + jADC = jADC + 1 + B_ADC(jADC,jADC) = e(p) + + enddo + +end subroutine EOM_CCSD diff --git a/src/MCQC/chem_to_phys_ERI.f90 b/src/MCQC/chem_to_phys_ERI.f90 new file mode 100644 index 0000000..9764bcb --- /dev/null +++ b/src/MCQC/chem_to_phys_ERI.f90 @@ -0,0 +1,33 @@ +subroutine chem_to_phys_ERI(nBas,ERI) + +! Antisymmetrize ERIs + + implicit none + +! Input variables + + integer,intent(in) :: nBas + double precision,intent(inout):: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: p,q,r,s + double precision,allocatable :: pERI(:,:,:,:) + + allocate(pERI(nBas,nBas,nBas,nBas)) + + do p=1,nBas + do q=1,nBas + do r=1,nBas + do s=1,nBas + + pERI(p,q,r,s) = ERI(p,r,q,s) + + enddo + enddo + enddo + enddo + + ERI(:,:,:,:) = pERI(:,:,:,:) + +end subroutine chem_to_phys_ERI diff --git a/src/MCQC/form_T.f90 b/src/MCQC/form_T.f90 new file mode 100644 index 0000000..dc66b63 --- /dev/null +++ b/src/MCQC/form_T.f90 @@ -0,0 +1,46 @@ +subroutine form_T(nO,nV,delta_OOOVVV,ub,ubb,EcCCT) + +! Compute (T) correction + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + + double precision,intent(in) :: delta_OOOVVV(nO,nO,nO,nV,nV,nV) + double precision,intent(in) :: ub(nO,nO,nO,nV,nV,nV) + double precision,intent(in) :: ubb(nO,nO,nO,nV,nV,nV) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + +! Output variables + + double precision,intent(out) :: EcCCT + + EcCCT = 0d0 + + do i=1,nO + do j=1,nO + do k=1,nO + do a=1,nV + do b=1,nV + do c=1,nV + + EcCCT = EcCCT & + + (ub(i,j,k,a,b,c) + ubb(i,j,k,a,b,c)) & + * ubb(i,j,k,a,b,c)/delta_OOOVVV(i,j,k,a,b,c) + + end do + end do + end do + end do + end do + end do + + EcCCT = - EcCCT/36d0 + +end subroutine form_T diff --git a/src/MCQC/form_X.f90 b/src/MCQC/form_X.f90 new file mode 100644 index 0000000..9f2eeb5 --- /dev/null +++ b/src/MCQC/form_X.f90 @@ -0,0 +1,92 @@ +subroutine form_X(nO,nV,OOVV,t2,X1,X2,X3,X4) + +! Form intermediate arrays X's in CCD + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + double precision,intent(in) :: t2(nO,nO,nV,nV) + double precision,intent(in) :: OOVV(nO,nO,nV,nV) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + +! Output variables + + double precision,intent(out) :: X1(nO,nO,nO,nO) + double precision,intent(out) :: X2(nV,nV) + double precision,intent(out) :: X3(nO,nO) + double precision,intent(out) :: X4(nO,nO,nV,nV) + +! Initialization + + X1(:,:,:,:) = 0d0 + X2(:,:) = 0d0 + X3(:,:) = 0d0 + X4(:,:,:,:) = 0d0 + +! Build X1 + + do k=1,nO + do l=1,nO + do i=1,nO + do j=1,nO + do c=1,nV + do d=1,nV + X1(k,l,i,j) = X1(k,l,i,j) + OOVV(k,l,c,d)*t2(i,j,c,d) + enddo + enddo + enddo + enddo + enddo + enddo + +! Build X2 + + do b=1,nV + do c=1,nV + do k=1,nO + do l=1,nO + do d=1,nV + X2(b,c) = X2(b,c) + OOVV(k,l,c,d)*t2(k,l,b,d) + enddo + enddo + enddo + enddo + enddo + +! Build X3 + + do k=1,nO + do j=1,nO + do l=1,nO + do c=1,nV + do d=1,nV + X3(k,j) = X3(k,j) + OOVV(k,l,c,d)*t2(j,l,c,d) + enddo + enddo + enddo + enddo + enddo + +! Build X4 + + do i=1,nO + do l=1,nO + do a=1,nV + do d=1,nV + do k=1,nO + do c=1,nV + X4(i,l,a,d) = X4(i,l,a,d) + OOVV(k,l,c,d)*t2(i,k,a,c) + enddo + enddo + enddo + enddo + enddo + enddo + +end subroutine form_X diff --git a/src/MCQC/form_abh.f90 b/src/MCQC/form_abh.f90 new file mode 100644 index 0000000..4b0efcc --- /dev/null +++ b/src/MCQC/form_abh.f90 @@ -0,0 +1,105 @@ +subroutine form_abh(nO,nV,OOOO,OVOO,OOVV,VVVV,VOVV,OVVO,OVVV,t1,tau,aoooo,bvvvv,hovvo) + +! Scuseria Eqs. (11),(12) and (13) + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + + double precision,intent(in) :: OOOO(nO,nO,nO,nO) + double precision,intent(in) :: OVOO(nO,nV,nO,nO) + double precision,intent(in) :: OOVV(nO,nO,nV,nV) + double precision,intent(in) :: OVVO(nO,nV,nV,nO) + double precision,intent(in) :: OVVV(nO,nV,nV,nV) + double precision,intent(in) :: VOVV(nV,nO,nV,nV) + double precision,intent(in) :: VVVV(nV,nV,nV,nV) + + double precision,intent(in) :: t1(nO,nV) + double precision,intent(in) :: tau(nO,nO,nV,nV) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + +! Output variables + + double precision,intent(out) :: aoooo(nO,nO,nO,nO) + double precision,intent(out) :: bvvvv(nV,nV,nV,nV) + double precision,intent(out) :: hovvo(nO,nV,nV,nO) + + aoooo(:,:,:,:) = OOOO(:,:,:,:) + + do i=1,nO + do j=1,nO + do k=1,nO + do l=1,nO + + do c=1,nV + aoooo(i,j,k,l) = aoooo(i,j,k,l) + OVOO(i,c,k,l)*t1(j,c) + end do + + do c=1,nV + aoooo(i,j,k,l) = aoooo(i,j,k,l) - OVOO(j,c,k,l)*t1(i,c) + end do + + do c=1,nV + do d=1,nV + aoooo(i,j,k,l) = aoooo(i,j,k,l) + OOVV(k,l,c,d)*tau(i,j,c,d) + end do + end do + + end do + end do + end do + end do + + bvvvv(:,:,:,:) = VVVV(:,:,:,:) + + do c=1,nV + do d=1,nV + do a=1,nV + do b=1,nV + + do k=1,nO + bvvvv(c,d,a,b) = bvvvv(c,d,a,b) - VOVV(a,k,c,d)*t1(k,b) + end do + + do k=1,nO + bvvvv(c,d,a,b) = bvvvv(c,d,a,b) + VOVV(b,k,c,d)*t1(k,a) + end do + + end do + end do + end do + end do + + hovvo(:,:,:,:) = OVVO(:,:,:,:) + + do i=1,nO + do c=1,nV + do a=1,nV + do k=1,nO + + do l=1,nO + hovvo(i,c,a,k) = hovvo(i,c,a,k) - OVOO(i,c,l,k)*t1(l,a) + end do + + do d=1,nV + hovvo(i,c,a,k) = hovvo(i,c,a,k) + OVVV(k,a,c,d)*t1(i,d) + end do + + do l=1,nO + do d=1,nV + hovvo(i,c,a,k) = hovvo(i,c,a,k) - OOVV(k,l,c,d)*tau(i,l,d,a) + end do + end do + + end do + end do + end do + end do + +end subroutine form_abh diff --git a/src/MCQC/form_delta_OOOVVV.f90 b/src/MCQC/form_delta_OOOVVV.f90 new file mode 100644 index 0000000..973b579 --- /dev/null +++ b/src/MCQC/form_delta_OOOVVV.f90 @@ -0,0 +1,37 @@ +subroutine form_delta_OOOVVV(nO,nV,eO,eV,delta) + +! Form energy denominator for CC + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + double precision,intent(in) :: eO(nO) + double precision,intent(in) :: eV(nV) + +! Local variables + + integer :: i,j,k,a,b,c + +! Output variables + + double precision,intent(out) :: delta(nO,nO,nO,nV,nV,nV) + + do i=1,nO + do j=1,nO + do k=1,nO + do a=1,nV + do b=1,nV + do c=1,nV + + delta(i,j,k,a,b,c) = eV(a) + eV(b) + eV(c) - eO(i) - eO(j) - eO(k) + + enddo + enddo + enddo + enddo + enddo + enddo + +end subroutine form_delta_OOOVVV diff --git a/src/MCQC/form_delta_OOVV.f90 b/src/MCQC/form_delta_OOVV.f90 new file mode 100644 index 0000000..15bb6d9 --- /dev/null +++ b/src/MCQC/form_delta_OOVV.f90 @@ -0,0 +1,33 @@ +subroutine form_delta_OOVV(nO,nV,eO,eV,delta) + +! Form energy denominator for CC + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + double precision,intent(in) :: eO(nO) + double precision,intent(in) :: eV(nV) + +! Local variables + + integer :: i,j,a,b + +! Output variables + + double precision,intent(out) :: delta(nO,nO,nV,nV) + + do i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + + delta(i,j,a,b) = eV(a) + eV(b) - eO(i) - eO(j) + + enddo + enddo + enddo + enddo + +end subroutine form_delta_OOVV diff --git a/src/MCQC/form_delta_OV.f90 b/src/MCQC/form_delta_OV.f90 new file mode 100644 index 0000000..e785a0c --- /dev/null +++ b/src/MCQC/form_delta_OV.f90 @@ -0,0 +1,27 @@ +subroutine form_delta_OV(nO,nV,eO,eV,delta) + +! Form energy denominator for CC + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + double precision,intent(in) :: eO(nO) + double precision,intent(in) :: eV(nV) + +! Local variables + + integer :: i,a + +! Output variables + + double precision,intent(out) :: delta(nO,nV) + + do i=1,nO + do a=1,nV + delta(i,a) = eV(a) - eO(i) + enddo + enddo + +end subroutine form_delta_OV diff --git a/src/MCQC/form_g.f90 b/src/MCQC/form_g.f90 new file mode 100644 index 0000000..ebb8427 --- /dev/null +++ b/src/MCQC/form_g.f90 @@ -0,0 +1,53 @@ +subroutine form_g(nO,nV,hvv,hoo,VOVV,OOOV,t1,gvv,goo) + +! Scuseria Eqs. (9), (10) + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + + double precision,intent(in) :: hvv(nV,nV) + double precision,intent(in) :: hoo(nO,nO) + + double precision,intent(in) :: VOVV(nV,nO,nV,nV) + double precision,intent(in) :: OOOV(nO,nO,nO,nV) + + double precision,intent(in) :: t1(nO,nV) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + +! Output variables + + double precision,intent(out) :: gvv(nV,nV) + double precision,intent(out) :: goo(nO,nO) + + gvv(:,:) = hvv(:,:) + + do c=1,nV + do a=1,nV + do k=1,nO + do d=1,nV + gvv(c,a) = gvv(c,a) + VOVV(a,k,c,d)*t1(k,d) + end do + end do + end do + end do + + goo(:,:) = hoo(:,:) + + do i=1,nO + do k=1,nO + do l=1,nO + do c=1,nV + goo(i,k) = goo(i,k) + OOOV(k,l,i,c)*t1(l,c) + end do + end do + end do + end do + +end subroutine form_g diff --git a/src/MCQC/form_h.f90 b/src/MCQC/form_h.f90 new file mode 100644 index 0000000..4bc16f1 --- /dev/null +++ b/src/MCQC/form_h.f90 @@ -0,0 +1,79 @@ +subroutine form_h(nO,nV,eO,eV,OOVV,t1,tau,hvv,hoo,hvo) + +! Scuseria Eqs. (5), (6) and (7) + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + + double precision,intent(in) :: eO(nO) + double precision,intent(in) :: eV(nV) + double precision,intent(in) :: OOVV(nO,nO,nV,nV) + + double precision,intent(in) :: t1(nO,nV) + double precision,intent(in) :: tau(nO,nO,nV,nV) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + +! Output variables + + double precision,intent(out) :: hvv(nV,nV) + double precision,intent(out) :: hoo(nO,nO) + double precision,intent(out) :: hvo(nV,nO) + + hvv(:,:) = 0d0 + + do b=1,nV + hvv(b,b) = eV(b) + do a=1,nV + do j=1,nO + do k=1,nO + do c=1,nV + + hvv(b,a) = hvv(b,a) - OOVV(j,k,b,c)*tau(j,k,a,c) + + end do + end do + end do + end do + end do + + hoo(:,:) = 0d0 + + do i=1,nO + hoo(i,i) = eO(i) + do j=1,nO + do k=1,nO + do b=1,nV + do c=1,nV + + hoo(i,j) = hoo(i,j) + OOVV(j,k,b,c)*tau(i,k,b,c) + + end do + end do + end do + end do + end do + + hvo(:,:) = 0d0 + + do b=1,nV + do j=1,nO + do k=1,nO + do c=1,nV + + hvo(b,j) = hvo(b,j) + OOVV(j,k,b,c)*t1(k,c) + + end do + end do + end do + end do + +! print*,'hvv',hvv + +end subroutine form_h diff --git a/src/MCQC/form_r1.f90 b/src/MCQC/form_r1.f90 new file mode 100644 index 0000000..ab1bd57 --- /dev/null +++ b/src/MCQC/form_r1.f90 @@ -0,0 +1,77 @@ +subroutine form_r1(nO,nV,OVVO,OVVV,OOOV,hvv,hoo,hvo,t1,t2,tau,r1) + +! Form tau in CCSD + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + + double precision,intent(in) :: OVVO(nO,nV,nV,nO) + double precision,intent(in) :: OVVV(nO,nV,nV,nV) + double precision,intent(in) :: OOOV(nO,nO,nO,nV) + + + double precision,intent(in) :: hvv(nV,nV) + double precision,intent(in) :: hoo(nO,nO) + double precision,intent(in) :: hvo(nV,nO) + + double precision,intent(in) :: t1(nO,nV) + double precision,intent(in) :: t2(nO,nO,nV,nV) + double precision,intent(in) :: tau(nO,nO,nV,nV) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + +! Output variables + + double precision,intent(out) :: r1(nO,nV) + + r1(:,:) = 0d0 + + do i=1,nO + do a=1,nV + + do b=1,nV + r1(i,a) = r1(i,a) + hvv(b,a)*t1(i,b) + end do + + do j=1,nO + r1(i,a) = r1(i,a) - hoo(i,j)*t1(j,a) + end do + + do j=1,nO + do b=1,nV + r1(i,a) = r1(i,a) + hvo(b,j)*(t2(i,j,a,b) + t1(i,b)*t1(j,a)) + end do + end do + + do j=1,nO + do b=1,nV + r1(i,a) = r1(i,a) + OVVO(i,b,a,j)*t1(j,b) + end do + end do + + do j=1,nO + do b=1,nV + do c=1,nV + r1(i,a) = r1(i,a) - OVVV(j,a,b,c)*tau(i,j,b,c) + end do + end do + end do + + do j=1,nO + do k=1,nO + do b=1,nV + r1(i,a) = r1(i,a) - OOOV(j,k,i,b)*tau(j,k,a,b) + end do + end do + end do + + end do + end do + +end subroutine form_r1 diff --git a/src/MCQC/form_r2.f90 b/src/MCQC/form_r2.f90 new file mode 100644 index 0000000..c57b8a8 --- /dev/null +++ b/src/MCQC/form_r2.f90 @@ -0,0 +1,139 @@ +subroutine form_r2(nO,nV,OOVV,OVOO,OVVV,OVVO,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau,r2) + +! Form tau in CCSD + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + + double precision,intent(in) :: OOVV(nO,nO,nV,nV) + double precision,intent(in) :: OVOO(nO,nV,nO,nO) + double precision,intent(in) :: OVVV(nO,nV,nV,nV) + double precision,intent(in) :: OVVO(nO,nV,nV,nO) + + double precision,intent(in) :: gvv(nV,nV) + double precision,intent(in) :: goo(nO,nO) + double precision,intent(in) :: aoooo(nO,nO,nO,nO) + double precision,intent(in) :: bvvvv(nV,nV,nV,nV) + double precision,intent(in) :: hovvo(nO,nV,nV,nO) + + double precision,intent(in) :: t1(nO,nV) + double precision,intent(in) :: t2(nO,nO,nV,nV) + double precision,intent(in) :: tau(nO,nO,nV,nV) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + +! Output variables + + double precision,intent(out) :: r2(nO,nO,nV,nV) + + r2(:,:,:,:) = OOVV(:,:,:,:) + + do i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + + do k=1,nO + do l=1,nO + r2(i,j,a,b) = r2(i,j,a,b) + aoooo(i,j,k,l)*tau(k,l,a,b) + end do + end do + + do c=1,nV + do d=1,nV + r2(i,j,a,b) = r2(i,j,a,b) + bvvvv(c,d,a,b)*tau(i,j,c,d) + end do + end do + + do c=1,nV + r2(i,j,a,b) = r2(i,j,a,b) + gvv(c,a)*t2(i,j,c,b) + end do + + do k=1,nO + r2(i,j,a,b) = r2(i,j,a,b) + OVOO(k,a,i,j)*t1(k,b) + end do + + do c=1,nV + r2(i,j,a,b) = r2(i,j,a,b) - gvv(c,b)*t2(i,j,c,a) + end do + + do k=1,nO + r2(i,j,a,b) = r2(i,j,a,b) - OVOO(k,b,i,j)*t1(k,a) + end do + + do k=1,nO + r2(i,j,a,b) = r2(i,j,a,b) - goo(i,k)*t2(k,j,a,b) + end do + + do c=1,nV + r2(i,j,a,b) = r2(i,j,a,b) + OVVV(j,c,b,a)*t1(i,c) + end do + + do k=1,nO + r2(i,j,a,b) = r2(i,j,a,b) + goo(j,k)*t2(k,i,a,b) + end do + + do c=1,nV + r2(i,j,a,b) = r2(i,j,a,b) - OVVV(i,c,b,a)*t1(j,c) + end do + + do k=1,nO + do c=1,nV + r2(i,j,a,b) = r2(i,j,a,b) + hovvo(i,c,a,k)*t2(j,k,b,c) + end do + end do + + do k=1,nO + do c=1,nV + r2(i,j,a,b) = r2(i,j,a,b) - OVVO(i,c,a,k)*t1(j,c)*t1(k,b) + end do + end do + + do k=1,nO + do c=1,nV + r2(i,j,a,b) = r2(i,j,a,b) - hovvo(j,c,a,k)*t2(i,k,b,c) + end do + end do + + do k=1,nO + do c=1,nV + r2(i,j,a,b) = r2(i,j,a,b) + OVVO(j,c,a,k)*t1(i,c)*t1(k,b) + end do + end do + + do k=1,nO + do c=1,nV + r2(i,j,a,b) = r2(i,j,a,b) - hovvo(i,c,b,k)*t2(j,k,a,c) + end do + end do + + do k=1,nO + do c=1,nV + r2(i,j,a,b) = r2(i,j,a,b) + OVVO(i,c,b,k)*t1(j,c)*t1(k,a) + end do + end do + + do k=1,nO + do c=1,nV + r2(i,j,a,b) = r2(i,j,a,b) + hovvo(j,c,b,k)*t2(i,k,a,c) + end do + end do + + do k=1,nO + do c=1,nV + r2(i,j,a,b) = r2(i,j,a,b) - OVVO(j,c,b,k)*t1(i,c)*t1(k,a) + end do + end do + + end do + end do + end do + end do + +end subroutine form_r2 diff --git a/src/MCQC/form_tau.f90 b/src/MCQC/form_tau.f90 new file mode 100644 index 0000000..d666383 --- /dev/null +++ b/src/MCQC/form_tau.f90 @@ -0,0 +1,34 @@ +subroutine form_tau(nO,nV,t1,t2,tau) + +! Form tau in CCSD + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + double precision,intent(in) :: t1(nO,nV) + double precision,intent(in) :: t2(nO,nO,nV,nV) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + +! Output variables + + double precision,intent(out) :: tau(nO,nO,nV,nV) + + do i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + + tau(i,j,a,b) = 0.5d0*t2(i,j,a,b) + t1(i,a)*t1(j,b) + + enddo + enddo + enddo + enddo + +end subroutine form_tau diff --git a/src/MCQC/form_ub.f90 b/src/MCQC/form_ub.f90 new file mode 100644 index 0000000..14f8d9c --- /dev/null +++ b/src/MCQC/form_ub.f90 @@ -0,0 +1,48 @@ +subroutine form_ub(nO,nV,OOVV,t1,ub) + +! Form 1st term in (T) correction + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + + double precision,intent(in) :: OOVV(nO,nO,nV,nV) + + double precision,intent(in) :: t1(nO,nV) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + +! Output variables + + double precision,intent(out) :: ub(nO,nO,nO,nV,nV,nV) + + do i=1,nO + do j=1,nO + do k=1,nO + do a=1,nV + do b=1,nV + do c=1,nV + + ub(i,j,k,a,b,c) = t1(i,a)*OOVV(j,k,b,c) & + + t1(i,b)*OOVV(j,k,c,a) & + + t1(i,c)*OOVV(j,k,a,b) & + + t1(j,a)*OOVV(k,i,b,c) & + + t1(j,b)*OOVV(k,i,c,a) & + + t1(j,c)*OOVV(k,i,a,b) & + + t1(k,a)*OOVV(i,j,b,c) & + + t1(k,b)*OOVV(i,j,c,a) & + + t1(k,c)*OOVV(i,j,a,b) + + end do + end do + end do + end do + end do + end do + +end subroutine form_ub diff --git a/src/MCQC/form_ubb.f90 b/src/MCQC/form_ubb.f90 new file mode 100644 index 0000000..5beaa6e --- /dev/null +++ b/src/MCQC/form_ubb.f90 @@ -0,0 +1,67 @@ +subroutine form_ubb(nO,nV,VVVO,VOOO,t2,ubb) + +! Form 2nd term in (T) correction + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + + double precision,intent(in) :: VVVO(nV,nV,nV,nO) + double precision,intent(in) :: VOOO(nV,nO,nO,nO) + + double precision,intent(in) :: t2(nO,nO,nV,nV) + +! Local variables + + integer :: i,j,k,l,m + integer :: a,b,c,d,e + +! Output variables + + double precision,intent(out) :: ubb(nO,nO,nO,nV,nV,nV) + + ubb(:,:,:,:,:,:) = 0d0 + + do i=1,nO + do j=1,nO + do k=1,nO + do a=1,nV + do b=1,nV + do c=1,nV + + do e=1,nV + ubb(i,j,k,a,b,c) = ubb(i,j,k,a,b,c) & + + t2(i,j,a,e)*VVVO(b,c,e,k) & + + t2(i,j,b,e)*VVVO(c,a,e,k) & + + t2(i,j,c,e)*VVVO(a,b,e,k) & + + t2(k,i,a,e)*VVVO(b,c,e,j) & + + t2(k,i,b,e)*VVVO(c,a,e,j) & + + t2(k,i,c,e)*VVVO(a,b,e,j) & + + t2(j,k,a,e)*VVVO(b,c,e,i) & + + t2(j,k,b,e)*VVVO(c,a,e,i) & + + t2(j,k,c,e)*VVVO(a,b,e,i) + end do + + do m=1,nO + ubb(i,j,k,a,b,c) = ubb(i,j,k,a,b,c) & + + t2(i,m,a,b)*VOOO(c,m,j,k) & + + t2(i,m,b,c)*VOOO(a,m,j,k) & + + t2(i,m,c,a)*VOOO(b,m,j,k) & + + t2(j,m,a,b)*VOOO(c,m,k,i) & + + t2(j,m,b,c)*VOOO(a,m,k,i) & + + t2(j,m,c,a)*VOOO(b,m,k,i) & + + t2(k,m,a,b)*VOOO(c,m,i,j) & + + t2(k,m,b,c)*VOOO(a,m,i,j) & + + t2(k,m,c,a)*VOOO(b,m,i,j) + end do + + end do + end do + end do + end do + end do + end do + +end subroutine form_ubb diff --git a/src/MCQC/form_v.f90 b/src/MCQC/form_v.f90 new file mode 100644 index 0000000..f589d31 --- /dev/null +++ b/src/MCQC/form_v.f90 @@ -0,0 +1,79 @@ +subroutine form_v(nO,nV,X1,X2,X3,X4,t2,v) + +! Form quadratic array in CCD + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + double precision,intent(in) :: t2(nO,nO,nV,nV) + double precision,intent(in) :: X1(nO,nO,nO,nO) + double precision,intent(in) :: X2(nV,nV) + double precision,intent(in) :: X3(nO,nO) + double precision,intent(in) :: X4(nO,nO,nV,nV) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + +! Output variables + + double precision,intent(out) :: v(nO,nO,nV,nV) + + v(:,:,:,:) = 0d0 + + do i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + do k=1,nO + do l=1,nO + v(i,j,a,b) = v(i,j,a,b) + 0.25d0*X1(k,l,i,j)*t2(k,l,a,b) + enddo + enddo + enddo + enddo + enddo + enddo + + do i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + do c=1,nV + v(i,j,a,b) = v(i,j,a,b) - 0.5d0*(X2(b,c)*t2(i,j,a,c) + X2(a,c)*t2(i,j,c,b)) + enddo + enddo + enddo + enddo + enddo + + do i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + do k=1,nO + v(i,j,a,b) = v(i,j,a,b) - 0.5d0*(X3(k,j)*t2(i,k,a,b) + X3(k,i)*t2(k,j,a,b)) + enddo + enddo + enddo + enddo + enddo + + do i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + do k=1,nO + do c=1,nV + v(i,j,a,b) = v(i,j,a,b) + (X4(i,k,a,c)*t2(j,k,b,c) + X4(i,k,b,c)*t2(k,j,a,c)) + enddo + enddo + enddo + enddo + enddo + enddo + +end subroutine form_v diff --git a/src/MCQC/spatial_to_spin_ERI.f90 b/src/MCQC/spatial_to_spin_ERI.f90 new file mode 100644 index 0000000..bdb2919 --- /dev/null +++ b/src/MCQC/spatial_to_spin_ERI.f90 @@ -0,0 +1,35 @@ +subroutine spatial_to_spin_ERI(nBas,ERI,nBas2,sERI) + +! Convert ERIs from spatial to spin orbitals + + implicit none + +! Input variables + + integer,intent(in) :: nBas,nBas2 + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: p,q,r,s + double precision,external :: Kronecker_delta + +! Output variables + + double precision,intent(out) :: sERI(nBas2,nBas2,nBas2,nBas2) + + do p=1,nBas2 + do q=1,nBas2 + do r=1,nBas2 + do s=1,nBas2 + + sERI(p,q,r,s) = Kronecker_delta(mod(p,2),mod(r,2)) & + * Kronecker_delta(mod(q,2),mod(s,2)) & + * ERI((p+1)/2,(q+1)/2,(r+1)/2,(s+1)/2) + + enddo + enddo + enddo + enddo + +end subroutine spatial_to_spin_ERI diff --git a/src/MCQC/spatial_to_spin_MO_energy.f90 b/src/MCQC/spatial_to_spin_MO_energy.f90 new file mode 100644 index 0000000..688dc1d --- /dev/null +++ b/src/MCQC/spatial_to_spin_MO_energy.f90 @@ -0,0 +1,29 @@ +subroutine spatial_to_spin_MO_energy(nBas,e,nBas2,se) + +! Convert ERIs from spatial to spin orbitals + + implicit none + +! Input variables + + integer,intent(in) :: nBas,nBas2 + double precision,intent(in) :: e(nBas) + +! Local variables + + integer :: p + +! Output variables + + double precision,intent(out) :: se(nBas2) + + do p=1,nBas2 + + se(p) = e((p+1)/2) + + enddo + +! print*,'MO energies in spinorbital basis' +! call matout(nBas2,1,se) + +end subroutine spatial_to_spin_MO_energy