mirror of
https://github.com/pfloos/quack
synced 2024-11-06 22:24:03 +01:00
BCCD start
This commit is contained in:
parent
e0a28f71bd
commit
e793ceaeca
@ -1,11 +1,11 @@
|
|||||||
# RHF UHF KS MOM
|
# RHF UHF KS MOM
|
||||||
F T F F
|
T F F F
|
||||||
# MP2* MP3 MP2-F12
|
# MP2* MP3 MP2-F12
|
||||||
F F F
|
F F F
|
||||||
# CCD DCD CCSD CCSD(T)
|
# CCD DCD CCSD CCSD(T)
|
||||||
F F F F
|
F F F F
|
||||||
# drCCD rCCD lCCD pCCD
|
# drCCD rCCD lCCD pCCD
|
||||||
F F F F
|
F F F T
|
||||||
# CIS* CIS(D) CID CISD
|
# CIS* CIS(D) CID CISD
|
||||||
F F F F
|
F F F F
|
||||||
# RPA* RPAx* ppRPA
|
# RPA* RPAx* ppRPA
|
||||||
@ -13,7 +13,7 @@
|
|||||||
# G0F2 evGF2 G0F3 evGF3
|
# G0F2 evGF2 G0F3 evGF3
|
||||||
F F F F
|
F F F F
|
||||||
# G0W0* evGW* qsGW*
|
# G0W0* evGW* qsGW*
|
||||||
F F T
|
F F F
|
||||||
# G0T0 evGT qsGT
|
# G0T0 evGT qsGT
|
||||||
F F F
|
F F F
|
||||||
# MCMP2
|
# MCMP2
|
||||||
|
|||||||
@ -1,4 +1,4 @@
|
|||||||
2
|
2
|
||||||
|
|
||||||
H 0.0 0.0 0.0
|
H 0.0 0.0 0.0
|
||||||
H 0.0 0.0 1.0
|
H 0.0 0.0 0.740848
|
||||||
|
|||||||
231
src/CC/BCCD.f90
Normal file
231
src/CC/BCCD.f90
Normal file
@ -0,0 +1,231 @@
|
|||||||
|
subroutine BCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
|
||||||
|
|
||||||
|
! Brueckner CCD module
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
integer,intent(in) :: maxSCF
|
||||||
|
integer,intent(in) :: max_diis
|
||||||
|
double precision,intent(in) :: thresh
|
||||||
|
|
||||||
|
integer,intent(in) :: nBasin
|
||||||
|
integer,intent(in) :: nCin
|
||||||
|
integer,intent(in) :: nOin
|
||||||
|
integer,intent(in) :: nVin
|
||||||
|
integer,intent(in) :: nRin
|
||||||
|
double precision,intent(in) :: ENuc,ERHF
|
||||||
|
double precision,intent(in) :: eHF(nBasin)
|
||||||
|
double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
integer :: nBas
|
||||||
|
integer :: nC
|
||||||
|
integer :: nO
|
||||||
|
integer :: nV
|
||||||
|
integer :: nR
|
||||||
|
integer :: nSCF
|
||||||
|
double precision :: Conv
|
||||||
|
double precision :: EcMP2,EcMP3,EcMP4
|
||||||
|
double precision :: ECCD,EcCCD
|
||||||
|
double precision,allocatable :: seHF(:)
|
||||||
|
double precision,allocatable :: sERI(:,:,:,:)
|
||||||
|
double precision,allocatable :: dbERI(:,:,:,:)
|
||||||
|
|
||||||
|
double precision,allocatable :: eO(:)
|
||||||
|
double precision,allocatable :: eV(:)
|
||||||
|
double precision,allocatable :: delta_OOVV(:,:,:,:)
|
||||||
|
|
||||||
|
double precision,allocatable :: OOOO(:,:,:,:)
|
||||||
|
double precision,allocatable :: OOVV(:,:,:,:)
|
||||||
|
double precision,allocatable :: OVOV(:,:,:,:)
|
||||||
|
double precision,allocatable :: VVVV(:,:,:,:)
|
||||||
|
|
||||||
|
double precision,allocatable :: X1(:,:,:,:)
|
||||||
|
double precision,allocatable :: X2(:,:)
|
||||||
|
double precision,allocatable :: X3(:,:)
|
||||||
|
double precision,allocatable :: X4(:,:,:,:)
|
||||||
|
|
||||||
|
double precision,allocatable :: u(:,:,:,:)
|
||||||
|
double precision,allocatable :: v(:,:,:,:)
|
||||||
|
|
||||||
|
double precision,allocatable :: r2(:,:,:,:)
|
||||||
|
double precision,allocatable :: t2(:,:,:,:)
|
||||||
|
|
||||||
|
integer :: n_diis,i,j,a,b
|
||||||
|
double precision :: rcond
|
||||||
|
double precision,allocatable :: error_diis(:,:)
|
||||||
|
double precision,allocatable :: t_diis(:,:)
|
||||||
|
|
||||||
|
! Hello world
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'**************************************'
|
||||||
|
write(*,*)'| BCCD calculation |'
|
||||||
|
write(*,*)'**************************************'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
! Spatial to spin orbitals
|
||||||
|
|
||||||
|
nBas = 2*nBasin
|
||||||
|
nC = 2*nCin
|
||||||
|
nO = 2*nOin
|
||||||
|
nV = 2*nVin
|
||||||
|
nR = 2*nRin
|
||||||
|
|
||||||
|
allocate(seHF(nBas),sERI(nBas,nBas,nBas,nBas))
|
||||||
|
|
||||||
|
call spatial_to_spin_MO_energy(nBasin,eHF,nBas,seHF)
|
||||||
|
call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI)
|
||||||
|
|
||||||
|
! Antysymmetrize ERIs
|
||||||
|
|
||||||
|
allocate(dbERI(nBas,nBas,nBas,nBas))
|
||||||
|
|
||||||
|
call antisymmetrize_ERI(2,nBas,sERI,dbERI)
|
||||||
|
|
||||||
|
deallocate(sERI)
|
||||||
|
|
||||||
|
! Form energy denominator
|
||||||
|
|
||||||
|
allocate(eO(nO-nC),eV(nV-nR))
|
||||||
|
allocate(delta_OOVV(nO-nC,nO-nC,nV-nR,nV-nR))
|
||||||
|
|
||||||
|
eO(:) = seHF(nC+1:nO)
|
||||||
|
eV(:) = seHF(nO+1:nBas-nR)
|
||||||
|
|
||||||
|
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||||
|
|
||||||
|
deallocate(seHF)
|
||||||
|
|
||||||
|
! Create integral batches
|
||||||
|
|
||||||
|
allocate(OOOO(nO-nC,nO-nC,nO-nC,nO-nC),OOVV(nO-nC,nO-nC,nV-nR,nV-nR), &
|
||||||
|
OVOV(nO-nC,nV-nR,nO-nC,nV-nR),VVVV(nV-nR,nV-nR,nV-nR,nV-nR))
|
||||||
|
|
||||||
|
OOOO(:,:,:,:) = dbERI(nC+1:nO ,nC+1:nO ,nC+1:nO ,nC+1:nO )
|
||||||
|
OOVV(:,:,:,:) = dbERI(nC+1:nO ,nC+1:nO ,nO+1:nBas-nR,nO+1:nBas-nR)
|
||||||
|
OVOV(:,:,:,:) = dbERI(nC+1:nO ,nO+1:nBas-nR,nC+1:nO ,nO+1:nBas-nR)
|
||||||
|
VVVV(:,:,:,:) = dbERI(nO+1:nBas-nR,nO+1:nBas-nR,nO+1:nBas-nR,nO+1:nBas-nR)
|
||||||
|
|
||||||
|
deallocate(dbERI)
|
||||||
|
|
||||||
|
! MP2 guess amplitudes
|
||||||
|
|
||||||
|
allocate(t2(nO-nC,nO-nC,nV-nR,nV-nR))
|
||||||
|
|
||||||
|
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||||
|
|
||||||
|
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcMP2)
|
||||||
|
EcMP4 = 0d0
|
||||||
|
|
||||||
|
! Memory allocation for DIIS
|
||||||
|
|
||||||
|
allocate(error_diis((nO-nR)**2*(nV-nR)**2,max_diis),t_diis((nO-nR)**2*(nV-nR)**2,max_diis))
|
||||||
|
|
||||||
|
! Initialization
|
||||||
|
|
||||||
|
allocate(r2(nO-nC,nO-nC,nV-nR,nV-nR),u(nO-nC,nO-nC,nV-nR,nV-nR),v(nO-nC,nO-nC,nV-nR,nV-nR))
|
||||||
|
allocate(X1(nO-nC,nO-nC,nO-nC,nO-nC),X2(nV-nR,nV-nR),X3(nO-nC,nO-nC),X4(nO-nC,nO-nC,nV-nR,nV-nR))
|
||||||
|
|
||||||
|
Conv = 1d0
|
||||||
|
nSCF = 0
|
||||||
|
|
||||||
|
n_diis = 0
|
||||||
|
t_diis(:,:) = 0d0
|
||||||
|
error_diis(:,:) = 0d0
|
||||||
|
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
! Main SCF loop
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'----------------------------------------------------'
|
||||||
|
write(*,*)'| BCCD calculation |'
|
||||||
|
write(*,*)'----------------------------------------------------'
|
||||||
|
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
|
||||||
|
'|','#','|','E(BCCD)','|','Ec(BCCD)','|','Conv','|'
|
||||||
|
write(*,*)'----------------------------------------------------'
|
||||||
|
|
||||||
|
do while(Conv > thresh .and. nSCF < maxSCF)
|
||||||
|
|
||||||
|
! Increment
|
||||||
|
|
||||||
|
nSCF = nSCF + 1
|
||||||
|
|
||||||
|
! Form linear array
|
||||||
|
|
||||||
|
call form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t2,u)
|
||||||
|
|
||||||
|
! Form interemediate arrays
|
||||||
|
|
||||||
|
call form_X(nC,nO,nV,nR,OOVV,t2,X1,X2,X3,X4)
|
||||||
|
|
||||||
|
! Form quadratic array
|
||||||
|
|
||||||
|
call form_v(nC,nO,nV,nR,X1,X2,X3,X4,t2,v)
|
||||||
|
|
||||||
|
! Compute residual
|
||||||
|
|
||||||
|
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + u(:,:,:,:) + v(:,:,:,:)
|
||||||
|
|
||||||
|
! Check convergence
|
||||||
|
|
||||||
|
Conv = maxval(abs(r2(:,:,:,:)))
|
||||||
|
|
||||||
|
! Update amplitudes
|
||||||
|
|
||||||
|
t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||||
|
|
||||||
|
! Compute correlation energy
|
||||||
|
|
||||||
|
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCCD)
|
||||||
|
|
||||||
|
if(nSCF == 1) call MP3_correlation_energy(nC,nO,nV,nR,OOVV,t2,v,delta_OOVV,EcMP3)
|
||||||
|
|
||||||
|
! Dump results
|
||||||
|
|
||||||
|
ECCD = ERHF + EcCCD
|
||||||
|
|
||||||
|
! DIIS extrapolation
|
||||||
|
|
||||||
|
n_diis = min(n_diis+1,max_diis)
|
||||||
|
call DIIS_extrapolation(rcond,(nO-nC)**2*(nV-nR)**2,(nO-nC)**2*(nV-nR)**2,n_diis,error_diis,t_diis,-r2/delta_OOVV,t2)
|
||||||
|
|
||||||
|
! Reset DIIS if required
|
||||||
|
|
||||||
|
if(abs(rcond) < 1d-15) n_diis = 0
|
||||||
|
|
||||||
|
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
|
||||||
|
'|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|'
|
||||||
|
|
||||||
|
enddo
|
||||||
|
write(*,*)'----------------------------------------------------'
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
! End of SCF loop
|
||||||
|
!------------------------------------------------------------------------
|
||||||
|
|
||||||
|
! Did it actually converge?
|
||||||
|
|
||||||
|
if(nSCF == maxSCF) then
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||||
|
write(*,*)' Convergence failed '
|
||||||
|
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
stop
|
||||||
|
|
||||||
|
endif
|
||||||
|
|
||||||
|
! Moller-Plesset energies
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP2) = ',EcMP2
|
||||||
|
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP3) = ',EcMP3
|
||||||
|
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP4-SDQ) = ',EcMP4
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
end subroutine BCCD
|
||||||
@ -97,7 +97,6 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
|||||||
EcCCD = EcCCD + OOVV(i,a)*t(i,a)
|
EcCCD = EcCCD + OOVV(i,a)*t(i,a)
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
print*,'Ec = ',EcCCD
|
|
||||||
|
|
||||||
! Memory allocation for DIIS
|
! Memory allocation for DIIS
|
||||||
|
|
||||||
|
|||||||
@ -79,9 +79,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,sing
|
|||||||
nOOs = nO*nO
|
nOOs = nO*nO
|
||||||
nVVs = nV*nV
|
nVVs = nV*nV
|
||||||
|
|
||||||
! nOOs = nO*(nO + 1)/2
|
|
||||||
! nVVs = nV*(nV + 1)/2
|
|
||||||
|
|
||||||
nOOt = nO*(nO - 1)/2
|
nOOt = nO*(nO - 1)/2
|
||||||
nVVt = nV*(nV - 1)/2
|
nVVt = nV*(nV - 1)/2
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user