From 5ab11e50bb85631deaa9a82ad56ea0b630a0c15f Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 11 Sep 2019 14:47:18 +0200 Subject: [PATCH] non-canonical CCSD --- input/basis | 60 +++++++++++-------- input/molecule | 4 +- input/weight | 60 +++++++++++-------- src/QuAcK/CCSD_Ec_nc.f90 | 56 ++++++++++++++++++ src/QuAcK/form_cF_nc.f90 | 107 ++++++++++++++++++++++++++++++++++ src/QuAcK/form_cW_nc.f90 | 114 +++++++++++++++++++++++++++++++++++++ src/QuAcK/form_r1_nc.f90 | 82 ++++++++++++++++++++++++++ src/QuAcK/form_r2_nc.f90 | 109 +++++++++++++++++++++++++++++++++++ src/QuAcK/form_tau_nc.f90 | 34 +++++++++++ src/QuAcK/form_taus_nc.f90 | 34 +++++++++++ 10 files changed, 608 insertions(+), 52 deletions(-) create mode 100644 src/QuAcK/CCSD_Ec_nc.f90 create mode 100644 src/QuAcK/form_cF_nc.f90 create mode 100644 src/QuAcK/form_cW_nc.f90 create mode 100644 src/QuAcK/form_r1_nc.f90 create mode 100644 src/QuAcK/form_r2_nc.f90 create mode 100644 src/QuAcK/form_tau_nc.f90 create mode 100644 src/QuAcK/form_taus_nc.f90 diff --git a/input/basis b/input/basis index eb58543..35be48a 100644 --- a/input/basis +++ b/input/basis @@ -1,29 +1,39 @@ -1 6 -S 8 1.00 - 17880.0000000 0.0007380 - 2683.0000000 0.0056770 - 611.5000000 0.0288830 - 173.5000000 0.1085400 - 56.6400000 0.2909070 - 20.4200000 0.4483240 - 7.8100000 0.2580260 - 1.6530000 0.0150630 -S 8 1.00 - 17880.0000000 -0.0001720 - 2683.0000000 -0.0013570 - 611.5000000 -0.0067370 - 173.5000000 -0.0276630 - 56.6400000 -0.0762080 - 20.4200000 -0.1752270 - 7.8100000 -0.1070380 - 1.6530000 0.5670500 +1 10 +S 9 1.00 + 6863.0000000 0.0002360 + 1030.0000000 0.0018260 + 234.7000000 0.0094520 + 66.5600000 0.0379570 + 21.6900000 0.1199650 + 7.7340000 0.2821620 + 2.9160000 0.4274040 + 1.1300000 0.2662780 + 0.1101000 -0.0072750 +S 9 1.00 + 6863.0000000 -0.0000430 + 1030.0000000 -0.0003330 + 234.7000000 -0.0017360 + 66.5600000 -0.0070120 + 21.6900000 -0.0231260 + 7.7340000 -0.0581380 + 2.9160000 -0.1145560 + 1.1300000 -0.1359080 + 0.1101000 0.5774410 S 1 1.00 - 0.4869000 1.0000000 + 0.2577000 1.0000000 +S 1 1.00 + 0.0440900 1.0000000 P 3 1.00 - 28.3900000 0.0460870 - 6.2700000 0.2401810 - 1.6950000 0.5087440 + 7.4360000 0.0107360 + 1.5770000 0.0628540 + 0.4352000 0.2481800 P 1 1.00 - 0.4317000 1.0000000 + 0.1438000 1.0000000 +P 1 1.00 + 0.0499400 1.0000000 D 1 1.00 - 2.2020000 1.0000000 + 0.3480000 1.0000000 +D 1 1.00 + 0.1803000 1.0000000 +F 1 1.00 + 0.3250000 1.0000000 diff --git a/input/molecule b/input/molecule index edeba31..6a6f6d1 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 5 5 0 0 + 1 2 2 0 0 # Znuc x y z - Ne 0.0 0.0 0.0 + Be 0.0 0.0 0.0 diff --git a/input/weight b/input/weight index eb58543..35be48a 100644 --- a/input/weight +++ b/input/weight @@ -1,29 +1,39 @@ -1 6 -S 8 1.00 - 17880.0000000 0.0007380 - 2683.0000000 0.0056770 - 611.5000000 0.0288830 - 173.5000000 0.1085400 - 56.6400000 0.2909070 - 20.4200000 0.4483240 - 7.8100000 0.2580260 - 1.6530000 0.0150630 -S 8 1.00 - 17880.0000000 -0.0001720 - 2683.0000000 -0.0013570 - 611.5000000 -0.0067370 - 173.5000000 -0.0276630 - 56.6400000 -0.0762080 - 20.4200000 -0.1752270 - 7.8100000 -0.1070380 - 1.6530000 0.5670500 +1 10 +S 9 1.00 + 6863.0000000 0.0002360 + 1030.0000000 0.0018260 + 234.7000000 0.0094520 + 66.5600000 0.0379570 + 21.6900000 0.1199650 + 7.7340000 0.2821620 + 2.9160000 0.4274040 + 1.1300000 0.2662780 + 0.1101000 -0.0072750 +S 9 1.00 + 6863.0000000 -0.0000430 + 1030.0000000 -0.0003330 + 234.7000000 -0.0017360 + 66.5600000 -0.0070120 + 21.6900000 -0.0231260 + 7.7340000 -0.0581380 + 2.9160000 -0.1145560 + 1.1300000 -0.1359080 + 0.1101000 0.5774410 S 1 1.00 - 0.4869000 1.0000000 + 0.2577000 1.0000000 +S 1 1.00 + 0.0440900 1.0000000 P 3 1.00 - 28.3900000 0.0460870 - 6.2700000 0.2401810 - 1.6950000 0.5087440 + 7.4360000 0.0107360 + 1.5770000 0.0628540 + 0.4352000 0.2481800 P 1 1.00 - 0.4317000 1.0000000 + 0.1438000 1.0000000 +P 1 1.00 + 0.0499400 1.0000000 D 1 1.00 - 2.2020000 1.0000000 + 0.3480000 1.0000000 +D 1 1.00 + 0.1803000 1.0000000 +F 1 1.00 + 0.3250000 1.0000000 diff --git a/src/QuAcK/CCSD_Ec_nc.f90 b/src/QuAcK/CCSD_Ec_nc.f90 new file mode 100644 index 0000000..f55f48f --- /dev/null +++ b/src/QuAcK/CCSD_Ec_nc.f90 @@ -0,0 +1,56 @@ +subroutine CCSD_Ec_nc(nO,nV,t1,t2,Fov,OOVV,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) + + double precision,intent(in) :: OOVV(nO,nO,nV,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,nO + + 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,nO + do b=1,nO + + 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/src/QuAcK/form_cF_nc.f90 b/src/QuAcK/form_cF_nc.f90 new file mode 100644 index 0000000..6ebe8d1 --- /dev/null +++ b/src/QuAcK/form_cF_nc.f90 @@ -0,0 +1,107 @@ +subroutine form_cF_nc(nO,nV,t1,taus,Foo,Fov,Fvv,OOOV,OOVV,OVVV,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) + + double precision,intent(in) :: OOOV(nO,nO,nO,nV) + double precision,intent(in) :: OOVV(nO,nO,nV,nV) + double precision,intent(in) :: OVVV(nO,nV,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/src/QuAcK/form_cW_nc.f90 b/src/QuAcK/form_cW_nc.f90 new file mode 100644 index 0000000..9c0c2a2 --- /dev/null +++ b/src/QuAcK/form_cW_nc.f90 @@ -0,0 +1,114 @@ +subroutine form_cW_nc(nO,nV,t1,t2,tau,OOOO,OOOV,OOVO,OOVV,OVVO,OVVV,VOVV,VVVV,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) + + double precision,intent(in) :: OOOO(nO,nO,nO,nO) + double precision,intent(in) :: OOOV(nO,nO,nO,nV) + double precision,intent(in) :: OOVO(nO,nO,nV,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) + +! 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/src/QuAcK/form_r1_nc.f90 b/src/QuAcK/form_r1_nc.f90 new file mode 100644 index 0000000..5b1337d --- /dev/null +++ b/src/QuAcK/form_r1_nc.f90 @@ -0,0 +1,82 @@ +subroutine form_r1_nc(nO,nV,t1,t2,delta_ov,Fov,cFoo,cFov,cFvv,OOVO,OVOV,OVVV,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) :: delta_ov(nO,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) + + double precision,intent(in) :: OOVO(nO,nO,nV,nO) + double precision,intent(in) :: OVOV(nO,nV,nO,nV) + double precision,intent(in) :: OVVV(nO,nV,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(e,a) + 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/src/QuAcK/form_r2_nc.f90 b/src/QuAcK/form_r2_nc.f90 new file mode 100644 index 0000000..9ca8b4d --- /dev/null +++ b/src/QuAcK/form_r2_nc.f90 @@ -0,0 +1,109 @@ +subroutine form_r2_nc(nO,nV,t1,t2,tau,delta_oovv,cFoo,cFov,cFvv,cWoooo,cWvvvv,cWovvo,OVOO,OOVV,OVVO,VVVO,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) :: delta_oovv(nO,nO,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) + + 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) :: VVVO(nV,nV,nV,nO) + +! 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/src/QuAcK/form_tau_nc.f90 b/src/QuAcK/form_tau_nc.f90 new file mode 100644 index 0000000..c7f4aaa --- /dev/null +++ b/src/QuAcK/form_tau_nc.f90 @@ -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 i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + + 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/src/QuAcK/form_taus_nc.f90 b/src/QuAcK/form_taus_nc.f90 new file mode 100644 index 0000000..b834c9f --- /dev/null +++ b/src/QuAcK/form_taus_nc.f90 @@ -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 i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + + 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