mirror of
https://gitlab.com/scemama/qp_plugins_scemama.git
synced 2024-09-01 08:43:46 +02:00
Compare commits
3 Commits
dd65b2677b
...
aee13ced42
Author | SHA1 | Date | |
---|---|---|---|
aee13ced42 | |||
c69525722a | |||
ab8ff0c2ba |
11
devel/cc/CC.irp.f
Normal file
11
devel/cc/CC.irp.f
Normal 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
211
devel/cc/CCSD.irp.f
Normal 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
35
devel/cc/CCSDT.irp.f
Normal 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
|
28
devel/cc/EZFIO.cfg
Normal file
28
devel/cc/EZFIO.cfg
Normal file
@ -0,0 +1,28 @@
|
|||||||
|
[energy]
|
||||||
|
type: Threshold
|
||||||
|
doc: Energy CC
|
||||||
|
interface: ezfio
|
||||||
|
|
||||||
|
[cc_guess]
|
||||||
|
type: integer
|
||||||
|
doc: 0: MP2 amplitudes, 1: read from files
|
||||||
|
interface: ezfio,provider,ocaml
|
||||||
|
default: 0
|
||||||
|
|
||||||
|
[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
22
devel/cc/Kronecker_delta.irp.f
Normal 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
71
devel/cc/MP2.irp.f
Normal 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
2
devel/cc/NEED
Normal file
@ -0,0 +1,2 @@
|
|||||||
|
hartree_fock
|
||||||
|
mo_two_e_ints
|
4
devel/cc/README.rst
Normal file
4
devel/cc/README.rst
Normal file
@ -0,0 +1,4 @@
|
|||||||
|
==
|
||||||
|
CC
|
||||||
|
==
|
||||||
|
|
84
devel/cc/amplitude_guess.irp.f
Normal file
84
devel/cc/amplitude_guess.irp.f
Normal file
@ -0,0 +1,84 @@
|
|||||||
|
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
|
||||||
|
if (cc_guess == 1) then
|
||||||
|
integer :: iunit
|
||||||
|
integer, external :: getunitandopen
|
||||||
|
|
||||||
|
character :: check
|
||||||
|
integer :: i, a
|
||||||
|
double precision :: amplitude
|
||||||
|
|
||||||
|
iunit = getunitandopen('t1','r')
|
||||||
|
read(iunit,*)
|
||||||
|
do
|
||||||
|
read(iunit,*,err=10) i, a, amplitude
|
||||||
|
i = 2*i-1
|
||||||
|
a = 2*a-1 - spin_occ_num
|
||||||
|
t1_guess(i,a) = amplitude
|
||||||
|
enddo
|
||||||
|
10 continue
|
||||||
|
do
|
||||||
|
read(iunit,*,end=20) i, a, amplitude
|
||||||
|
i = 2*i
|
||||||
|
a = 2*a - spin_occ_num
|
||||||
|
t1_guess(i,a) = amplitude
|
||||||
|
enddo
|
||||||
|
20 continue
|
||||||
|
close(iunit)
|
||||||
|
endif
|
||||||
|
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(:,:,:,:)
|
||||||
|
|
||||||
|
if (cc_guess == 1) then
|
||||||
|
integer :: iunit
|
||||||
|
integer, external :: getunitandopen
|
||||||
|
character :: check
|
||||||
|
integer :: i, j, a, b
|
||||||
|
double precision :: amplitude
|
||||||
|
|
||||||
|
iunit = getunitandopen('t2','r')
|
||||||
|
read(iunit,*)
|
||||||
|
do
|
||||||
|
read(iunit,*,err=10) i, j, a, b, amplitude
|
||||||
|
i = 2*i-1
|
||||||
|
j = 2*j-1
|
||||||
|
a = 2*a-1 - spin_occ_num
|
||||||
|
b = 2*b-1 - spin_occ_num
|
||||||
|
t2_guess(i,j,a,b) = amplitude
|
||||||
|
enddo
|
||||||
|
10 continue
|
||||||
|
do
|
||||||
|
read(iunit,*,err=20) i, j, a, b, amplitude
|
||||||
|
i = 2*i
|
||||||
|
j = 2*j
|
||||||
|
a = 2*a - spin_occ_num
|
||||||
|
b = 2*b - spin_occ_num
|
||||||
|
t2_guess(i,j,a,b) = amplitude
|
||||||
|
enddo
|
||||||
|
20 continue
|
||||||
|
do
|
||||||
|
read(iunit,*,end=30) i, j, a, b, amplitude
|
||||||
|
i = 2*i-1
|
||||||
|
j = 2*j
|
||||||
|
a = 2*a-1 - spin_occ_num
|
||||||
|
b = 2*b - spin_occ_num
|
||||||
|
t2_guess(i,j,a,b) = amplitude
|
||||||
|
enddo
|
||||||
|
30 continue
|
||||||
|
close(iunit)
|
||||||
|
endif
|
||||||
|
|
||||||
|
END_PROVIDER
|
46
devel/cc/antisymmetrize_ERI.irp.f
Normal file
46
devel/cc/antisymmetrize_ERI.irp.f
Normal 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
70
devel/cc/denominators.irp.f
Normal 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
20
devel/cc/dimensions.irp.f
Normal 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
45
devel/cc/form_T.irp.f
Normal 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
97
devel/cc/form_abh.irp.f
Normal 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
50
devel/cc/form_g.irp.f
Normal 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
75
devel/cc/form_h.irp.f
Normal 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
72
devel/cc/form_r1.irp.f
Normal 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
134
devel/cc/form_r2.irp.f
Normal 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
34
devel/cc/form_tau.irp.f
Normal 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
46
devel/cc/form_ub.irp.f
Normal 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
64
devel/cc/form_ubb.irp.f
Normal 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
79
devel/cc/integral_batches.irp.f
Normal 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
8
devel/cc/save_energy.irp.f
Normal 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
52
devel/cc/spatial_to_spin_ERI.irp.f
Normal 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
13
devel/cc/spatial_to_spin_MO_energy.irp.f
Normal 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
|
||||||
|
|
212
stable/utilities/extract_amplitudes.irp.f
Normal file
212
stable/utilities/extract_amplitudes.irp.f
Normal file
@ -0,0 +1,212 @@
|
|||||||
|
program extract_amplitudes
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Print the T1 and T2 amplitudes into a file.
|
||||||
|
END_DOC
|
||||||
|
|
||||||
|
read_wf = .True.
|
||||||
|
TOUCH read_wf
|
||||||
|
call run
|
||||||
|
end
|
||||||
|
|
||||||
|
subroutine run
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Compute t1 and T2 amplitudes
|
||||||
|
END_DOC
|
||||||
|
double precision, allocatable :: t1_a(:,:), t1_b(:,:)
|
||||||
|
double precision, allocatable :: t2_aa(:,:,:,:)
|
||||||
|
double precision, allocatable :: t2_ab(:,:,:,:)
|
||||||
|
double precision, allocatable :: t2_bb(:,:,:,:)
|
||||||
|
double precision :: phase, norm
|
||||||
|
integer :: exc(0:2,2,2), degree
|
||||||
|
integer :: h1,h2,p1,p2,s1,s2, istate
|
||||||
|
integer :: i,j,k,l,a,b,c,d
|
||||||
|
|
||||||
|
istate = 1
|
||||||
|
|
||||||
|
allocate ( &
|
||||||
|
t1_a(elec_alpha_num,mo_num), &
|
||||||
|
t1_b(elec_alpha_num,mo_num), &
|
||||||
|
t2_aa(elec_alpha_num,elec_alpha_num,mo_num,mo_num), &
|
||||||
|
t2_ab(elec_alpha_num,elec_alpha_num,mo_num,mo_num), &
|
||||||
|
t2_bb(elec_alpha_num,elec_alpha_num,mo_num,mo_num) )
|
||||||
|
|
||||||
|
t1_a = 0.d0
|
||||||
|
t1_b = 0.d0
|
||||||
|
t2_aa = 0.d0
|
||||||
|
t2_ab = 0.d0
|
||||||
|
t2_bb = 0.d0
|
||||||
|
|
||||||
|
norm = 1.d0 / psi_coef_sorted(1,istate)
|
||||||
|
|
||||||
|
do k=1,N_det
|
||||||
|
|
||||||
|
call get_excitation_degree(hf_bitmask,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 decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||||
|
|
||||||
|
if (h1 > elec_alpha_num) then
|
||||||
|
print *, irp_here, "h1 > elec_alpha_num"
|
||||||
|
endif
|
||||||
|
if (p1 <= elec_alpha_num) then
|
||||||
|
print *, irp_here, "p1 <= elec_alpha_num"
|
||||||
|
endif
|
||||||
|
|
||||||
|
if (s1 == 1) then
|
||||||
|
t1_a(h1,p1) += phase * psi_coef(k,istate) * norm
|
||||||
|
else if (s1 == 2) then
|
||||||
|
t1_b(h1,p1) += phase * psi_coef(k,istate) * norm
|
||||||
|
else
|
||||||
|
print *, irp_here, ': Bug!', s1, s2
|
||||||
|
endif
|
||||||
|
|
||||||
|
|
||||||
|
case (2)
|
||||||
|
call get_excitation(hf_bitmask,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
|
||||||
|
print *, irp_here, "h1 > elec_alpha_num"
|
||||||
|
endif
|
||||||
|
if (p1 <= elec_alpha_num) then
|
||||||
|
print *, irp_here, "p1 <= elec_alpha_num"
|
||||||
|
endif
|
||||||
|
if (h2 > elec_alpha_num) then
|
||||||
|
print *, irp_here, "h2 > elec_alpha_num"
|
||||||
|
endif
|
||||||
|
if (p2 <= elec_alpha_num) then
|
||||||
|
print *, irp_here, "p2 <= elec_alpha_num"
|
||||||
|
endif
|
||||||
|
|
||||||
|
if ( (s1 == 1).and.(s2 == 1) ) then
|
||||||
|
t2_aa(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm
|
||||||
|
else if ( (s1 == 1).and.(s2 == 2) ) then
|
||||||
|
t2_ab(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm
|
||||||
|
else if ( (s1 == 2).and.(s2 == 1) ) then
|
||||||
|
t2_ab(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm
|
||||||
|
else if ( (s1 == 2).and.(s2 == 2) ) then
|
||||||
|
t2_bb(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm
|
||||||
|
else
|
||||||
|
print *, irp_here, ': Bug!'
|
||||||
|
endif
|
||||||
|
|
||||||
|
end select
|
||||||
|
|
||||||
|
enddo
|
||||||
|
|
||||||
|
do b=elec_alpha_num+1, mo_num
|
||||||
|
do a=elec_alpha_num+1, mo_num
|
||||||
|
do j=1,elec_alpha_num
|
||||||
|
do i=1,elec_alpha_num
|
||||||
|
t2_aa(i,j,a,b) -= t1_a(i,a)*t1_a(j,b)
|
||||||
|
t2_bb(i,j,a,b) -= t1_b(i,a)*t1_b(j,b)
|
||||||
|
t2_ab(i,j,a,b) -= t1_a(i,a)*t1_b(j,b)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
integer :: iunit, getunitandopen
|
||||||
|
|
||||||
|
iunit = getunitandopen('t1','w')
|
||||||
|
|
||||||
|
write(iunit,*) '# T1_a'
|
||||||
|
do i=1,elec_alpha_num
|
||||||
|
do a=elec_alpha_num+1, mo_num
|
||||||
|
if (t1_a(i,a) /= 0.d0) then
|
||||||
|
write(iunit,'(I4,X,I4,X,E20.10)') i, a, t1_a(i,a)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
write(iunit,*) '# T1_b'
|
||||||
|
do i=1,elec_alpha_num
|
||||||
|
do a=elec_alpha_num+1, mo_num
|
||||||
|
if (t1_b(i,a) /= 0.d0) then
|
||||||
|
write(iunit,'(I4,X,I4,X,E20.10)') i, a, t1_b(i,a)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
iunit = getunitandopen('t2','w')
|
||||||
|
|
||||||
|
write(iunit,*) '# T2_aa'
|
||||||
|
do i=1,elec_alpha_num
|
||||||
|
do j=1,elec_alpha_num
|
||||||
|
do a=elec_alpha_num+1, mo_num
|
||||||
|
do b=elec_alpha_num+1, mo_num
|
||||||
|
if (t2_aa(i,j,a,b) /= 0.d0) then
|
||||||
|
write(iunit,'(4(I4,X),E20.10)') i, j, a, b, t2_aa(i,j,a,b)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
write(iunit,*) '# T2_bb'
|
||||||
|
do i=1,elec_alpha_num
|
||||||
|
do j=1,elec_alpha_num
|
||||||
|
do a=elec_alpha_num+1, mo_num
|
||||||
|
do b=elec_alpha_num+1, mo_num
|
||||||
|
if (t2_bb(i,j,a,b) /= 0.d0) then
|
||||||
|
write(iunit,'(4(I4,X),E20.10)') i, j, a, b, t2_bb(i,j,a,b)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
write(iunit,*) '# T2_ab'
|
||||||
|
do i=1,elec_alpha_num
|
||||||
|
do j=1,elec_alpha_num
|
||||||
|
do a=elec_alpha_num+1, mo_num
|
||||||
|
do b=elec_alpha_num+1, mo_num
|
||||||
|
if (t2_ab(i,j,a,b) /= 0.d0) then
|
||||||
|
write(iunit,'(4(I4,X),E20.10)') i, j, a, b, t2_ab(i,j,a,b)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
double precision :: 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 i=1,elec_alpha_num
|
||||||
|
! The t1 terms are zero because HF
|
||||||
|
e_cor = e_cor + ( &
|
||||||
|
t2_aa(i,j,a,b) + t2_bb(i,j,a,b) + &
|
||||||
|
t1_a(i,a) * t1_a(j,b) + &
|
||||||
|
t1_b(i,a) * t1_b(j,b) ) * ( &
|
||||||
|
get_two_e_integral(i,j,a,b,mo_integrals_map) - &
|
||||||
|
get_two_e_integral(i,j,b,a,mo_integrals_map) )
|
||||||
|
e_cor = e_cor + ( &
|
||||||
|
t2_ab(i,j,a,b) + &
|
||||||
|
t1_a(i,a) * t1_b(j,b) ) * &
|
||||||
|
get_two_e_integral(i,j,a,b,mo_integrals_map)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
e_cor = e_cor
|
||||||
|
print '(A,F15.10)', 'HF energy: ', hf_energy
|
||||||
|
print '(A,F15.10)', 'corr energy: ', e_cor
|
||||||
|
print '(A,F15.10)', 'total energy: ', hf_energy + e_cor
|
||||||
|
|
||||||
|
deallocate(t1_a,t1_b,t2_aa,t2_ab,t2_bb)
|
||||||
|
|
||||||
|
end
|
Loading…
Reference in New Issue
Block a user