From c69525722a0c150f3f7f067c4b6140d9e54f4348 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 9 Sep 2019 16:51:15 +0200 Subject: [PATCH] fork CC from mveril and pfloos --- devel/cc/CC.irp.f | 11 ++ devel/cc/CCSD.irp.f | 211 +++++++++++++++++++++++ devel/cc/CCSDT.irp.f | 35 ++++ devel/cc/EZFIO.cfg | 22 +++ devel/cc/Kronecker_delta.irp.f | 22 +++ devel/cc/MP2.irp.f | 71 ++++++++ devel/cc/NEED | 2 + devel/cc/README.rst | 4 + devel/cc/amplitude_guess.irp.f | 19 ++ devel/cc/antisymmetrize_ERI.irp.f | 46 +++++ devel/cc/denominators.irp.f | 70 ++++++++ devel/cc/dimensions.irp.f | 20 +++ devel/cc/form_T.irp.f | 45 +++++ devel/cc/form_abh.irp.f | 97 +++++++++++ devel/cc/form_g.irp.f | 50 ++++++ devel/cc/form_h.irp.f | 75 ++++++++ devel/cc/form_r1.irp.f | 72 ++++++++ devel/cc/form_r2.irp.f | 134 ++++++++++++++ devel/cc/form_tau.irp.f | 34 ++++ devel/cc/form_ub.irp.f | 46 +++++ devel/cc/form_ubb.irp.f | 64 +++++++ devel/cc/integral_batches.irp.f | 79 +++++++++ devel/cc/save_energy.irp.f | 8 + devel/cc/spatial_to_spin_ERI.irp.f | 52 ++++++ devel/cc/spatial_to_spin_MO_energy.irp.f | 13 ++ 25 files changed, 1302 insertions(+) create mode 100644 devel/cc/CC.irp.f create mode 100644 devel/cc/CCSD.irp.f create mode 100644 devel/cc/CCSDT.irp.f create mode 100644 devel/cc/EZFIO.cfg create mode 100644 devel/cc/Kronecker_delta.irp.f create mode 100644 devel/cc/MP2.irp.f create mode 100644 devel/cc/NEED create mode 100644 devel/cc/README.rst create mode 100644 devel/cc/amplitude_guess.irp.f create mode 100644 devel/cc/antisymmetrize_ERI.irp.f create mode 100644 devel/cc/denominators.irp.f create mode 100644 devel/cc/dimensions.irp.f create mode 100644 devel/cc/form_T.irp.f create mode 100644 devel/cc/form_abh.irp.f create mode 100644 devel/cc/form_g.irp.f create mode 100644 devel/cc/form_h.irp.f create mode 100644 devel/cc/form_r1.irp.f create mode 100644 devel/cc/form_r2.irp.f create mode 100644 devel/cc/form_tau.irp.f create mode 100644 devel/cc/form_ub.irp.f create mode 100644 devel/cc/form_ubb.irp.f create mode 100644 devel/cc/integral_batches.irp.f create mode 100644 devel/cc/save_energy.irp.f create mode 100644 devel/cc/spatial_to_spin_ERI.irp.f create mode 100644 devel/cc/spatial_to_spin_MO_energy.irp.f diff --git a/devel/cc/CC.irp.f b/devel/cc/CC.irp.f new file mode 100644 index 0000000..0710677 --- /dev/null +++ b/devel/cc/CC.irp.f @@ -0,0 +1,11 @@ +program CC + implicit none + BEGIN_DOC + ! CC + ! Coupled cluster + ! + ! This program perform coupled cluster calculation + ! CCSD or CCSD(T) + END_DOC + call CCSD +end diff --git a/devel/cc/CCSD.irp.f b/devel/cc/CCSD.irp.f new file mode 100644 index 0000000..4c57b22 --- /dev/null +++ b/devel/cc/CCSD.irp.f @@ -0,0 +1,211 @@ +subroutine CCSD + +! CCSD module + + implicit none + +! Input variables (provided by IRP) + + integer :: maxscf + double precision :: thresh + + logical :: doCCSDT + integer :: nBas,nEl + double precision :: ERHF + +! Local variables + + integer :: p,q,r,s + 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 :: get_two_e_integral,u_dot_v + + 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(*,*) + +! IRP init + provide cc_mode + maxSCF=cc_n_it_max + thresh=cc_thresh + + doCCSDT = trim(cc_mode)=='CCSD(T)' + + nBas=mo_num + nEl=elec_num + ERHF=hf_energy + +! Spatial to spin orbitals + + nBas2 = spin_mo_num + + +! Define occupied and virtual spaces + + nO = spin_occ_num + nV = spin_vir_num + +! Guess amplitudes + + allocate(t1(nO,nV),t2(nO,nO,nV,nV),tau(nO,nO,nV,nV)) + + t1(:,:) = t1_guess(:,:) + t2(:,:,:,:) = t2_guess(:,:,:,:) + + call form_tau(nO,nV,t1,t2,tau) + + EcMP2 = 0.5d0*u_dot_v(pack(OOVV,.true.),pack(tau,.true.),size(OOVV)) + 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,t1,tau,hvv,hoo,hvo) + +! Scuseria Eqs. (9), (10), (11), (12) and (13) + + call form_g(nO,nV,hvv,hoo,t1,gvv,goo) + + call form_abh(nO,nV,t1,tau,aoooo,bvvvv,hovvo) + +! Compute residuals + + call form_r1(nO,nV,hvv,hoo,hvo,t1,t2,tau,r1) + + call form_r2(nO,nV,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*u_dot_v(pack(OOVV,.true.),pack(tau,.true.),size(OOVV)) + +! 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,'|',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, & + gvv,goo, & + aoooo,bvvvv,hovvo, & + tau, & + r1,r2) + +!------------------------------------------------------------------------ +! (T) correction +!------------------------------------------------------------------------ + if(doCCSDT) then + write(*,*) "Starting (T) calculation" +! call cpu_time(start_CCSDT)V + call CCSDT(nO,nV,t1,t2,EcCCT) +! call cpu_time(end_CCSDT) + call write_time(6) + +! 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(*,*) + + call save_energy(ECCSD + EcCCT) + + else + + call save_energy(ECCSD) + + end if + +end subroutine CCSD diff --git a/devel/cc/CCSDT.irp.f b/devel/cc/CCSDT.irp.f new file mode 100644 index 0000000..700848d --- /dev/null +++ b/devel/cc/CCSDT.irp.f @@ -0,0 +1,35 @@ +subroutine CCSDT(nO,nV,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) :: t1(nO,nV) + double precision,intent(in) :: t2(nO,nO,nV,nV) + +! Local variables + + double precision,allocatable :: ub(:,:,:,:,:,:) + double precision,allocatable :: ubb(:,:,:,:,:,:) + +! Output variables + + double precision,intent(out) :: EcCCT + +! Memory allocation + + allocate(ub(nO,nO,nO,nV,nV,nV),ubb(nO,nO,nO,nV,nV,nV)) + +! Form CCSD(T) quantities + + call form_ub(nO,nV,t1,ub) + + call form_ubb(nO,nV,t2,ubb) + + call form_T(nO,nV,ub,ubb,EcCCT) + +end subroutine CCSDT diff --git a/devel/cc/EZFIO.cfg b/devel/cc/EZFIO.cfg new file mode 100644 index 0000000..32891b4 --- /dev/null +++ b/devel/cc/EZFIO.cfg @@ -0,0 +1,22 @@ +[energy] +type: Threshold +doc: Energy CC +interface: ezfio + +[cc_mode] +type: character*(32) +doc: mode CCD CCSD CCSD(T). +interface: ezfio,provider,ocaml +default:CCSD + +[cc_thresh] +type: Threshold +doc: Threshold on the convergence of the Coupled-Cluster residual. +interface: ezfio,provider,ocaml +default: 1.e-10 + +[cc_n_it_max] +type: Strictly_positive_int +doc: Maximum number of Coupled-Cluster iterations. +interface: ezfio,provider,ocaml +default: 500 diff --git a/devel/cc/Kronecker_delta.irp.f b/devel/cc/Kronecker_delta.irp.f new file mode 100644 index 0000000..beac77d --- /dev/null +++ b/devel/cc/Kronecker_delta.irp.f @@ -0,0 +1,22 @@ +!------------------------------------------------------------------------ +function Kronecker_delta(i,j) result(delta) + +! Kronecker Delta + + implicit none + +! Input variables + + integer,intent(in) :: i,j + +! Output variables + + double precision :: delta + + if(i == j) then + delta = 1d0 + else + delta = 0d0 + endif + +end function Kronecker_delta diff --git a/devel/cc/MP2.irp.f b/devel/cc/MP2.irp.f new file mode 100644 index 0000000..3e38149 --- /dev/null +++ b/devel/cc/MP2.irp.f @@ -0,0 +1,71 @@ +subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2) + +! Perform third-order Moller-Plesset calculation + + implicit none + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nV,nR + double precision,intent(in) :: ENuc,EHF + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),e(nBas) + +! Local variables + + integer :: i,j,a,b + double precision :: eps,E2a,E2b + +! Output variables + + double precision,intent(out) :: EcMP2(3) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| Moller-Plesset second-order calculation |' + write(*,*)'************************************************' + write(*,*) + +! Compute MP2 energy + + E2a = 0d0 + E2b = 0d0 + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + + eps = e(i) + e(j) - e(a) - e(b) + +! Second-order ring diagram + + E2a = E2a + ERI(i,j,a,b)*ERI(i,j,a,b)/eps + +! Second-order exchange diagram + + E2b = E2b + ERI(i,j,a,b)*ERI(i,j,b,a)/eps + + enddo + enddo + enddo + enddo + + EcMP2(2) = 2d0*E2a + EcMP2(3) = -E2b + EcMP2(1) = EcMP2(2) + EcMP2(3) + + write(*,*) + write(*,'(A32)') '-----------------------' + write(*,'(A32)') ' MP2 calculation ' + write(*,'(A32)') '-----------------------' + write(*,'(A32,1X,F16.10)') ' MP2 correlation energy',EcMP2(1) + write(*,'(A32,1X,F16.10)') ' Direct part ',EcMP2(2) + write(*,'(A32,1X,F16.10)') ' Exchange part ',EcMP2(3) + write(*,'(A32)') '-----------------------' + write(*,'(A32,1X,F16.10)') ' MP2 electronic energy',EHF + EcMP2(1) + write(*,'(A32,1X,F16.10)') ' MP2 total energy',ENuc + EHF + EcMP2(1) + write(*,'(A32)') '-----------------------' + write(*,*) + +end subroutine MP2 diff --git a/devel/cc/NEED b/devel/cc/NEED new file mode 100644 index 0000000..5e2e753 --- /dev/null +++ b/devel/cc/NEED @@ -0,0 +1,2 @@ +hartree_fock +mo_two_e_ints diff --git a/devel/cc/README.rst b/devel/cc/README.rst new file mode 100644 index 0000000..eba5f26 --- /dev/null +++ b/devel/cc/README.rst @@ -0,0 +1,4 @@ +== +CC +== + diff --git a/devel/cc/amplitude_guess.irp.f b/devel/cc/amplitude_guess.irp.f new file mode 100644 index 0000000..678e298 --- /dev/null +++ b/devel/cc/amplitude_guess.irp.f @@ -0,0 +1,19 @@ +BEGIN_PROVIDER [ double precision, t1_guess, (spin_occ_num,spin_vir_num) ] + implicit none + BEGIN_DOC + ! Guess amplitudes for single excitations + END_DOC + + t1_guess(:,:) = 0d0 +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, t2_guess, (spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num) ] + implicit none + BEGIN_DOC + ! Guess amplitudes for double excitations + END_DOC + + t2_guess(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + +END_PROVIDER diff --git a/devel/cc/antisymmetrize_ERI.irp.f b/devel/cc/antisymmetrize_ERI.irp.f new file mode 100644 index 0000000..96e196d --- /dev/null +++ b/devel/cc/antisymmetrize_ERI.irp.f @@ -0,0 +1,46 @@ +subroutine antisymmetrize_ERI(ispin,nBas,ERI,db_ERI) + +! Antisymmetrize ERIs + + implicit none + +! Input variables + + integer,intent(in) :: ispin,nBas + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: i,j,k,l + +! Output variables + + double precision,intent(out) :: db_ERI(nBas,nBas,nBas,nBas) + + if(ispin == 1) then + + do l=1,nBas + do k=1,nBas + do j=1,nBas + do i=1,nBas + db_ERI(i,j,k,l) = 2d0*ERI(i,j,k,l) - ERI(i,j,l,k) + enddo + enddo + enddo + enddo + + elseif(ispin == 2) then + + do l=1,nBas + do k=1,nBas + do j=1,nBas + do i=1,nBas + db_ERI(i,j,k,l) = ERI(i,j,k,l) - ERI(i,j,l,k) + enddo + enddo + enddo + enddo + + endif + +end subroutine antisymmetrize_ERI diff --git a/devel/cc/denominators.irp.f b/devel/cc/denominators.irp.f new file mode 100644 index 0000000..b191940 --- /dev/null +++ b/devel/cc/denominators.irp.f @@ -0,0 +1,70 @@ + BEGIN_PROVIDER [ double precision, eO, (1:spin_occ_num) ] +&BEGIN_PROVIDER [ double precision, eV, (1:spin_vir_num) ] + implicit none + BEGIN_DOC + ! Eigenvalues of the spin-orbitals + END_DOC + eO(:) = spin_fock_matrix_diag_mo(1:spin_occ_num) + eV(:) = spin_fock_matrix_diag_mo(spin_occ_num+1:spin_mo_num) +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, delta_OV, (spin_occ_num, spin_vir_num) ] + implicit none + BEGIN_DOC + ! Energy denominator for CC + END_DOC + integer :: i, a + + do a=1,spin_vir_num + do i=1,spin_occ_num + delta_OV(i,a) = eV(a) - eO(i) + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, delta_OOVV, (spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num) ] + implicit none + BEGIN_DOC + ! Energy denominator for CC + END_DOC + integer :: i,j,a,b + + do b=1,spin_vir_num + do a=1,spin_vir_num + do j=1,spin_occ_num + do i=1,spin_occ_num + delta_OOVV(i,j,a,b) = eV(a) - eO(i) + eV(b) - eO(j) + enddo + enddo + enddo + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, delta_OOOVVV, (spin_occ_num,spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num,spin_vir_num) ] + implicit none + BEGIN_DOC + ! Energy denominator for CC + END_DOC + integer :: i,j,k,a,b,c + + do c=1,spin_vir_num + do b=1,spin_vir_num + do a=1,spin_vir_num + do k=1,spin_occ_num + do j=1,spin_occ_num + do i=1,spin_occ_num + delta_OOOVVV(i,j,k,a,b,c) = eV(a) - eO(i) + eV(b) - eO(j) + eV(c) - eO(k) + enddo + enddo + enddo + enddo + enddo + enddo +END_PROVIDER + + diff --git a/devel/cc/dimensions.irp.f b/devel/cc/dimensions.irp.f new file mode 100644 index 0000000..24cd58c --- /dev/null +++ b/devel/cc/dimensions.irp.f @@ -0,0 +1,20 @@ +BEGIN_PROVIDER [ integer, spin_mo_num ] + implicit none + BEGIN_DOC + ! Number of spin-orbitals + END_DOC + spin_mo_num = mo_num * 2 +END_PROVIDER + + + BEGIN_PROVIDER [ integer, spin_occ_num ] +&BEGIN_PROVIDER [ integer, spin_vir_num ] + implicit none + BEGIN_DOC + ! Number of occupied/virtual spin-orbitals + END_DOC + spin_occ_num = elec_num + spin_vir_num = spin_mo_num - spin_occ_num +END_PROVIDER + + diff --git a/devel/cc/form_T.irp.f b/devel/cc/form_T.irp.f new file mode 100644 index 0000000..0ace502 --- /dev/null +++ b/devel/cc/form_T.irp.f @@ -0,0 +1,45 @@ +subroutine form_T(nO,nV,ub,ubb,EcCCT) + +! Compute (T) correction + + implicit none + +! Input variables + + integer,intent(in) :: nO,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 c=1,nV + do b=1,nV + do a=1,nV + do k=1,nO + do j=1,nO + do i=1,nO + + 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/devel/cc/form_abh.irp.f b/devel/cc/form_abh.irp.f new file mode 100644 index 0000000..85f09c9 --- /dev/null +++ b/devel/cc/form_abh.irp.f @@ -0,0 +1,97 @@ +subroutine form_abh(nO,nV,t1,tau,aoooo,bvvvv,hovvo) + +! Scuseria Eqs. (11),(12) and (13) + + implicit none + +! Input variables + + integer,intent(in) :: nO,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 l=1,nO + do k=1,nO + do j=1,nO + do i=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 d=1,nV + do c=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 b=1,nV + do a=1,nV + do d=1,nV + do c=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 k=1,nO + do a=1,nV + do c=1,nV + do i=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 d=1,nV + do l=1,nO + 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/devel/cc/form_g.irp.f b/devel/cc/form_g.irp.f new file mode 100644 index 0000000..b78d11a --- /dev/null +++ b/devel/cc/form_g.irp.f @@ -0,0 +1,50 @@ +subroutine form_g(nO,nV,hvv,hoo,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) :: 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 a=1,nV + do c=1,nV + do d=1,nV + do k=1,nO + 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 k=1,nO + do i=1,nO + do c=1,nV + do l=1,nO + 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/devel/cc/form_h.irp.f b/devel/cc/form_h.irp.f new file mode 100644 index 0000000..fde9c4c --- /dev/null +++ b/devel/cc/form_h.irp.f @@ -0,0 +1,75 @@ +subroutine form_h(nO,nV,t1,tau,hvv,hoo,hvo) + +! Scuseria Eqs. (5), (6) and (7) + + implicit none + +! Input variables + + integer,intent(in) :: nO,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 j=1,nO + do b=1,nV + do c=1,nV + do k=1,nO + + 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/devel/cc/form_r1.irp.f b/devel/cc/form_r1.irp.f new file mode 100644 index 0000000..4e00050 --- /dev/null +++ b/devel/cc/form_r1.irp.f @@ -0,0 +1,72 @@ +subroutine form_r1(nO,nV,hvv,hoo,hvo,t1,t2,tau,r1) + +! Form tau in CCSD + + 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) :: 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 a=1,nV + do i=1,nO + + 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 b=1,nV + do j=1,nO + r1(i,a) = r1(i,a) + OVVO(i,b,a,j)*t1(j,b) + end do + end do + + do c=1,nV + do b=1,nV + do j=1,nO + r1(i,a) = r1(i,a) - OVVV(j,a,b,c)*tau(i,j,b,c) + end do + end do + end do + + do b=1,nV + do k=1,nO + do j=1,nO + 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/devel/cc/form_r2.irp.f b/devel/cc/form_r2.irp.f new file mode 100644 index 0000000..f869300 --- /dev/null +++ b/devel/cc/form_r2.irp.f @@ -0,0 +1,134 @@ +subroutine form_r2(nO,nV,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) :: 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 b=1,nV + do a=1,nV + do j=1,nO + do i=1,nO + + 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 d=1,nV + do c=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 c=1,nV + do k=1,nO + 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 c=1,nV + do k=1,nO + 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 c=1,nV + do k=1,nO + 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 c=1,nV + do k=1,nO + 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 c=1,nV + do k=1,nO + 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/devel/cc/form_tau.irp.f b/devel/cc/form_tau.irp.f new file mode 100644 index 0000000..c7f1cd8 --- /dev/null +++ b/devel/cc/form_tau.irp.f @@ -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 b=1,nV + do a=1,nV + do j=1,nO + do i=1,nO + + 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/devel/cc/form_ub.irp.f b/devel/cc/form_ub.irp.f new file mode 100644 index 0000000..008bfe9 --- /dev/null +++ b/devel/cc/form_ub.irp.f @@ -0,0 +1,46 @@ +subroutine form_ub(nO,nV,t1,ub) + +! Form 1st term in (T) correction + + implicit none + +! Input variables + + integer,intent(in) :: 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) :: ub(nO,nO,nO,nV,nV,nV) + + do c=1,nV + do b=1,nV + do a=1,nV + do k=1,nO + do j=1,nO + do i=1,nO + + 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/devel/cc/form_ubb.irp.f b/devel/cc/form_ubb.irp.f new file mode 100644 index 0000000..985dc81 --- /dev/null +++ b/devel/cc/form_ubb.irp.f @@ -0,0 +1,64 @@ +subroutine form_ubb(nO,nV,t2,ubb) + +! Form 2nd term in (T) correction + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + + 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 c=1,nV + do b=1,nV + do a=1,nV + do k=1,nO + do j=1,nO + do i=1,nO + + 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/devel/cc/integral_batches.irp.f b/devel/cc/integral_batches.irp.f new file mode 100644 index 0000000..50dc42d --- /dev/null +++ b/devel/cc/integral_batches.irp.f @@ -0,0 +1,79 @@ +BEGIN_PROVIDER [ double precision, OOOO, (spin_occ_num,spin_occ_num,spin_occ_num,spin_occ_num) ] + implicit none + BEGIN_DOC + END_DOC + OOOO(:,:,:,:) = dbERI( 1:spin_occ_num , 1:spin_occ_num , 1:spin_occ_num , 1:spin_occ_num ) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, OOOV, (spin_occ_num,spin_occ_num,spin_occ_num,spin_vir_num) ] + implicit none + BEGIN_DOC + END_DOC + OOOV(:,:,:,:) = dbERI( 1:spin_occ_num , 1:spin_occ_num , 1:spin_occ_num ,spin_occ_num+1:spin_mo_num) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, OVOO, (spin_occ_num,spin_vir_num,spin_occ_num,spin_occ_num) ] + implicit none + BEGIN_DOC + END_DOC + OVOO(:,:,:,:) = dbERI( 1:spin_occ_num ,spin_occ_num+1:spin_mo_num, 1:spin_occ_num , 1:spin_occ_num ) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, VOOO, (spin_vir_num,spin_occ_num,spin_occ_num,spin_occ_num) ] + implicit none + BEGIN_DOC + END_DOC + VOOO(:,:,:,:) = dbERI(spin_occ_num+1:spin_mo_num, 1:spin_occ_num , 1:spin_occ_num , 1:spin_occ_num ) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, OOVV, (spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num) ] + implicit none + BEGIN_DOC + END_DOC + OOVV(:,:,:,:) = dbERI( 1:spin_occ_num , 1:spin_occ_num ,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, OVVO, (spin_occ_num,spin_vir_num,spin_vir_num,spin_occ_num) ] + implicit none + BEGIN_DOC + END_DOC + OVVO(:,:,:,:) = dbERI( 1:spin_occ_num ,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num, 1:spin_occ_num ) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, OVVV, (spin_occ_num,spin_vir_num,spin_vir_num,spin_vir_num) ] + implicit none + BEGIN_DOC + END_DOC + OVVV(:,:,:,:) = dbERI( 1:spin_occ_num ,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, VOVV, (spin_vir_num,spin_occ_num,spin_vir_num,spin_vir_num) ] + implicit none + BEGIN_DOC + END_DOC + VOVV(:,:,:,:) = dbERI(spin_occ_num+1:spin_mo_num, 1:spin_occ_num ,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, VVVO, (spin_vir_num,spin_vir_num,spin_vir_num,spin_occ_num) ] + implicit none + BEGIN_DOC + END_DOC + VVVO(:,:,:,:) = dbERI(spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num, 1:spin_occ_num ) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, VVVV, (spin_vir_num,spin_vir_num,spin_vir_num,spin_vir_num) ] + implicit none + BEGIN_DOC + END_DOC + VVVV(:,:,:,:) = dbERI(spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num,spin_occ_num+1:spin_mo_num) +END_PROVIDER + diff --git a/devel/cc/save_energy.irp.f b/devel/cc/save_energy.irp.f new file mode 100644 index 0000000..82aad00 --- /dev/null +++ b/devel/cc/save_energy.irp.f @@ -0,0 +1,8 @@ +subroutine save_energy(E) + implicit none + BEGIN_DOC +! Saves the energy in |EZFIO|. + END_DOC + double precision, intent(in) :: E + call ezfio_set_cc_energy(E) +end diff --git a/devel/cc/spatial_to_spin_ERI.irp.f b/devel/cc/spatial_to_spin_ERI.irp.f new file mode 100644 index 0000000..e09c0da --- /dev/null +++ b/devel/cc/spatial_to_spin_ERI.irp.f @@ -0,0 +1,52 @@ +BEGIN_PROVIDER [ double precision, dbERI, (spin_mo_num,spin_mo_num,spin_mo_num,spin_mo_num) ] + implicit none + BEGIN_DOC + ! Anti-symmetrized Electron repulsion integrals in spin-orbital basis + END_DOC + +! Local variables + + integer :: p,q,r,s + double precision,external :: Kronecker_delta + +! Output variables + + double precision, external :: get_two_e_integral + + double precision :: pqrs, pqsr + + PROVIDE mo_two_e_integrals_in_map + + do s=1,spin_mo_num + do r=1,spin_mo_num + do q=1,spin_mo_num + do p=1,spin_mo_num + + pqrs = Kronecker_delta(mod(p,2),mod(r,2)) & + * Kronecker_delta(mod(q,2),mod(s,2)) & + * get_two_e_integral( & + (p+1)/2, & + (q+1)/2, & + (r+1)/2, & + (s+1)/2, & + mo_two_e_integrals_in_map) + + pqsr = Kronecker_delta(mod(p,2),mod(s,2)) & + * Kronecker_delta(mod(q,2),mod(r,2)) & + * get_two_e_integral( & + (p+1)/2, & + (q+1)/2, & + (s+1)/2, & + (r+1)/2, & + mo_two_e_integrals_in_map) + + dbERI(p,q,r,s) = pqrs - pqsr + enddo + enddo + enddo + enddo + + +END_PROVIDER + + diff --git a/devel/cc/spatial_to_spin_MO_energy.irp.f b/devel/cc/spatial_to_spin_MO_energy.irp.f new file mode 100644 index 0000000..0ce0576 --- /dev/null +++ b/devel/cc/spatial_to_spin_MO_energy.irp.f @@ -0,0 +1,13 @@ +BEGIN_PROVIDER [ double precision, spin_fock_matrix_diag_mo, (spin_mo_num) ] + implicit none + BEGIN_DOC + ! Diagonal of the Fock matrix in the spin-orbital basis + END_DOC + integer :: p + + do p=1,spin_mo_num + spin_fock_matrix_diag_mo(p) = fock_matrix_diag_mo((p+1)/2) + enddo + +END_PROVIDER +