mirror of
https://github.com/pfloos/quack
synced 2025-01-03 10:05:59 +01:00
adding routines for MBPT vs CC
This commit is contained in:
parent
69bab74d81
commit
96621bd038
@ -3,13 +3,13 @@
|
||||
# MP2* MP3 MP2-F12
|
||||
F F F
|
||||
# CCD pCCD DCD CCSD CCSD(T)
|
||||
T F F F F
|
||||
F F F F F
|
||||
# drCCD rCCD crCCD lCCD
|
||||
F F F F
|
||||
F T F F
|
||||
# CIS* CIS(D) CID CISD FCI
|
||||
F F F F F
|
||||
# RPA* RPAx* crRPA ppRPA
|
||||
F T F F
|
||||
F F F F
|
||||
# G0F2* evGF2* qsGF2* G0F3 evGF3
|
||||
F F F F F
|
||||
# G0W0* evGW* qsGW* ufG0W0 ufGW
|
||||
|
||||
@ -3,9 +3,9 @@
|
||||
# MP:
|
||||
|
||||
# CC: maxSCF thresh DIIS n_diis
|
||||
64 0.00001 T 5
|
||||
64 0.0000001 T 5
|
||||
# spin: TDA singlet triplet spin_conserved spin_flip
|
||||
T T T T T
|
||||
F T T T T
|
||||
# GF: maxSCF thresh DIIS n_diis lin eta renorm reg
|
||||
256 0.00001 T 5 T 0.0 3 F
|
||||
# GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 reg
|
||||
@ -13,8 +13,8 @@
|
||||
# GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg
|
||||
10 0.00001 T 5 T 0.0 T F
|
||||
# ACFDT: AC Kx XBS
|
||||
F F T
|
||||
F T T
|
||||
# BSE: BSE dBSE dTDA evDyn ppBSE
|
||||
T F T F F
|
||||
F F T F F
|
||||
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
|
||||
1000000 100000 10 0.3 10000 1234 T
|
||||
|
||||
@ -62,7 +62,7 @@ subroutine CCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ER
|
||||
double precision,allocatable :: error_diis(:,:)
|
||||
double precision,allocatable :: t_diis(:,:)
|
||||
|
||||
logical :: do_EE_EOM_CC_1h1p = .true.
|
||||
logical :: do_EE_EOM_CC_1h1p = .false.
|
||||
logical :: do_EA_EOM_CC_1p = .false.
|
||||
logical :: do_IP_EOM_CC_1h = .false.
|
||||
logical :: do_DEA_EOM_CC_2p = .false.
|
||||
|
||||
@ -28,7 +28,9 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
|
||||
double precision,allocatable :: Wovvo(:,:,:,:)
|
||||
double precision,allocatable :: H(:,:)
|
||||
double precision,allocatable :: Om(:)
|
||||
double precision,allocatable :: Z(:,:)
|
||||
|
||||
integer,allocatable :: order(:)
|
||||
|
||||
! Hello world
|
||||
|
||||
@ -44,7 +46,9 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(Foo(nO,nO),Fvv(nV,nV),Wovvo(nO,nV,nV,nO),H(nS,nS),Om(nS))
|
||||
allocate(Foo(nO,nO),Fvv(nV,nV),Wovvo(nO,nV,nV,nO),H(nS,nS),Om(nS),Z(nS,nS))
|
||||
allocate(order(nS))
|
||||
|
||||
|
||||
! Form one-body terms
|
||||
|
||||
@ -57,7 +61,7 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
|
||||
do j=1,nO-nC
|
||||
do c=1,nV-nR
|
||||
|
||||
Fvv(a,b) = Fvv(a,b) - 0.5d0*OOVV(i,j,b,c)*t(i,j,a,c)
|
||||
! Fvv(a,b) = Fvv(a,b) - 0.5d0*OOVV(i,j,b,c)*t(i,j,a,c)
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -75,7 +79,7 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
Foo(i,j) = Foo(i,j) + 0.5d0*OOVV(i,k,a,b)*t(j,k,a,b)
|
||||
! Foo(i,j) = Foo(i,j) + 0.5d0*OOVV(i,k,a,b)*t(j,k,a,b)
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -128,10 +132,19 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
|
||||
|
||||
! Diagonalize EOM Hamiltonian
|
||||
|
||||
if(nS > 0) call diagonalize_matrix(nS,H,Om)
|
||||
if(nS > 0) then
|
||||
|
||||
! Dump results
|
||||
call diagonalize_general_matrix(nS,H,Om,Z)
|
||||
|
||||
call print_excitation('EE-EOM-CCD ',3,nS,Om)
|
||||
do ia=1,nS
|
||||
order(ia) = ia
|
||||
end do
|
||||
|
||||
call quick_sort(Om,order,nS)
|
||||
call set_order(Z,order,nS,nS)
|
||||
|
||||
call print_excitation('EE-EOM-CCD ',3,nS,Om)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine EE_EOM_CCD_1h1p
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF)
|
||||
subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF,eGW)
|
||||
|
||||
! Ring CCD module
|
||||
|
||||
@ -16,8 +16,10 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E
|
||||
integer,intent(in) :: nOin
|
||||
integer,intent(in) :: nVin
|
||||
integer,intent(in) :: nRin
|
||||
double precision,intent(in) :: ENuc,ERHF
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: ERHF
|
||||
double precision,intent(in) :: eHF(nBasin)
|
||||
double precision,intent(in) :: eGW(nBasin)
|
||||
double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin)
|
||||
|
||||
! Local variables
|
||||
@ -32,6 +34,7 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E
|
||||
double precision :: EcMP2
|
||||
double precision :: ECCD,EcCCD
|
||||
double precision,allocatable :: seHF(:)
|
||||
double precision,allocatable :: seGW(:)
|
||||
double precision,allocatable :: sERI(:,:,:,:)
|
||||
double precision,allocatable :: dbERI(:,:,:,:)
|
||||
|
||||
@ -42,14 +45,16 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E
|
||||
double precision,allocatable :: OOVV(:,:,:,:)
|
||||
double precision,allocatable :: OVVO(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: r2(:,:,:,:)
|
||||
double precision,allocatable :: t2(:,:,:,:)
|
||||
double precision,allocatable :: r(:,:,:,:)
|
||||
double precision,allocatable :: t(:,:,:,:)
|
||||
|
||||
integer :: n_diis
|
||||
double precision :: rcond
|
||||
double precision,allocatable :: error_diis(:,:)
|
||||
double precision,allocatable :: t_diis(:,:)
|
||||
|
||||
logical :: do_EE_EOM_CC_1h1p = .true.
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
@ -66,11 +71,16 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E
|
||||
nV = 2*nVin
|
||||
nR = 2*nRin
|
||||
|
||||
allocate(seHF(nBas),sERI(nBas,nBas,nBas,nBas))
|
||||
allocate(seHF(nBas),seGW(nBas),sERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call spatial_to_spin_MO_energy(nBasin,eHF,nBas,seHF)
|
||||
call spatial_to_spin_MO_energy(nBasin,eGW,nBas,seGW)
|
||||
call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI)
|
||||
|
||||
! call soRPAx(.false.,nBas,nC,nO,nV,nR,nO*nV,ENuc,ERHF,sERI,seHF)
|
||||
|
||||
call soBSE(.false.,.false.,0.0,nBas,nC,nO,nV,nR,nO*nV,sERI,seHF,seGW)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
allocate(dbERI(nBas,nBas,nBas,nBas))
|
||||
@ -92,12 +102,12 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E
|
||||
allocate(eO(nO),eV(nV))
|
||||
allocate(delta_OOVV(nO,nO,nV,nV))
|
||||
|
||||
eO(:) = seHF(1:nO)
|
||||
eV(:) = seHF(nO+1:nBas)
|
||||
eO(:) = seGW(1:nO)
|
||||
eV(:) = seGW(nO+1:nBas)
|
||||
|
||||
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta_OOVV)
|
||||
|
||||
deallocate(seHF)
|
||||
! deallocate(seHF,seGW)
|
||||
|
||||
! Create integral batches
|
||||
|
||||
@ -110,11 +120,11 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E
|
||||
|
||||
! MP2 guess amplitudes
|
||||
|
||||
allocate(t2(nO,nO,nV,nV))
|
||||
allocate(t(nO,nO,nV,nV))
|
||||
|
||||
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
t(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcMP2)
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcMP2)
|
||||
|
||||
! Memory allocation for DIIS
|
||||
|
||||
@ -122,7 +132,7 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E
|
||||
|
||||
! Initialization
|
||||
|
||||
allocate(r2(nO,nO,nV,nV))
|
||||
allocate(r(nO,nO,nV,nV))
|
||||
|
||||
Conv = 1d0
|
||||
nSCF = 0
|
||||
@ -150,21 +160,21 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E
|
||||
|
||||
! Compute residual
|
||||
|
||||
call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
|
||||
call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t,r)
|
||||
|
||||
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
|
||||
r(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t(:,:,:,:) + r(:,:,:,:)
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = maxval(abs(r2(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR)))
|
||||
Conv = maxval(abs(r(nC+1:nO,nC+1:nO,1:nV-nR,1:nV-nR)))
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
t(:,:,:,:) = t(:,:,:,:) - r(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
! Compute correlation energy
|
||||
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCCD)
|
||||
call CCD_correlation_energy(nC,nO,nV,nR,OOVV,t,EcCCD)
|
||||
|
||||
! Dump results
|
||||
|
||||
@ -173,7 +183,7 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E
|
||||
! DIIS extrapolation
|
||||
|
||||
n_diis = min(n_diis+1,max_diis)
|
||||
call DIIS_extrapolation(rcond,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,error_diis,t_diis,-r2/delta_OOVV,t2)
|
||||
call DIIS_extrapolation(rcond,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,error_diis,t_diis,-r/delta_OOVV,t)
|
||||
|
||||
! Reset DIIS if required
|
||||
|
||||
@ -211,4 +221,19 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*)'----------------------------------------------------'
|
||||
! write(*,*)' ring CCD amplitudes '
|
||||
! write(*,*)'----------------------------------------------------'
|
||||
! call matout(nO*nO,nV*nV,t)
|
||||
! write(*,*)
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! EOM section
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! EE-EOM-CCD (1h1p)
|
||||
|
||||
if(do_EE_EOM_CC_1h1p) call EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
|
||||
|
||||
end subroutine rCCD
|
||||
|
||||
122
src/GW/soBSE.f90
Normal file
122
src/GW/soBSE.f90
Normal file
@ -0,0 +1,122 @@
|
||||
subroutine soBSE(TDA_W,TDA,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,eGW)
|
||||
|
||||
! Compute the Bethe-Salpeter excitation energies
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: TDA_W
|
||||
logical,intent(in) :: TDA
|
||||
|
||||
double precision,intent(in) :: eta
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
double precision,intent(in) :: eW(nBas)
|
||||
double precision,intent(in) :: eGW(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: ispin
|
||||
integer :: isp_W
|
||||
|
||||
double precision :: EcRPA
|
||||
double precision,allocatable :: OmRPA(:)
|
||||
double precision,allocatable :: XpY_RPA(:,:)
|
||||
double precision,allocatable :: XmY_RPA(:,:)
|
||||
double precision,allocatable :: rho_RPA(:,:,:)
|
||||
|
||||
double precision,allocatable :: OmBSE(:)
|
||||
double precision,allocatable :: XpY_BSE(:,:)
|
||||
double precision,allocatable :: XmY_BSE(:,:)
|
||||
|
||||
double precision,allocatable :: WA_sta(:,:)
|
||||
double precision,allocatable :: WB_sta(:,:)
|
||||
|
||||
double precision,allocatable :: X(:,:)
|
||||
double precision,allocatable :: Y(:,:)
|
||||
double precision,allocatable :: Xinv(:,:)
|
||||
double precision,allocatable :: t(:,:,:,:)
|
||||
|
||||
integer ::i,a,j,b,k,c,ia,jb,kc
|
||||
|
||||
double precision :: EcBSE
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), &
|
||||
WA_sta(nS,nS),WB_sta(nS,nS),OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS))
|
||||
|
||||
!---------------------------------
|
||||
! Compute (singlet) RPA screening
|
||||
!---------------------------------
|
||||
|
||||
isp_W = 3
|
||||
EcRPA = 0d0
|
||||
|
||||
call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
|
||||
EcRPA,OmRPA,XpY_RPA,XmY_RPA)
|
||||
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
|
||||
|
||||
call static_screening_WA_so(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA_sta)
|
||||
call static_screening_WB_so(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB_sta)
|
||||
|
||||
|
||||
ispin = 3
|
||||
EcBSE = 0d0
|
||||
|
||||
! Compute BSE excitation energies
|
||||
|
||||
call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,WA_sta,WB_sta, &
|
||||
EcBSE,OmBSE,XpY_BSE,XmY_BSE)
|
||||
call print_excitation('soBSE@GW ',ispin,nS,OmBSE)
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy =',0.5d0*EcBSE,' au'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
! allocate(X(nS,nS),Y(nS,nS),Xinv(nS,nS),t(nO,nO,nV,nV))
|
||||
|
||||
! X(:,:) = transpose(0.5d0*(XpY_BSE(:,:) + XmY_BSE(:,:)))
|
||||
! Y(:,:) = transpose(0.5d0*(XpY_BSE(:,:) - XmY_BSE(:,:)))
|
||||
|
||||
! call matout(nS,nS,matmul(X,transpose(X))-matmul(Y,transpose(Y)))
|
||||
|
||||
! call inverse_matrix(nS,X,Xinv)
|
||||
|
||||
! t = 0d0
|
||||
! ia = 0
|
||||
! do i=1,nO
|
||||
! do a=1,nV
|
||||
! ia = ia + 1
|
||||
|
||||
! jb = 0
|
||||
! do j=1,nO
|
||||
! do b=1,nV
|
||||
! jb = jb + 1
|
||||
|
||||
! kc = 0
|
||||
! do k=1,nO
|
||||
! do c=1,nV
|
||||
! kc = kc + 1
|
||||
|
||||
! t(i,j,a,b) = t(i,j,a,b) + Y(ia,kc)*Xinv(kc,jb)
|
||||
|
||||
! end do
|
||||
! end do
|
||||
! end do
|
||||
! end do
|
||||
! end do
|
||||
! end do
|
||||
|
||||
! call matout(nO*nO,nV*nV,t)
|
||||
|
||||
end subroutine soBSE
|
||||
53
src/GW/static_screening_WA_so.f90
Normal file
53
src/GW/static_screening_WA_so.f90
Normal file
@ -0,0 +1,53 @@
|
||||
subroutine static_screening_WA_so(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WA)
|
||||
|
||||
! Compute the OOVV block of the static screening W for the resonant block
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
|
||||
double precision,intent(in) :: eta
|
||||
double precision,intent(in) :: lambda
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: Omega(nS)
|
||||
double precision,intent(in) :: rho(nBas,nBas,nS)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: chi
|
||||
double precision :: eps
|
||||
integer :: i,j,a,b,ia,jb,kc
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: WA(nS,nS)
|
||||
|
||||
! Initialize
|
||||
|
||||
WA(:,:) = 0d0
|
||||
|
||||
ia = 0
|
||||
do i=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
ia = ia + 1
|
||||
jb = 0
|
||||
do j=nC+1,nO
|
||||
do b=nO+1,nBas-nR
|
||||
jb = jb + 1
|
||||
|
||||
chi = 0d0
|
||||
do kc=1,nS
|
||||
eps = Omega(kc)**2 + eta**2
|
||||
chi = chi + rho(i,j,kc)*rho(a,b,kc)*Omega(kc)/eps
|
||||
enddo
|
||||
|
||||
WA(ia,jb) = WA(ia,jb) + lambda*ERI(i,b,j,a) - 2d0*lambda*chi
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine static_screening_WA_so
|
||||
58
src/GW/static_screening_WB_so.f90
Normal file
58
src/GW/static_screening_WB_so.f90
Normal file
@ -0,0 +1,58 @@
|
||||
subroutine static_screening_WB_so(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WB)
|
||||
|
||||
! Compute the static screening W for the coupling block
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
double precision,intent(in) :: eta
|
||||
double precision,intent(in) :: lambda
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: Omega(nS)
|
||||
double precision,intent(in) :: rho(nBas,nBas,nS)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: chi
|
||||
double precision :: eps
|
||||
integer :: i,j,a,b,ia,jb,kc
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: WB(nS,nS)
|
||||
|
||||
! Initialize
|
||||
|
||||
WB(:,:) = 0d0
|
||||
|
||||
ia = 0
|
||||
do i=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
ia = ia + 1
|
||||
jb = 0
|
||||
do j=nC+1,nO
|
||||
do b=nO+1,nBas-nR
|
||||
jb = jb + 1
|
||||
|
||||
chi = 0d0
|
||||
do kc=1,nS
|
||||
eps = Omega(kc)**2 + eta**2
|
||||
chi = chi + rho(i,b,kc)*rho(a,j,kc)*Omega(kc)/eps
|
||||
enddo
|
||||
|
||||
WB(ia,jb) = WB(ia,jb) + lambda*ERI(i,j,b,a) - 2d0*lambda*chi
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine static_screening_WB_so
|
||||
@ -662,21 +662,22 @@ program QuAcK
|
||||
|
||||
if(do_rCCD) then
|
||||
|
||||
call cpu_time(start_CCD)
|
||||
call rCCD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call cpu_time(end_CCD)
|
||||
! call cpu_time(start_CCD)
|
||||
! call rCCD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,eHF)
|
||||
! call cpu_time(end_CCD)
|
||||
|
||||
! t_CCD = end_CCD - start_CCD
|
||||
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD = ',t_CCD,' seconds'
|
||||
! write(*,*)
|
||||
|
||||
t_CCD = end_CCD - start_CCD
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ring CCD = ',t_CCD,' seconds'
|
||||
write(*,*)
|
||||
call cpu_time(start_CCD)
|
||||
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, &
|
||||
linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0)
|
||||
call rCCD(.true.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eG0W0)
|
||||
call rCCD(.true.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,eG0W0)
|
||||
call cpu_time(end_CCD)
|
||||
|
||||
t_CCD = end_CCD - start_CCD
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ring CCD = ',t_CCD,' seconds'
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD@BSE = ',t_CCD,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
@ -829,7 +830,7 @@ program QuAcK
|
||||
|
||||
else
|
||||
|
||||
call RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
|
||||
call RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
|
||||
|
||||
end if
|
||||
call cpu_time(end_RPA)
|
||||
|
||||
113
src/RPA/soRPAx.f90
Normal file
113
src/RPA/soRPAx.f90
Normal file
@ -0,0 +1,113 @@
|
||||
subroutine soRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
|
||||
|
||||
! Perform random phase approximation calculation with exchange (aka TDHF) in the
|
||||
! spinorbital basis
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
include 'quadrature.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: TDA
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
integer,intent(in) :: nS
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: ispin
|
||||
double precision,allocatable :: Omega(:)
|
||||
double precision,allocatable :: XpY(:,:)
|
||||
double precision,allocatable :: XmY(:,:)
|
||||
double precision,allocatable :: X(:,:)
|
||||
double precision,allocatable :: Y(:,:)
|
||||
double precision,allocatable :: Xinv(:,:)
|
||||
double precision,allocatable :: t(:,:,:,:)
|
||||
|
||||
double precision :: EcRPAx
|
||||
|
||||
integer ::i,a,j,b,k,c,ia,jb,kc
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'***********************************************************'
|
||||
write(*,*)'| Random phase approximation calculation with exchange |'
|
||||
write(*,*)'***********************************************************'
|
||||
write(*,*)
|
||||
|
||||
! TDA
|
||||
|
||||
if(TDA) then
|
||||
write(*,*) 'Tamm-Dancoff approximation activated!'
|
||||
write(*,*) ' => RPAx + TDA = CIS '
|
||||
write(*,*)
|
||||
end if
|
||||
|
||||
! Initialization
|
||||
|
||||
EcRPAx = 0d0
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(Omega(nS),XpY(nS,nS),XmY(nS,nS))
|
||||
|
||||
ispin = 3
|
||||
|
||||
call linear_response(ispin,.false.,TDA,0d0,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,EcRPAx,Omega,XpY,XmY)
|
||||
call print_excitation('soRPAx@HF ',ispin,nS,Omega)
|
||||
|
||||
EcRPAx = 0.5d0*EcRPAx
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@RPAx correlation energy =',EcRPAx
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@RPAx total energy =',ENuc + ERHF + EcRPAx
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
! allocate(X(nS,nS),Y(nS,nS),Xinv(nS,nS),t(nO,nO,nV,nV))
|
||||
|
||||
! X(:,:) = transpose(0.5d0*(XpY(:,:) + XmY(:,:)))
|
||||
! Y(:,:) = transpose(0.5d0*(XpY(:,:) - XmY(:,:)))
|
||||
|
||||
! call matout(nS,nS,matmul(transpose(X),X)-matmul(transpose(Y),Y))
|
||||
|
||||
! call inverse_matrix(nS,X,Xinv)
|
||||
|
||||
! t = 0d0
|
||||
! ia = 0
|
||||
! do i=1,nO
|
||||
! do a=1,nV
|
||||
! ia = ia + 1
|
||||
|
||||
! jb = 0
|
||||
! do j=1,nO
|
||||
! do b=1,nV
|
||||
! jb = jb + 1
|
||||
|
||||
! kc = 0
|
||||
! do k=1,nO
|
||||
! do c=1,nV
|
||||
! kc = kc + 1
|
||||
|
||||
! t(i,j,a,b) = t(i,j,a,b) + Y(ia,kc)*Xinv(kc,jb)
|
||||
|
||||
! end do
|
||||
! end do
|
||||
! end do
|
||||
! end do
|
||||
! end do
|
||||
! end do
|
||||
|
||||
! call matout(nO*nO,nV*nV,t)
|
||||
|
||||
end subroutine soRPAx
|
||||
@ -15,8 +15,8 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
|
||||
! Local variables
|
||||
|
||||
integer :: pq,ab,ij
|
||||
integer :: deg1,ab_start,ab_end
|
||||
integer :: deg2,ij_start,ij_end
|
||||
! integer :: deg1,ab_start,ab_end
|
||||
! integer :: deg2,ij_start,ij_end
|
||||
double precision,allocatable :: M(:,:)
|
||||
double precision,allocatable :: Z1(:,:)
|
||||
double precision,allocatable :: Z2(:,:)
|
||||
|
||||
Loading…
Reference in New Issue
Block a user