From baa350efb16dcd542da874968399fb972c2db2a5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 11 Sep 2019 17:09:30 +0200 Subject: [PATCH] General CCSD --- devel/cc/.gitignore | 5 + devel/cc/CCSD.irp.f | 75 +++++++------ devel/cc/CCSD_Ec_nc.irp.f | 54 +++++++++ devel/cc/NEED | 1 + devel/cc/amplitude_guess.irp.f | 3 +- devel/cc/denominators.irp.f | 7 +- devel/cc/form_cF_nc.irp.f | 102 +++++++++++++++++ devel/cc/form_cW_nc.irp.f | 105 ++++++++++++++++++ devel/cc/form_r1_nc.irp.f | 77 +++++++++++++ devel/cc/form_r2_nc.irp.f | 102 +++++++++++++++++ devel/cc/form_tau_nc.irp.f | 34 ++++++ devel/cc/form_taus_nc.irp.f | 34 ++++++ devel/cc/integral_batches.irp.f | 122 +++++++++++++++++++-- devel/cc/spatial_to_spin_MO.irp.f | 120 ++++++++++++++++++++ devel/cc/spatial_to_spin_MO_energy.irp.f | 13 --- stable/amplitudes/extract_amplitudes.irp.f | 31 ++++-- 16 files changed, 810 insertions(+), 75 deletions(-) create mode 100644 devel/cc/.gitignore create mode 100644 devel/cc/CCSD_Ec_nc.irp.f create mode 100644 devel/cc/form_cF_nc.irp.f create mode 100644 devel/cc/form_cW_nc.irp.f create mode 100644 devel/cc/form_r1_nc.irp.f create mode 100644 devel/cc/form_r2_nc.irp.f create mode 100644 devel/cc/form_tau_nc.irp.f create mode 100644 devel/cc/form_taus_nc.irp.f create mode 100644 devel/cc/spatial_to_spin_MO.irp.f delete mode 100644 devel/cc/spatial_to_spin_MO_energy.irp.f diff --git a/devel/cc/.gitignore b/devel/cc/.gitignore new file mode 100644 index 0000000..7ac9fbf --- /dev/null +++ b/devel/cc/.gitignore @@ -0,0 +1,5 @@ +IRPF90_temp/ +IRPF90_man/ +irpf90.make +irpf90_entities +tags \ No newline at end of file diff --git a/devel/cc/CCSD.irp.f b/devel/cc/CCSD.irp.f index 4c57b22..20fde0c 100644 --- a/devel/cc/CCSD.irp.f +++ b/devel/cc/CCSD.irp.f @@ -27,14 +27,13 @@ subroutine CCSD 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 :: cFvv(:,:) + double precision,allocatable :: cFoo(:,:) + double precision,allocatable :: cFov(:,:) + + double precision,allocatable :: cWoooo(:,:,:,:) + double precision,allocatable :: cWvvvv(:,:,:,:) + double precision,allocatable :: cWovvo(:,:,:,:) double precision,allocatable :: r1(:,:) double precision,allocatable :: r2(:,:,:,:) @@ -42,6 +41,7 @@ subroutine CCSD double precision,allocatable :: t1(:,:) double precision,allocatable :: t2(:,:,:,:) double precision,allocatable :: tau(:,:,:,:) + double precision,allocatable :: taus(:,:,:,:) ! Hello world @@ -74,26 +74,27 @@ subroutine CCSD ! Guess amplitudes - allocate(t1(nO,nV),t2(nO,nO,nV,nV),tau(nO,nO,nV,nV)) + allocate(t1(nO,nV),t2(nO,nO,nV,nV),tau(nO,nO,nV,nV),taus(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), & + allocate(cFvv(nV,nV),cFoo(nO,nO),cFov(nO,nV), & + cWoooo(nO,nO,nO,nO),cWvvvv(nV,nV,nV,nV),cWovvo(nO,nV,nV,nO), & r1(nO,nV),r2(nO,nO,nV,nV)) Conv = 1d0 nSCF = 0 + call form_taus_nc(nO,nV,t1,t2,taus) + call form_tau_nc (nO,nV,t1,t2,tau) + + EcMP2 = 0.25d0*u_dot_v(OOVV,tau,size(OOVV)) + write(*,'(1X,A10,1X,F10.6)') 'Ec(MP2) = ',EcMP2 + write(*,'(1X,A10,1X,F10.6)') 'E (MP2) = ',EcMP2 + ERHF + !------------------------------------------------------------------------ ! Main SCF loop !------------------------------------------------------------------------ @@ -107,25 +108,26 @@ subroutine CCSD do while(Conv > thresh .and. nSCF < maxSCF) -! Increment +! Excrement nSCF = nSCF + 1 -! Scuseria Eqs. (5), (6) and (7) + call form_cf_nc (nO,nV,t1,taus, & + spin_fock_matrix_mo_oo, & + spin_fock_matrix_mo_ov, & + spin_fock_matrix_mo_vv, & + cFoo,cFov,cFvv) - 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) + call form_cw_nc (nO,nV,t1,t2,tau, & + cWoooo,cWovvo,cWvvvv) ! Compute residuals - call form_r1(nO,nV,hvv,hoo,hvo,t1,t2,tau,r1) + call form_r1_nc(nO,nV,t1,t2,spin_fock_matrix_mo_ov, & + cFoo,cFov,cFvv,r1) - call form_r2(nO,nV,gvv,goo,aoooo,bvvvv,hovvo,t1,t2,tau,r2) + call form_r2_nc(nO,nV,t1,t2,tau,cFoo,cFov,cFvv, & + cWoooo,cWvvvv,cWovvo,r2) ! Check convergence @@ -136,11 +138,12 @@ subroutine CCSD t1(:,:) = t1(:,:) - r1(:,:) /delta_OV (:,:) t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:) - call form_tau(nO,nV,t1,t2,tau) - + call form_taus_nc(nO,nV,t1,t2,taus) + call form_tau_nc (nO,nV,t1,t2,tau) + ! Compute correlation energy - EcCCSD = 0.5d0*u_dot_v(pack(OOVV,.true.),pack(tau,.true.),size(OOVV)) + call CCSD_Ec_nc(nO,nV,t1,t2,spin_fock_matrix_mo_ov,EcCCSD) ! Dump results @@ -171,11 +174,11 @@ subroutine CCSD ! Deallocate memory - deallocate(hvv,hoo,hvo, & - gvv,goo, & - aoooo,bvvvv,hovvo, & - tau, & - r1,r2) + deallocate( & + cFvv,cFoo,cFov, & + cWoooo,cWvvvv,cWovvo, & + tau,taus, & + r1,r2) !------------------------------------------------------------------------ ! (T) correction diff --git a/devel/cc/CCSD_Ec_nc.irp.f b/devel/cc/CCSD_Ec_nc.irp.f new file mode 100644 index 0000000..2d57fc9 --- /dev/null +++ b/devel/cc/CCSD_Ec_nc.irp.f @@ -0,0 +1,54 @@ +subroutine CCSD_Ec_nc(nO,nV,t1,t2,Fov,EcCCSD) + +! Compute the CCSD correlatio energy in non-conanical form + + 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) + + double precision,intent(in) :: Fov(nO,nV) + +! Local variables + + integer :: i,j,a,b + +! Output variables + + double precision,intent(out) :: EcCCSD + +! Compute CCSD correlation energy + + EcCCSD = 0d0 + +! Singles contribution + + do i=1,nO + do a=1,nV + + EcCCSD = EcCCSD + Fov(i,a)*t1(i,a) + + end do + end do + +! Doubles contribution + + do i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + + EcCCSD = EcCCSD & + + 0.5d0*OOVV(i,j,a,b)*t1(i,a)*t1(j,b) & + + 0.25d0*OOVV(i,j,a,b)*t2(i,j,a,b) + + end do + end do + end do + end do + +end subroutine CCSD_Ec_nc diff --git a/devel/cc/NEED b/devel/cc/NEED index 5e2e753..9360ba0 100644 --- a/devel/cc/NEED +++ b/devel/cc/NEED @@ -1,2 +1,3 @@ hartree_fock +determinants mo_two_e_ints diff --git a/devel/cc/amplitude_guess.irp.f b/devel/cc/amplitude_guess.irp.f index a305466..06ffa72 100644 --- a/devel/cc/amplitude_guess.irp.f +++ b/devel/cc/amplitude_guess.irp.f @@ -43,9 +43,10 @@ BEGIN_PROVIDER [ double precision, t2_guess, (spin_occ_num,spin_occ_num,spin_vir ! Guess amplitudes for double excitations END_DOC - t2_guess(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + t2_guess(:,:,:,:) = OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) if (cc_guess == 1) then + t2_guess(:,:,:,:) = 0.d0 integer :: iunit integer, external :: getunitandopen character :: check diff --git a/devel/cc/denominators.irp.f b/devel/cc/denominators.irp.f index b191940..d605086 100644 --- a/devel/cc/denominators.irp.f +++ b/devel/cc/denominators.irp.f @@ -19,7 +19,8 @@ BEGIN_PROVIDER [ double precision, delta_OV, (spin_occ_num, spin_vir_num) ] do a=1,spin_vir_num do i=1,spin_occ_num - delta_OV(i,a) = eV(a) - eO(i) + delta_OV(i,a) = eO(i) - eV(a) +! delta_OV(i,a) = min(eO(i) - eV(a), -1.d-3) enddo enddo @@ -37,7 +38,7 @@ BEGIN_PROVIDER [ double precision, delta_OOVV, (spin_occ_num,spin_occ_num,spin_v 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) + delta_OOVV(i,j,a,b) = delta_OV(i,a) + delta_OV(j,b) enddo enddo enddo @@ -58,7 +59,7 @@ BEGIN_PROVIDER [ double precision, delta_OOOVVV, (spin_occ_num,spin_occ_num,spin 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) + delta_OOOVVV(i,j,k,a,b,c) = delta_OV(i,a) + delta_OV(j,b) + delta_OV(k,c) enddo enddo enddo diff --git a/devel/cc/form_cF_nc.irp.f b/devel/cc/form_cF_nc.irp.f new file mode 100644 index 0000000..31fc64f --- /dev/null +++ b/devel/cc/form_cF_nc.irp.f @@ -0,0 +1,102 @@ +subroutine form_cF_nc(nO,nV,t1,taus,Foo,Fov,Fvv,cFoo,cFov,cFvv) + +! Compute F terms in CCSD + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + + double precision,intent(in) :: t1(nO,nV) + double precision,intent(in) :: taus(nO,nO,nV,nV) + + double precision,intent(in) :: Foo(nO,nO) + double precision,intent(in) :: Fov(nO,nV) + double precision,intent(in) :: Fvv(nV,nV) + +! Local variables + + integer :: i,j,m,n + integer :: a,b,e,f + double precision,external :: Kronecker_Delta + +! Output variables + + double precision,intent(out) :: cFoo(nO,nO) + double precision,intent(out) :: cFov(nO,nV) + double precision,intent(out) :: cFvv(nV,nV) + +! Occupied-occupied block + + do m=1,nO + do i=1,nO + + cFoo(m,i) = (1d0 - Kronecker_delta(m,i))*Foo(m,i) + + do e=1,nV + cFoo(m,i) = cFoo(m,i) + 0.5d0*t1(i,e)*Fov(m,e) + end do + + do e=1,nV + do n=1,nO + cFoo(m,i) = cFoo(m,i) + t1(n,e)*OOOV(m,n,i,e) + end do + end do + + do e=1,nV + do n=1,nO + do f=1,nV + cFoo(m,i) = cFoo(m,i) + 0.5d0*taus(i,n,e,f)*OOVV(m,n,e,f) + end do + end do + end do + + end do + end do + +! Occupied-virtual block + + cFov(:,:) = Fov(:,:) + + do m=1,nO + do e=1,nV + + do n=1,nO + do f=1,nV + cFov(m,e) = cFov(m,e) + t1(n,f)*OOVV(m,n,e,f) + end do + end do + + end do + end do + +! Virtual-virtual block + + do a=1,nV + do e=1,nV + + cFvv(a,e) = (1d0 - Kronecker_delta(a,e))*Fvv(a,e) + + do m=1,nO + cFvv(a,e) = cFvv(a,e) - 0.5d0*t1(m,a)*Fov(m,e) + end do + + do m=1,nO + do f=1,nV + cFvv(a,e) = cFvv(a,e) + t1(m,f)*OVVV(m,a,f,e) + end do + end do + + do m=1,nO + do n=1,nO + do f=1,nV + cFvv(a,e) = cFvv(a,e) - 0.5d0*taus(m,n,a,f)*OOVV(m,n,e,f) + end do + end do + end do + + end do + end do + +end subroutine form_cF_nc diff --git a/devel/cc/form_cW_nc.irp.f b/devel/cc/form_cW_nc.irp.f new file mode 100644 index 0000000..cbf5d30 --- /dev/null +++ b/devel/cc/form_cW_nc.irp.f @@ -0,0 +1,105 @@ +subroutine form_cW_nc(nO,nV,t1,t2,tau,cWoooo,cWovvo,cWvvvv) + +! Compute W terms 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) + double precision,intent(in) :: tau(nO,nO,nV,nV) + +! Local variables + + integer :: i,j,m,n + integer :: a,b,e,f + double precision,external :: Kronecker_Delta + +! Output variables + + double precision,intent(out) :: cWoooo(nO,nO,nO,nO) + double precision,intent(out) :: cWovvo(nO,nV,nV,nO) + double precision,intent(out) :: cWvvvv(nV,nV,nV,nV) + +! OOOO block + + cWoooo(:,:,:,:) = OOOO(:,:,:,:) + + do m=1,nO + do n=1,nO + do i=1,nO + do j=1,nO + + do e=1,nV + cWoooo(m,n,i,j) = cWoooo(m,n,i,j) + t1(j,e)*OOOV(m,n,i,e) - t1(i,e)*OOOV(m,n,j,e) + end do + + do e=1,nV + do f=1,nV + cWoooo(m,n,i,j) = cWoooo(m,n,i,j) + 0.25d0*tau(i,j,e,f)*OOVV(m,n,e,f) + end do + end do + + end do + end do + end do + end do + +! OVVO block + + cWovvo(:,:,:,:) = OVVO(:,:,:,:) + + do m=1,nO + do b=1,nV + do e=1,nV + do j=1,nO + + do f=1,nV + cWovvo(m,b,e,j) = cWovvo(m,b,e,j) + t1(j,f)*OVVV(m,b,e,f) + end do + + do n=1,nO + cWovvo(m,b,e,j) = cWovvo(m,b,e,j) - t1(n,b)*OOVO(m,n,e,j) + end do + + do n=1,nO + do f=1,nV + cWovvo(m,b,e,j) = cWovvo(m,b,e,j) & + - ( 0.5d0*t2(j,n,f,b) + t1(j,f)*t1(n,b) )*OOVV(m,n,e,f) + end do + end do + + end do + end do + end do + end do + +! VVVV block + + cWvvvv(:,:,:,:) = VVVV(:,:,:,:) + + do a=1,nV + do b=1,nV + do e=1,nV + do f=1,nV + + do m=1,nO + cWvvvv(a,b,e,f) = cWvvvv(a,b,e,f) - t1(m,b)*VOVV(a,m,e,f) + t1(m,a)*VOVV(b,m,e,f) + end do + + do m=1,nO + do n=1,nO + cWvvvv(a,b,e,f) = cWvvvv(a,b,e,f) + 0.25d0*tau(m,n,a,b)*OOVV(m,n,e,f) + end do + end do + + end do + end do + end do + end do + + +end subroutine form_cW_nc diff --git a/devel/cc/form_r1_nc.irp.f b/devel/cc/form_r1_nc.irp.f new file mode 100644 index 0000000..82ab1f2 --- /dev/null +++ b/devel/cc/form_r1_nc.irp.f @@ -0,0 +1,77 @@ +subroutine form_r1_nc(nO,nV,t1,t2,Fov,cFoo,cFov,cFvv,r1) + +! Form residues for t1 in non-canonical 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) + + double precision,intent(in) :: Fov(nO,nV) + + double precision,intent(in) :: cFoo(nO,nO) + double precision,intent(in) :: cFov(nO,nV) + double precision,intent(in) :: cFvv(nV,nV) + +! Local variables + + integer :: i,j,m,n + integer :: a,b,e,f + +! Output variables + + double precision,intent(out) :: r1(nO,nV) + + r1(:,:) = Fov(:,:) + + do i=1,nO + do a=1,nV + + do e=1,nV + r1(i,a) = r1(i,a) + t1(i,e)*cFvv(a,e) + end do + + do m=1,nO + r1(i,a) = r1(i,a) - t1(m,a)*cFoo(m,i) + end do + + do m=1,nO + do e=1,nV + r1(i,a) = r1(i,a) + t2(i,m,a,e)*cFov(m,e) + end do + end do + + do n=1,nO + do f=1,nV + r1(i,a) = r1(i,a) - t1(n,f)*OVOV(n,a,i,f) + end do + end do + + do m=1,nO + do e=1,nV + do f=1,nV + r1(i,a) = r1(i,a) - 0.5d0*t2(i,m,e,f)*OVVV(m,a,e,f) + end do + end do + end do + + do m=1,nO + do n=1,nO + do e=1,nV + r1(i,a) = r1(i,a) - 0.5d0*t2(m,n,a,e)*OOVO(n,m,e,i) + end do + end do + end do + + end do + end do + +! Final expression for t1 residue + + r1(:,:) = delta_ov(:,:)*t1(:,:) - r1(:,:) + +end subroutine form_r1_nc diff --git a/devel/cc/form_r2_nc.irp.f b/devel/cc/form_r2_nc.irp.f new file mode 100644 index 0000000..2df576d --- /dev/null +++ b/devel/cc/form_r2_nc.irp.f @@ -0,0 +1,102 @@ +subroutine form_r2_nc(nO,nV,t1,t2,tau,cFoo,cFov,cFvv,cWoooo,cWvvvv,cWovvo,r2) + +! Form t2 residues in non-canonical CCSD + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + + double precision,intent(in) :: cFoo(nO,nO) + double precision,intent(in) :: cFov(nO,nV) + double precision,intent(in) :: cFvv(nV,nV) + + double precision,intent(in) :: cWoooo(nO,nO,nO,nO) + double precision,intent(in) :: cWvvvv(nV,nV,nV,nV) + double precision,intent(in) :: cWovvo(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,m,n + integer :: a,b,e,f + +! 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 e=1,nV + r2(i,j,a,b) = r2(i,j,a,b) + t2(i,j,a,e)*cFvv(b,e) + r2(i,j,a,b) = r2(i,j,a,b) - t2(i,j,b,e)*cFvv(a,e) + end do + + do m=1,nO + do e=1,nV + r2(i,j,a,b) = r2(i,j,a,b) - 0.5d0*t2(i,j,a,e)*t1(m,b)*cFov(m,e) + r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*t2(i,j,b,e)*t1(m,a)*cFov(m,e) + end do + end do + + do m=1,nO + r2(i,j,a,b) = r2(i,j,a,b) - t2(i,m,a,b)*cFoo(m,j) + r2(i,j,a,b) = r2(i,j,a,b) + t2(j,m,a,b)*cFoo(m,i) + end do + + do m=1,nO + do e=1,nV + r2(i,j,a,b) = r2(i,j,a,b) - 0.5d0*t2(i,m,a,b)*t1(j,e)*cFov(m,e) + r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*t2(j,m,a,b)*t1(i,e)*cFov(m,e) + end do + end do + + do m=1,nO + do n=1,nO + r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*tau(m,n,a,b)*cWoooo(m,n,i,j) + end do + end do + + do e=1,nV + do f=1,nV + r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*tau(i,j,e,f)*cWvvvv(a,b,e,f) + end do + end do + + do m=1,nO + do e=1,nV + r2(i,j,a,b) = r2(i,j,a,b) & + + t2(i,m,a,e)*cWovvo(m,b,e,j) - t1(i,e)*t1(m,a)*OVVO(m,b,e,j) & + - t2(j,m,a,e)*cWovvo(m,b,e,i) + t1(j,e)*t1(m,a)*OVVO(m,b,e,i) & + - t2(i,m,b,e)*cWovvo(m,a,e,j) + t1(i,e)*t1(m,b)*OVVO(m,a,e,j) & + + t2(j,m,b,e)*cWovvo(m,a,e,i) - t1(j,e)*t1(m,b)*OVVO(m,a,e,i) + end do + end do + + do e=1,nV + r2(i,j,a,b) = r2(i,j,a,b) + t1(i,e)*VVVO(a,b,e,j) - t1(j,e)*VVVO(a,b,e,i) + end do + + do m=1,nO + r2(i,j,a,b) = r2(i,j,a,b) - t1(m,a)*OVOO(m,b,i,j) + t1(m,b)*OVOO(m,a,i,j) + end do + + end do + end do + end do + end do + +! Final expression of the t2 residue + + r2(:,:,:,:) = delta_oovv(:,:,:,:)*t2(:,:,:,:) - r2(:,:,:,:) + +end subroutine form_r2_nc diff --git a/devel/cc/form_tau_nc.irp.f b/devel/cc/form_tau_nc.irp.f new file mode 100644 index 0000000..ea378ef --- /dev/null +++ b/devel/cc/form_tau_nc.irp.f @@ -0,0 +1,34 @@ +subroutine form_tau_nc(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) = t2(i,j,a,b) + t1(i,a)*t1(j,b) - t1(i,b)*t1(j,a) + + enddo + enddo + enddo + enddo + +end subroutine form_tau_nc diff --git a/devel/cc/form_taus_nc.irp.f b/devel/cc/form_taus_nc.irp.f new file mode 100644 index 0000000..f2a3545 --- /dev/null +++ b/devel/cc/form_taus_nc.irp.f @@ -0,0 +1,34 @@ +subroutine form_taus_nc(nO,nV,t1,t2,taus) + +! 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) :: taus(nO,nO,nV,nV) + + do b=1,nV + do a=1,nV + do j=1,nO + do i=1,nO + + taus(i,j,a,b) = t2(i,j,a,b) + 0.5d0*(t1(i,a)*t1(j,b) - t1(i,b)*t1(j,a)) + + enddo + enddo + enddo + enddo + +end subroutine form_taus_nc diff --git a/devel/cc/integral_batches.irp.f b/devel/cc/integral_batches.irp.f index 50dc42d..65bfeae 100644 --- a/devel/cc/integral_batches.irp.f +++ b/devel/cc/integral_batches.irp.f @@ -2,7 +2,11 @@ BEGIN_PROVIDER [ double precision, OOOO, (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 ) + OOOO(:,:,:,:) = dbERI( & + 1:spin_occ_num, & + 1:spin_occ_num, & + 1:spin_occ_num, & + 1:spin_occ_num ) END_PROVIDER @@ -10,7 +14,23 @@ BEGIN_PROVIDER [ double precision, OOOV, (spin_occ_num,spin_occ_num,spin_occ_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) + 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, OOVO, (spin_occ_num,spin_occ_num,spin_vir_num,spin_occ_num) ] + implicit none + BEGIN_DOC + END_DOC + OOVO(:,:,:,:) = dbERI( & + 1:spin_occ_num, & + 1:spin_occ_num, & + spin_occ_num+1:spin_mo_num, & + 1:spin_occ_num) END_PROVIDER @@ -18,7 +38,11 @@ BEGIN_PROVIDER [ double precision, OVOO, (spin_occ_num,spin_vir_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 ) + OVOO(:,:,:,:) = dbERI( & + 1:spin_occ_num, & + spin_occ_num+1:spin_mo_num, & + 1:spin_occ_num, & + 1:spin_occ_num) END_PROVIDER @@ -26,7 +50,11 @@ BEGIN_PROVIDER [ double precision, VOOO, (spin_vir_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 ) + VOOO(:,:,:,:) = dbERI( & + spin_occ_num+1:spin_mo_num, & + 1:spin_occ_num, & + 1:spin_occ_num, & + 1:spin_occ_num) END_PROVIDER @@ -34,7 +62,23 @@ BEGIN_PROVIDER [ double precision, OOVV, (spin_occ_num,spin_occ_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) + 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, OVOV, (spin_occ_num,spin_vir_num,spin_occ_num,spin_vir_num) ] + implicit none + BEGIN_DOC + END_DOC + OVOV(:,:,:,:) = dbERI( & + 1:spin_occ_num, & + spin_occ_num+1:spin_mo_num, & + 1:spin_occ_num, & + spin_occ_num+1:spin_mo_num) END_PROVIDER @@ -42,7 +86,35 @@ BEGIN_PROVIDER [ double precision, OVVO, (spin_occ_num,spin_vir_num,spin_vir_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 ) + 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, VOVO, (spin_vir_num,spin_occ_num,spin_vir_num,spin_occ_num) ] + implicit none + BEGIN_DOC + END_DOC + VOVO(:,:,:,:) = dbERI( & + spin_occ_num+1:spin_mo_num, & + 1:spin_occ_num, & + spin_occ_num+1:spin_mo_num, & + 1:spin_occ_num) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, VVOO, (spin_vir_num,spin_vir_num,spin_occ_num,spin_occ_num) ] + implicit none + BEGIN_DOC + END_DOC + VVOO(:,:,:,:) = dbERI( & + spin_occ_num+1:spin_mo_num, & + spin_occ_num+1:spin_mo_num, & + 1:spin_occ_num, & + 1:spin_occ_num) END_PROVIDER @@ -50,15 +122,34 @@ BEGIN_PROVIDER [ double precision, OVVV, (spin_occ_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) + 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) + 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, VVOV, (spin_vir_num,spin_vir_num,spin_occ_num,spin_vir_num) ] + implicit none + BEGIN_DOC + END_DOC + VVOV(:,:,:,:) = dbERI( & + spin_occ_num+1:spin_mo_num, & + spin_occ_num+1:spin_mo_num, & + 1:spin_occ_num, & + spin_occ_num+1:spin_mo_num) END_PROVIDER @@ -66,14 +157,21 @@ BEGIN_PROVIDER [ double precision, VVVO, (spin_vir_num,spin_vir_num,spin_vir_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 ) + 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) + 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/spatial_to_spin_MO.irp.f b/devel/cc/spatial_to_spin_MO.irp.f new file mode 100644 index 0000000..6342f59 --- /dev/null +++ b/devel/cc/spatial_to_spin_MO.irp.f @@ -0,0 +1,120 @@ +BEGIN_PROVIDER [ double precision, spin_fock_matrix_mo, (spin_mo_num, spin_mo_num) ] + implicit none + BEGIN_DOC + ! Fock matrix in the spin-orbital basis + END_DOC + integer :: p,q + double precision :: F(mo_num,mo_num) + + spin_fock_matrix_mo = 0.d0 + call get_fock_matrix_alpha(hf_bitmask,F) + do q=1,spin_mo_num,2 + do p=1,spin_mo_num,2 + spin_fock_matrix_mo(p,q) = F((p+1)/2, (q+1)/2) + enddo + enddo + + call get_fock_matrix_beta (hf_bitmask,F) + do q=2,spin_mo_num,2 + do p=2,spin_mo_num,2 + spin_fock_matrix_mo(p,q) = F((p+1)/2, (q+1)/2) + enddo + enddo + +END_PROVIDER + +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) = spin_fock_matrix_mo(p,p) + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, spin_fock_matrix_mo_oo, (spin_occ_num, spin_occ_num) ] + implicit none + BEGIN_DOC + ! Occupied-Occupied block of the Fock matrix in the spin-orbital basis + END_DOC + integer :: p,q + + do q=1,spin_occ_num + do p=1,spin_occ_num + spin_fock_matrix_mo_oo(p,q) = spin_fock_matrix_mo(p,q) + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, spin_fock_matrix_mo_vv, (spin_vir_num, spin_vir_num) ] + implicit none + BEGIN_DOC + ! Virtual-Virtual block of the Fock matrix in the spin-orbital basis + END_DOC + integer :: p,q + + do q=1,spin_vir_num + do p=1,spin_vir_num + spin_fock_matrix_mo_vv(p,q) = spin_fock_matrix_mo(p+spin_occ_num, q+spin_occ_num) + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, spin_fock_matrix_mo_ov, (spin_occ_num, spin_vir_num) ] + implicit none + BEGIN_DOC + ! Occupied-Virtual block of the Fock matrix in the spin-orbital basis + END_DOC + integer :: p,q + + do q=1,spin_vir_num + do p=1,spin_occ_num + spin_fock_matrix_mo_ov(p,q) = spin_fock_matrix_mo(p, q+spin_occ_num) + enddo + enddo + +END_PROVIDER + + + +subroutine get_fock_matrix_alpha(det,F) + implicit none + BEGIN_DOC +! Returns the alpha Fock matrix in MO basis associated with the determinant given as input + END_DOC + integer(bit_kind), intent(in) :: det(N_int,2) + double precision, intent(out) :: F(mo_num,mo_num) + + integer :: i,j,k + + F(:,:) = fock_operator_closed_shell_ref_bitmask(:,:) + +end + + + +subroutine get_fock_matrix_beta(det,F) + implicit none + BEGIN_DOC +! Returns the beta Fock matrix in MO basis associated with the determinant given as input + END_DOC + integer(bit_kind), intent(in) :: det(N_int,2) + double precision, intent(out) :: F(mo_num,mo_num) + + integer :: i,j,k + + F(:,:) = fock_operator_closed_shell_ref_bitmask(:,:) + +end + + + diff --git a/devel/cc/spatial_to_spin_MO_energy.irp.f b/devel/cc/spatial_to_spin_MO_energy.irp.f deleted file mode 100644 index 0ce0576..0000000 --- a/devel/cc/spatial_to_spin_MO_energy.irp.f +++ /dev/null @@ -1,13 +0,0 @@ -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 - diff --git a/stable/amplitudes/extract_amplitudes.irp.f b/stable/amplitudes/extract_amplitudes.irp.f index fe0e99a..cb08637 100644 --- a/stable/amplitudes/extract_amplitudes.irp.f +++ b/stable/amplitudes/extract_amplitudes.irp.f @@ -43,16 +43,20 @@ subroutine run t2_ab = 0.d0 t2_bb = 0.d0 - norm = 1.d0 / psi_coef_sorted(1,istate) + integer :: i_ref + + i_ref = maxloc( dabs(psi_coef(:,istate)), dim=1 ) + + norm = 1.d0 / psi_coef(i_ref,istate) do k=1,N_det - call get_excitation_degree(hf_bitmask,psi_det(1,1,k),degree,N_int) + call get_excitation_degree(psi_det(1,1,i_ref),psi_det(1,1,k),degree,N_int) select case (degree) case (1) - call get_excitation(hf_bitmask,psi_det(1,1,k),exc,degree,phase,N_int) + call get_excitation(psi_det(1,1,i_ref),psi_det(1,1,k),exc,degree,phase,N_int) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) if (h1 > elec_alpha_num) then @@ -72,7 +76,7 @@ subroutine run case (2) - call get_excitation(hf_bitmask,psi_det(1,1,k),exc,degree,phase,N_int) + call get_excitation(psi_det(1,1,i_ref),psi_det(1,1,k),exc,degree,phase,N_int) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) if (h1 > elec_alpha_num) then @@ -200,15 +204,19 @@ subroutine run enddo - double precision :: e_cor + double precision :: E0, e_cor double precision, external :: get_two_e_integral e_cor = 0.d0 do b=elec_alpha_num+1, mo_num - do a=elec_alpha_num+1, mo_num - do j=1,elec_alpha_num + do j=1,elec_alpha_num +! TODO : singles contributions +! do k=1,mo_num +! e_cor = e_cor + t1_a(j,b) * get_two_e_integral(k,j,k,b,mo_integrals_map) +! e_cor = e_cor + t1_b(j,b) * get_two_e_integral(k,j,k,b,mo_integrals_map) +! enddo + do a=elec_alpha_num+1, mo_num do i=1,elec_alpha_num - ! The t1 terms are zero because HF e_cor = e_cor + 0.25d0*( & t2_aa(i,j,a,b) + t2_bb(i,j,a,b) + & t1_a(i,a) * t1_a(j,b) + & @@ -227,9 +235,12 @@ subroutine run enddo enddo e_cor = e_cor - print '(A,F15.10)', 'HF energy: ', hf_energy + + double precision, external :: diag_h_mat_elem + E0 = diag_h_mat_elem(psi_det(1,1,i_ref),N_int) + nuclear_repulsion + print '(A,F15.10)', 'E0 : ', E0 print '(A,F15.10)', 'corr energy: ', e_cor - print '(A,F15.10)', 'total energy: ', hf_energy + e_cor + print '(A,F15.10)', 'total energy: ', E0 + e_cor deallocate(t1_a,t1_b,t2_aa,t2_ab,t2_bb)