scemama/qp_plugins_scemama
1
0
Fork 0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-11-07 14:43:41 +01:00
Code Issues Releases Wiki Activity

fork CC from mveril and pfloos

Browse Source
This commit is contained in:
Anthony Scemama 2019-09-09 16:51:15 +02:00
parent ab8ff0c2ba
commit c69525722a
25 changed files with 1302 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
devel/cc/CC.irp.f Normal file
View File

@ -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

211
devel/cc/CCSD.irp.f Normal file
View File

@ -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

35
devel/cc/CCSDT.irp.f Normal file
View File

@ -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

22
devel/cc/EZFIO.cfg Normal file
View File

@ -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

22
devel/cc/Kronecker_delta.irp.f Normal file
View File

@ -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

71
devel/cc/MP2.irp.f Normal file
View File

@ -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

2
devel/cc/NEED Normal file
View File

@ -0,0 +1,2 @@
hartree_fock
mo_two_e_ints

4
devel/cc/README.rst Normal file
View File

@ -0,0 +1,4 @@
==
CC
==

19
devel/cc/amplitude_guess.irp.f Normal file
View File

@ -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

46
devel/cc/antisymmetrize_ERI.irp.f Normal file
View File

@ -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

70
devel/cc/denominators.irp.f Normal file
View File

@ -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

20
devel/cc/dimensions.irp.f Normal file
View File

@ -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

45
devel/cc/form_T.irp.f Normal file
View File

@ -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

97
devel/cc/form_abh.irp.f Normal file
View File

@ -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

50
devel/cc/form_g.irp.f Normal file
View File

@ -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

75
devel/cc/form_h.irp.f Normal file
View File

@ -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

72
devel/cc/form_r1.irp.f Normal file
View File

@ -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

134
devel/cc/form_r2.irp.f Normal file
View File

@ -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

34
devel/cc/form_tau.irp.f Normal file
View File

@ -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

46
devel/cc/form_ub.irp.f Normal file
View File

@ -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

64
devel/cc/form_ubb.irp.f Normal file
View File

@ -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

79
devel/cc/integral_batches.irp.f Normal file
View File

@ -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

8
devel/cc/save_energy.irp.f Normal file
View File

@ -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

52
devel/cc/spatial_to_spin_ERI.irp.f Normal file
View File

@ -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

13
devel/cc/spatial_to_spin_MO_energy.irp.f Normal file
View File

@ -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