mirror of
https://github.com/pfloos/quack
synced 2025-01-03 10:05:59 +01:00
post-GHF methods
This commit is contained in:
Pierre-Francois Loos
parent
7578468081
commit
5b8c0d542c
@ -1,5 +1,5 @@
|
||||
# RHF UHF GHF ROHF
|
||||
F F T F
|
||||
F T F F
|
||||
# MP2* MP3
|
||||
F F
|
||||
# CCD pCCD DCD CCSD CCSD(T)
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
# HF: maxSCF thresh DIIS guess mix_guess level_shift stability
|
||||
5000 0.0000001 5 1 0.0 0.0 T
|
||||
10000 0.000001 5 2 0.0 0.0 T
|
||||
# MP: reg
|
||||
F
|
||||
# CC: maxSCF thresh DIIS
|
||||
|
||||
90
src/AOtoMO/AOtoMO_integral_transform_GHF.f90
Normal file
90
src/AOtoMO/AOtoMO_integral_transform_GHF.f90
Normal file
@ -0,0 +1,90 @@
|
||||
subroutine AOtoMO_integral_transform_GHF(nBas,nBas2,c1,c2,c3,c4,ERI_AO_basis,ERI_MO_basis)
|
||||
|
||||
! AO to MO transformation of two-electron integrals via the semi-direct O(N^5) algorithm
|
||||
! bra and ket are the spin of (bra1 bra2|ket1 ket2)
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nBas2
|
||||
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: c1(nBas,nBas2)
|
||||
double precision,intent(in) :: c2(nBas,nBas2)
|
||||
double precision,intent(in) :: c3(nBas,nBas2)
|
||||
double precision,intent(in) :: c4(nBas,nBas2)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision,allocatable :: scr(:,:,:,:)
|
||||
integer :: mu,nu,la,si,i,j,k,l
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: ERI_MO_basis(nBas2,nBas2,nBas2,nBas2)
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(scr(nBas2,nBas2,nBas2,nBas2))
|
||||
|
||||
! Four-index transform via semi-direct O(N^5) algorithm
|
||||
|
||||
scr(:,:,:,:) = 0d0
|
||||
|
||||
do l=1,nBas2
|
||||
do si=1,nBas
|
||||
do la=1,nBas
|
||||
do nu=1,nBas
|
||||
do mu=1,nBas
|
||||
scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c4(si,l)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
ERI_MO_basis(:,:,:,:) = 0d0
|
||||
|
||||
do l=1,nBas2
|
||||
do la=1,nBas
|
||||
do nu=1,nBas
|
||||
do i=1,nBas2
|
||||
do mu=1,nBas
|
||||
ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c1(mu,i)*scr(mu,nu,la,l)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
scr(:,:,:,:) = 0d0
|
||||
|
||||
do l=1,nBas2
|
||||
do k=1,nBas2
|
||||
do la=1,nBas
|
||||
do nu=1,nBas
|
||||
do i=1,nBas2
|
||||
scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c3(la,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
ERI_MO_basis(:,:,:,:) = 0d0
|
||||
|
||||
do l=1,nBas2
|
||||
do k=1,nBas2
|
||||
do j=1,nBas2
|
||||
do i=1,nBas2
|
||||
do nu=1,nBas
|
||||
ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c2(nu,j)*scr(i,nu,k,l)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine
|
||||
231
src/GW/GG0W0.f90
Normal file
231
src/GW/GG0W0.f90
Normal file
@ -0,0 +1,231 @@
|
||||
subroutine GG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
|
||||
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
|
||||
! Perform G0W0 calculation
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
include 'quadrature.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: doACFDT
|
||||
logical,intent(in) :: exchange_kernel
|
||||
logical,intent(in) :: doXBS
|
||||
logical,intent(in) :: dophBSE
|
||||
logical,intent(in) :: dophBSE2
|
||||
logical,intent(in) :: doppBSE
|
||||
logical,intent(in) :: TDA_W
|
||||
logical,intent(in) :: TDA
|
||||
logical,intent(in) :: dBSE
|
||||
logical,intent(in) :: dTDA
|
||||
logical,intent(in) :: linearize
|
||||
double precision,intent(in) :: eta
|
||||
logical,intent(in) :: regularize
|
||||
|
||||
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) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
logical :: print_W = .true.
|
||||
logical :: dRPA
|
||||
integer :: ispin
|
||||
double precision :: EcRPA
|
||||
double precision :: EcBSE(nspin)
|
||||
double precision :: EcGM
|
||||
double precision,allocatable :: Aph(:,:)
|
||||
double precision,allocatable :: Bph(:,:)
|
||||
double precision,allocatable :: SigC(:)
|
||||
double precision,allocatable :: Z(:)
|
||||
double precision,allocatable :: Om(:)
|
||||
double precision,allocatable :: XpY(:,:)
|
||||
double precision,allocatable :: XmY(:,:)
|
||||
double precision,allocatable :: rho(:,:,:)
|
||||
|
||||
double precision,allocatable :: eGWlin(:)
|
||||
double precision,allocatable :: eGW(:)
|
||||
|
||||
! Output variables
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'************************************************'
|
||||
write(*,*)'| One-shot G0W0 calculation |'
|
||||
write(*,*)'************************************************'
|
||||
write(*,*)
|
||||
|
||||
! Initialization
|
||||
|
||||
dRPA = .true.
|
||||
EcRPA = 0d0
|
||||
|
||||
! TDA for W
|
||||
|
||||
if(TDA_W) then
|
||||
write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
|
||||
write(*,*)
|
||||
end if
|
||||
|
||||
! TDA
|
||||
|
||||
if(TDA) then
|
||||
write(*,*) 'Tamm-Dancoff approximation activated!'
|
||||
write(*,*)
|
||||
end if
|
||||
|
||||
! Spin manifold
|
||||
|
||||
ispin = 3
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(Aph(nS,nS),Bph(nS,nS),SigC(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), &
|
||||
eGW(nBas),eGWlin(nBas))
|
||||
|
||||
!-------------------!
|
||||
! Compute screening !
|
||||
!-------------------!
|
||||
|
||||
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
|
||||
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||
|
||||
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||
|
||||
if(print_W) call print_excitation_energies('phRPA@GHF',ispin,nS,Om)
|
||||
|
||||
!--------------------------!
|
||||
! Compute spectral weights !
|
||||
!--------------------------!
|
||||
|
||||
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
|
||||
|
||||
!------------------------!
|
||||
! Compute GW self-energy !
|
||||
!------------------------!
|
||||
|
||||
if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eHF,Om,rho)
|
||||
|
||||
call GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
|
||||
|
||||
!-----------------------------------!
|
||||
! Solve the quasi-particle equation !
|
||||
!-----------------------------------!
|
||||
|
||||
! Linearized or graphical solution?
|
||||
|
||||
eGWlin(:) = eHF(:) + Z(:)*SigC(:)
|
||||
|
||||
if(linearize) then
|
||||
|
||||
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
|
||||
write(*,*)
|
||||
|
||||
eGW(:) = eGWlin(:)
|
||||
|
||||
else
|
||||
|
||||
write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** '
|
||||
write(*,*)
|
||||
|
||||
call GGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eHF,eGW,Z)
|
||||
|
||||
end if
|
||||
|
||||
! call GW_plot_self_energy(nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
|
||||
|
||||
! Compute the RPA correlation energy
|
||||
|
||||
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
|
||||
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||
|
||||
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||
|
||||
!--------------!
|
||||
! Dump results !
|
||||
!--------------!
|
||||
|
||||
call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
|
||||
|
||||
! Deallocate memory
|
||||
|
||||
deallocate(SigC,Z,Om,XpY,XmY,rho)
|
||||
|
||||
! Perform BSE calculation
|
||||
|
||||
! if(dophBSE) then
|
||||
|
||||
! call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
|
||||
|
||||
! if(exchange_kernel) then
|
||||
!
|
||||
! EcBSE(1) = 0.5d0*EcBSE(1)
|
||||
! EcBSE(2) = 1.5d0*EcBSE(2)
|
||||
!
|
||||
! end if
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
|
||||
! Compute the BSE correlation energy via the adiabatic connection
|
||||
|
||||
! if(doACFDT) then
|
||||
|
||||
! write(*,*) '-------------------------------------------------------------'
|
||||
! write(*,*) ' Adiabatic connection version of BSE@G0W0 correlation energy '
|
||||
! write(*,*) '-------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
|
||||
! if(doXBS) then
|
||||
|
||||
! write(*,*) '*** scaled screening version (XBS) ***'
|
||||
! write(*,*)
|
||||
|
||||
! end if
|
||||
|
||||
! call GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcBSE)
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (triplet) =',EcBSE(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
|
||||
! end if
|
||||
|
||||
! end if
|
||||
|
||||
! if(doppBSE) then
|
||||
|
||||
! call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcBSE(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au'
|
||||
! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au'
|
||||
! write(*,*)'-------------------------------------------------------------------------------'
|
||||
! write(*,*)
|
||||
|
||||
! end if
|
||||
|
||||
end subroutine
|
||||
119
src/GW/GGW.f90
Normal file
119
src/GW/GGW.f90
Normal file
@ -0,0 +1,119 @@
|
||||
subroutine GGW(doG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT, &
|
||||
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA, &
|
||||
linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
|
||||
ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF)
|
||||
|
||||
! GW module
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical :: doG0W0
|
||||
logical :: doevGW
|
||||
logical :: doqsGW
|
||||
logical :: doufG0W0
|
||||
logical :: doufGW
|
||||
logical :: doSRGqsGW
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
logical,intent(in) :: doACFDT
|
||||
logical,intent(in) :: exchange_kernel
|
||||
logical,intent(in) :: doXBS
|
||||
logical,intent(in) :: dophBSE
|
||||
logical,intent(in) :: dophBSE2
|
||||
logical,intent(in) :: TDA_W
|
||||
logical,intent(in) :: TDA
|
||||
logical,intent(in) :: dBSE
|
||||
logical,intent(in) :: dTDA
|
||||
logical,intent(in) :: doppBSE
|
||||
logical,intent(in) :: linearize
|
||||
double precision,intent(in) :: eta
|
||||
logical,intent(in) :: regularize
|
||||
|
||||
integer,intent(in) :: nNuc
|
||||
double precision,intent(in) :: ZNuc(nNuc)
|
||||
double precision,intent(in) :: rNuc(nNuc,ncart)
|
||||
double precision,intent(in) :: ENuc
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC(nspin)
|
||||
integer,intent(in) :: nO(nspin)
|
||||
integer,intent(in) :: nV(nspin)
|
||||
integer,intent(in) :: nR(nspin)
|
||||
integer,intent(in) :: nS(nspin)
|
||||
|
||||
double precision,intent(in) :: EHF
|
||||
double precision,intent(in) :: epsHF(nBas)
|
||||
double precision,intent(in) :: cHF(nBas,nBas)
|
||||
double precision,intent(in) :: PHF(nBas,nBas)
|
||||
double precision,intent(in) :: S(nBas,nBas)
|
||||
double precision,intent(in) :: T(nBas,nBas)
|
||||
double precision,intent(in) :: V(nBas,nBas)
|
||||
double precision,intent(in) :: Hc(nBas,nBas)
|
||||
double precision,intent(in) :: X(nBas,nBas)
|
||||
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
|
||||
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: start_GW ,end_GW ,t_GW
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform G0W0 calculatiom
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doG0W0) then
|
||||
|
||||
call wall_time(start_GW)
|
||||
call GG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
|
||||
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
|
||||
call wall_time(end_GW)
|
||||
|
||||
t_GW = end_GW - start_GW
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_GW,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform evGW calculation
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doevGW) then
|
||||
|
||||
call wall_time(start_GW)
|
||||
! call evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
|
||||
! linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
|
||||
call wall_time(end_GW)
|
||||
|
||||
t_GW = end_GW - start_GW
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGW = ',t_GW,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Perform qsGW calculation
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doqsGW) then
|
||||
|
||||
call wall_time(start_GW)
|
||||
! call qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
|
||||
! eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, &
|
||||
! dipole_int_AO,dipole_int,PHF,cHF,epsHF)
|
||||
call wall_time(end_GW)
|
||||
|
||||
t_GW = end_GW - start_GW
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_GW,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
82
src/GW/GGW_QP_graph.f90
Normal file
82
src/GW/GGW_QP_graph.f90
Normal file
@ -0,0 +1,82 @@
|
||||
subroutine GGW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eOld,eGW,Z)
|
||||
|
||||
! Compute the graphical solution of the QP equation
|
||||
|
||||
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) :: eHF(nBas)
|
||||
double precision,intent(in) :: Om(nS)
|
||||
double precision,intent(in) :: rho(nBas,nBas,nS)
|
||||
|
||||
double precision,intent(in) :: eGWlin(nBas)
|
||||
double precision,intent(in) :: eOld(nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: p
|
||||
integer :: nIt
|
||||
integer,parameter :: maxIt = 64
|
||||
double precision,parameter :: thresh = 1d-6
|
||||
double precision,external :: GGW_SigC,GGW_dSigC
|
||||
double precision :: SigC,dSigC
|
||||
double precision :: f,df
|
||||
double precision :: w
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: eGW(nBas)
|
||||
double precision,intent(out) :: Z(nBas)
|
||||
|
||||
! Run Newton's algorithm to find the root
|
||||
|
||||
write(*,*)'-----------------------------------------------------'
|
||||
write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','e_GWlin (eV)','e_GW (eV)','Z'
|
||||
write(*,*)'-----------------------------------------------------'
|
||||
|
||||
do p=nC+1,nBas-nR
|
||||
|
||||
w = eGWlin(p)
|
||||
nIt = 0
|
||||
f = 1d0
|
||||
|
||||
do while (abs(f) > thresh .and. nIt < maxIt)
|
||||
|
||||
nIt = nIt + 1
|
||||
|
||||
SigC = GGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho)
|
||||
dSigC = GGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,eOld,Om,rho)
|
||||
f = w - eHF(p) - SigC
|
||||
df = 1d0/(1d0 - dSigC)
|
||||
|
||||
w = w - df*f
|
||||
|
||||
end do
|
||||
|
||||
if(nIt == maxIt) then
|
||||
|
||||
eGW(p) = eGWlin(p)
|
||||
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,eGWlin(p)*HaToeV,eGW(p)*HaToeV,Z(p),'Cvg Failed!'
|
||||
|
||||
else
|
||||
|
||||
eGW(p) = w
|
||||
Z(p) = df
|
||||
|
||||
write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,eGWlin(p)*HaToeV,eGW(p)*HaToeV,Z(p)
|
||||
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
end subroutine
|
||||
52
src/GW/GGW_SigC.f90
Normal file
52
src/GW/GGW_SigC.f90
Normal file
@ -0,0 +1,52 @@
|
||||
double precision function GGW_SigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
|
||||
|
||||
! Compute diagonal of the correlation part of the self-energy
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: p
|
||||
double precision,intent(in) :: w
|
||||
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) :: e(nBas)
|
||||
double precision,intent(in) :: Om(nS)
|
||||
double precision,intent(in) :: rho(nBas,nBas,nS)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,a,m
|
||||
double precision :: num,eps
|
||||
|
||||
! Initialize
|
||||
|
||||
GGW_SigC = 0d0
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do i=nC+1,nO
|
||||
do m=1,nS
|
||||
eps = w - e(i) + Om(m)
|
||||
num = rho(p,i,m)**2
|
||||
GGW_SigC = GGW_SigC + num*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do a=nO+1,nBas-nR
|
||||
do m=1,nS
|
||||
eps = w - e(a) - Om(m)
|
||||
num = rho(p,a,m)**2
|
||||
GGW_SigC = GGW_SigC + num*eps/(eps**2 + eta**2)
|
||||
end do
|
||||
end do
|
||||
|
||||
end function
|
||||
52
src/GW/GGW_dSigC.f90
Normal file
52
src/GW/GGW_dSigC.f90
Normal file
@ -0,0 +1,52 @@
|
||||
double precision function GGW_dSigC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Om,rho)
|
||||
|
||||
! Compute the derivative of the correlation part of the self-energy
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: p
|
||||
double precision,intent(in) :: w
|
||||
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) :: e(nBas)
|
||||
double precision,intent(in) :: Om(nS)
|
||||
double precision,intent(in) :: rho(nBas,nBas,nS)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,a,m
|
||||
double precision :: num,eps
|
||||
|
||||
! Initialize
|
||||
|
||||
GGW_dSigC = 0d0
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do i=nC+1,nO
|
||||
do m=1,nS
|
||||
eps = w - e(i) + Om(m)
|
||||
num = rho(p,i,m)**2
|
||||
GGW_dSigC = GGW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do a=nO+1,nBas-nR
|
||||
do m=1,nS
|
||||
eps = w - e(a) - Om(m)
|
||||
num = rho(p,a,m)**2
|
||||
GGW_dSigC = GGW_dSigC - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
end do
|
||||
end do
|
||||
|
||||
end function
|
||||
90
src/GW/GGW_self_energy_diag.f90
Normal file
90
src/GW/GGW_self_energy_diag.f90
Normal file
@ -0,0 +1,90 @@
|
||||
subroutine GGW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
|
||||
|
||||
! Compute diagonal of the correlation part of the self-energy and the renormalization factor
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
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) :: e(nBas)
|
||||
double precision,intent(in) :: Om(nS)
|
||||
double precision,intent(in) :: rho(nBas,nBas,nS)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,a,p,m
|
||||
double precision :: num,eps
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: Sig(nBas)
|
||||
double precision,intent(out) :: Z(nBas)
|
||||
double precision,intent(out) :: EcGM
|
||||
|
||||
! Initialize
|
||||
|
||||
Sig(:) = 0d0
|
||||
Z(:) = 0d0
|
||||
|
||||
!----------------!
|
||||
! GW self-energy !
|
||||
!----------------!
|
||||
|
||||
! Occupied part of the correlation self-energy
|
||||
|
||||
do p=nC+1,nBas-nR
|
||||
do i=nC+1,nO
|
||||
do m=1,nS
|
||||
|
||||
eps = e(p) - e(i) + Om(m)
|
||||
num = rho(p,i,m)**2
|
||||
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
|
||||
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Virtual part of the correlation self-energy
|
||||
|
||||
do p=nC+1,nBas-nR
|
||||
do a=nO+1,nBas-nR
|
||||
do m=1,nS
|
||||
|
||||
eps = e(p) - e(a) - Om(m)
|
||||
num = rho(p,a,m)**2
|
||||
Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2)
|
||||
Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Galitskii-Migdal correlation energy
|
||||
|
||||
EcGM = 0d0
|
||||
do i=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
do m=1,nS
|
||||
|
||||
eps = e(a) - e(i) + Om(m)
|
||||
num = rho(a,i,m)**2
|
||||
EcGM = EcGM - num*eps/(eps**2 + eta**2)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
! Compute renormalization factor from derivative
|
||||
|
||||
Z(:) = 1d0/(1d0 - Z(:))
|
||||
|
||||
end subroutine
|
||||
@ -56,6 +56,7 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
|
||||
double precision,allocatable :: F(:,:)
|
||||
double precision,allocatable :: Fp(:,:)
|
||||
double precision,allocatable :: Cp(:,:)
|
||||
double precision,allocatable :: H(:,:)
|
||||
double precision,allocatable :: S(:,:)
|
||||
double precision,allocatable :: X(:,:)
|
||||
double precision,allocatable :: err(:,:)
|
||||
@ -96,16 +97,9 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
|
||||
Faa(nBas,nBas),Fab(nBas,nBas),Fba(nBas,nBas),Fbb(nBas,nBas), &
|
||||
Paa(nBas,nBas),Pab(nBas,nBas),Pba(nBas,nBas),Pbb(nBas,nBas), &
|
||||
F(nBas2,nBas2),Fp(nBas2,nBas2),Cp(nBas2,nBas2), &
|
||||
S(nBas2,nBas2),X(nBas2,nBas2),err(nBas2,nBas2), &
|
||||
H(nBas2,nBas2),S(nBas2,nBas2),X(nBas2,nBas2),err(nBas2,nBas2), &
|
||||
err_diis(nBas2Sq,max_diis),F_diis(nBas2Sq,max_diis))
|
||||
|
||||
! Guess coefficients and demsity matrices
|
||||
|
||||
! do ispin=1,nspin
|
||||
! call mo_guess(nBas,guess_type,S,Hc,X,c(:,:,ispin))
|
||||
! P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
|
||||
! end do
|
||||
|
||||
! Initialization
|
||||
|
||||
nSCF = 0
|
||||
@ -127,6 +121,23 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
|
||||
X( 1:nBas , 1:nBas ) = Or(1:nBas,1:nBas)
|
||||
X(nBas+1:nBas2,nBas+1:nBas2) = Or(1:nBas,1:nBas)
|
||||
|
||||
! Construct super orthogonalization matrix
|
||||
|
||||
H( : , : ) = 0d0
|
||||
H( 1:nBas , 1:nBas ) = Hc(1:nBas,1:nBas)
|
||||
H(nBas+1:nBas2,nBas+1:nBas2) = Hc(1:nBas,1:nBas)
|
||||
|
||||
! Guess coefficients and density matrices
|
||||
|
||||
call mo_guess(nBas2,guess_type,S,H,X,C)
|
||||
|
||||
P(:,:) = matmul(C(:,1:nOcc),transpose(C(:,1:nOcc)))
|
||||
|
||||
Paa(:,:) = P(1:nBas,1:nBas)
|
||||
Pab(:,:) = P(1:nBas,nBas+1:nBas2)
|
||||
Pba(:,:) = P(nBas+1:nBas2,1:nBas)
|
||||
Pbb(:,:) = P(nBas+1:nBAs2,nBAs+1:nBas2)
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
@ -146,8 +157,6 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
|
||||
! Build individual Hartree matrices
|
||||
|
||||
call Hartree_matrix_AO_basis(nBas,Paa,ERI,Jaa)
|
||||
! call Hartree_matrix_AO_basis(nBas,Pab,ERI,Jab)
|
||||
! call Hartree_matrix_AO_basis(nBas,Pba,ERI,Jba)
|
||||
call Hartree_matrix_AO_basis(nBas,Pbb,ERI,Jbb)
|
||||
|
||||
! Compute individual exchange matrices
|
||||
@ -167,8 +176,8 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
|
||||
! Build super Fock matrix
|
||||
|
||||
F( 1:nBas , 1:nBas ) = Faa(1:nBas,1:nBas)
|
||||
F(nBas+1:nBas2, 1:nBas ) = Fab(1:nBas,1:nBas)
|
||||
F( 1:nBas ,nBas+1:nBas2) = Fba(1:nBas,1:nBas)
|
||||
F( 1:nBas ,nBas+1:nBas2) = Fab(1:nBas,1:nBas)
|
||||
F(nBas+1:nBas2, 1:nBas ) = Fba(1:nBas,1:nBas)
|
||||
F(nBas+1:nBas2,nBas+1:nBas2) = Fbb(1:nBas,1:nBas)
|
||||
|
||||
! Check convergence
|
||||
@ -217,11 +226,12 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
|
||||
! if(nSCF == 1) call mix_guess(nBas,nO,mix,c)
|
||||
|
||||
! Compute individual density matrices
|
||||
P(:,:) = matmul(C(:,1:nOcc),transpose(C(:,1:nOcc)))
|
||||
|
||||
Paa(:,:) = matmul(Ca(:,1:nOcc),transpose(Ca(:,1:nOcc)))
|
||||
Pab(:,:) = matmul(Ca(:,1:nOcc),transpose(Cb(:,1:nOcc)))
|
||||
Pba(:,:) = matmul(Cb(:,1:nOcc),transpose(Ca(:,1:nOcc)))
|
||||
Pbb(:,:) = matmul(Cb(:,1:nOcc),transpose(Cb(:,1:nOcc)))
|
||||
Paa(:,:) = P(1:nBas,1:nBas)
|
||||
Pab(:,:) = P(1:nBas,nBas+1:nBas2)
|
||||
Pba(:,:) = P(nBas+1:nBas2,1:nBas)
|
||||
Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2)
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute UHF energy
|
||||
@ -243,19 +253,19 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
|
||||
|
||||
! Hartree energy: 16 terms?
|
||||
|
||||
EJaaaa = +0.5d0*trace_matrix(nBas,matmul(Paa,Jaa))
|
||||
EJaabb = +0.5d0*trace_matrix(nBas,matmul(Paa,Jbb))
|
||||
EJbbaa = +0.5d0*trace_matrix(nBas,matmul(Pbb,Jaa))
|
||||
EJbbbb = +0.5d0*trace_matrix(nBas,matmul(Pbb,Jbb))
|
||||
EJaaaa = 0.5d0*trace_matrix(nBas,matmul(Paa,Jaa))
|
||||
EJaabb = 0.5d0*trace_matrix(nBas,matmul(Paa,Jbb))
|
||||
EJbbaa = 0.5d0*trace_matrix(nBas,matmul(Pbb,Jaa))
|
||||
EJbbbb = 0.5d0*trace_matrix(nBas,matmul(Pbb,Jbb))
|
||||
|
||||
EJ = EJaaaa + EJaabb + EJbbaa + EJbbbb
|
||||
|
||||
! Exchange energy
|
||||
|
||||
Exaaaa = -0.5d0*trace_matrix(nBas,matmul(Paa,Kaa))
|
||||
Exabba = -0.5d0*trace_matrix(nBas,matmul(Pab,Kba))
|
||||
Exbaab = -0.5d0*trace_matrix(nBas,matmul(Pba,Kab))
|
||||
Exbbbb = -0.5d0*trace_matrix(nBas,matmul(Pbb,Kbb))
|
||||
Exaaaa = 0.5d0*trace_matrix(nBas,matmul(Paa,Kaa))
|
||||
Exabba = 0.5d0*trace_matrix(nBas,matmul(Pab,Kba))
|
||||
Exbaab = 0.5d0*trace_matrix(nBas,matmul(Pba,Kab))
|
||||
Exbbbb = 0.5d0*trace_matrix(nBas,matmul(Pbb,Kbb))
|
||||
|
||||
Ex = Exaaaa + Exabba + Exbaab + Exbbbb
|
||||
|
||||
|
||||
109
src/HF/GHF_stability.f90
Normal file
109
src/HF/GHF_stability.f90
Normal file
@ -0,0 +1,109 @@
|
||||
subroutine GHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
|
||||
|
||||
! Perform a stability analysis of the GHF solution
|
||||
|
||||
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) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer,parameter :: maxS = 20
|
||||
integer :: ia
|
||||
integer :: ispin
|
||||
|
||||
double precision,allocatable :: A(:,:)
|
||||
double precision,allocatable :: B(:,:)
|
||||
double precision,allocatable :: AB(:,:)
|
||||
double precision,allocatable :: Om(:)
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(A(nS,nS),B(nS,nS),AB(nS,nS),Om(nS))
|
||||
|
||||
!-------------------------------------------------------------!
|
||||
! Stability analysis: Real GHF -> Real GHF
|
||||
!-------------------------------------------------------------!
|
||||
|
||||
ispin = 3
|
||||
|
||||
call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
|
||||
call phLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
|
||||
|
||||
AB(:,:) = A(:,:) + B(:,:)
|
||||
|
||||
call diagonalize_matrix(nS,AB,Om)
|
||||
|
||||
write(*,*)'-------------------------------------------------------------'
|
||||
write(*,*)'| Stability analysis: Real GHF -> Real GHF |'
|
||||
write(*,*)'-------------------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
|
||||
'|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'
|
||||
write(*,*)'-------------------------------------------------------------'
|
||||
do ia=1,min(nS,maxS)
|
||||
write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') &
|
||||
'|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|'
|
||||
enddo
|
||||
write(*,*)'-------------------------------------------------------------'
|
||||
|
||||
if(minval(Om(:)) < 0d0) then
|
||||
|
||||
write(*,'(1X,A40,1X)') 'Too bad, GHF solution is unstable!'
|
||||
write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om(1),' au'
|
||||
|
||||
else
|
||||
|
||||
write(*,'(1X,A40,1X)') 'Well done, GHF solution is stable!'
|
||||
write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om(1),' au'
|
||||
|
||||
end if
|
||||
write(*,*)'-------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
!-------------------------------------------------------------!
|
||||
! Stability analysis: Real GHF -> Complex GHF
|
||||
!-------------------------------------------------------------!
|
||||
|
||||
AB(:,:) = A(:,:) - B(:,:)
|
||||
|
||||
call diagonalize_matrix(nS,AB,Om)
|
||||
|
||||
write(*,*)'-------------------------------------------------------------'
|
||||
write(*,*)'| Stability analysis: Real GHF -> Complex GHF |'
|
||||
write(*,*)'-------------------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
|
||||
'|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'
|
||||
write(*,*)'-------------------------------------------------------------'
|
||||
do ia=1,min(nS,maxS)
|
||||
write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') &
|
||||
'|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|'
|
||||
enddo
|
||||
write(*,*)'-------------------------------------------------------------'
|
||||
|
||||
if(minval(Om(:)) < 0d0) then
|
||||
|
||||
write(*,'(1X,A40,1X)') 'Too bad, GHF solution is unstable!'
|
||||
write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om(1),' au'
|
||||
|
||||
else
|
||||
|
||||
write(*,'(1X,A40,1X)') 'Well done, GHF solution is stable!'
|
||||
write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om(1),' au'
|
||||
|
||||
end if
|
||||
write(*,*)'-------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
|
||||
end subroutine
|
||||
@ -28,6 +28,8 @@ subroutine mo_guess(nBas,guess_type,S,Hc,X,c)
|
||||
|
||||
call random_number(c)
|
||||
|
||||
c(:,:) = 2d0*c(:,:) - 1d0
|
||||
|
||||
else
|
||||
|
||||
print*,'Wrong guess option'
|
||||
|
||||
45
src/MP/GMP.f90
Normal file
45
src/MP/GMP.f90
Normal file
@ -0,0 +1,45 @@
|
||||
subroutine GMP(doMP2,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
|
||||
|
||||
! Moller-Plesset module
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: doMP2
|
||||
|
||||
logical,intent(in) :: regularize
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: EHF
|
||||
double precision,intent(in) :: epsHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: start_MP ,end_MP ,t_MP
|
||||
|
||||
! Output variables
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute MP2 energy
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doMP2) then
|
||||
|
||||
call wall_time(start_MP)
|
||||
call GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
|
||||
call wall_time(end_MP)
|
||||
|
||||
t_MP = end_MP - start_MP
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
174
src/MP/GMP2.f90
Normal file
174
src/MP/GMP2.f90
Normal file
@ -0,0 +1,174 @@
|
||||
subroutine GMP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e)
|
||||
|
||||
! Perform second-order Moller-Plesset calculation with and without regularizers
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: regularize
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC
|
||||
integer,intent(in) :: nO
|
||||
integer,intent(in) :: nV
|
||||
integer,intent(in) :: nR
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: EHF
|
||||
double precision,intent(in) :: e(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,j,a,b
|
||||
|
||||
double precision :: kappa,sigm1,sigm2
|
||||
double precision :: Dijab
|
||||
double precision :: fs,fs2,fk
|
||||
|
||||
double precision :: E2d,E2ds,E2ds2,E2dk
|
||||
double precision :: E2x,E2xs,E2xs2,E2xk
|
||||
|
||||
double precision :: EcsMP2,Ecs2MP2,EckMP2
|
||||
double precision :: EcMP2
|
||||
|
||||
! Output variables
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'************************************************'
|
||||
write(*,*)'| Moller-Plesset second-order calculation |'
|
||||
write(*,*)'************************************************'
|
||||
write(*,*)
|
||||
|
||||
!---------------------------------------------!
|
||||
! Parameters for regularized MP2 calculations !
|
||||
!---------------------------------------------!
|
||||
|
||||
kappa = 1.1d0
|
||||
sigm1 = 0.7d0
|
||||
sigm2 = 0.4d0
|
||||
|
||||
!--------------------------------------------------!
|
||||
! Compute conventinal and regularized MP2 energies !
|
||||
!--------------------------------------------------!
|
||||
|
||||
E2d = 0d0
|
||||
E2ds = 0d0
|
||||
E2ds2 = 0d0
|
||||
E2dk = 0d0
|
||||
|
||||
E2x = 0d0
|
||||
E2xs = 0d0
|
||||
E2xs2 = 0d0
|
||||
E2xk = 0d0
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=nO+1,nBas-nR
|
||||
do b=nO+1,nBas-nR
|
||||
|
||||
Dijab = e(a) + e(b) - e(i) - e(j)
|
||||
|
||||
! Second-order ring diagram
|
||||
|
||||
fs = (1d0 - exp(-sigm1*Dijab))/Dijab
|
||||
fs2 = (1d0 - exp(-sigm2*Dijab*Dijab))/Dijab
|
||||
fk = (1d0 - exp(-kappa*Dijab))**2/Dijab
|
||||
|
||||
E2d = E2d - 0.25d0*(ERI(i,j,a,b)-ERI(i,j,b,a))**2/Dijab
|
||||
E2ds = E2ds - ERI(i,j,a,b)*ERI(i,j,a,b)*fs
|
||||
E2ds2 = E2ds2 - ERI(i,j,a,b)*ERI(i,j,a,b)*fs2
|
||||
E2dk = E2dk - ERI(i,j,a,b)*ERI(i,j,a,b)*fk
|
||||
|
||||
! Second-order exchange diagram
|
||||
|
||||
E2x = E2x - ERI(i,j,a,b)*ERI(i,j,b,a)/Dijab
|
||||
E2x = 0d0
|
||||
E2xs = E2xs - ERI(i,j,a,b)*ERI(i,j,b,a)*fs
|
||||
E2xs2 = E2xs2 - ERI(i,j,a,b)*ERI(i,j,b,a)*fs2
|
||||
E2xk = E2xk - ERI(i,j,a,b)*ERI(i,j,b,a)*fk
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
EcMP2 = E2d - E2x
|
||||
EcsMP2 = E2ds - E2xs
|
||||
Ecs2MP2 = E2ds2 - E2xs2
|
||||
EckMP2 = E2dk - E2xk
|
||||
|
||||
!------------!
|
||||
! MP2 energy !
|
||||
!------------!
|
||||
|
||||
write(*,*)
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,'(A32)') ' GMP2 calculation '
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,'(A32,1X,F16.10)') ' GMP2 correlation energy = ',EcMP2
|
||||
write(*,'(A32,1X,F16.10)') ' Direct part = ',E2d
|
||||
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2x
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EHF + EcMP2
|
||||
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EHF + EcMP2
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,*)
|
||||
|
||||
if(regularize) then
|
||||
|
||||
!-------------------!
|
||||
! sigma1-MP2 energy !
|
||||
!-------------------!
|
||||
|
||||
write(*,*)
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,'(A32)') ' sigma-GMP2 calculation '
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,'(A32,1X,F16.10)') ' GMP2 correlation energy = ',EcsMP2
|
||||
write(*,'(A32,1X,F16.10)') ' Direct part = ',E2ds
|
||||
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xs
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EHF + EcsMP2
|
||||
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EHF + EcsMP2
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,*)
|
||||
|
||||
!--------------------!
|
||||
! sigma^2-MP2 energy !
|
||||
!--------------------!
|
||||
|
||||
write(*,*)
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,'(A32)') ' sigma^2-GMP2 calculation '
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,'(A32,1X,F16.10)') ' GMP2 correlation energy = ',Ecs2MP2
|
||||
write(*,'(A32,1X,F16.10)') ' Direct part = ',E2ds2
|
||||
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xs2
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EHF + Ecs2MP2
|
||||
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EHF + Ecs2MP2
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,*)
|
||||
|
||||
!------------------!
|
||||
! kappa-MP2 energy !
|
||||
!------------------!
|
||||
|
||||
write(*,*)
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,'(A32)') ' kappa-GMP2 calculation '
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,'(A32,1X,F16.10)') ' GMP2 correlation energy = ',EckMP2
|
||||
write(*,'(A32,1X,F16.10)') ' Direct part = ',E2dk
|
||||
write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xk
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,'(A32,1X,F16.10)') ' GMP2 electronic energy = ',EHF + EckMP2
|
||||
write(*,'(A32,1X,F16.10)') ' GMP2 total energy = ',ENuc + EHF + EckMP2
|
||||
write(*,'(A32)') '---------------------------'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
@ -72,6 +72,8 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
|
||||
double precision :: EHF
|
||||
double precision,allocatable :: dipole_int_MO(:,:,:)
|
||||
double precision,allocatable :: ERI_MO(:,:,:,:)
|
||||
double precision,allocatable :: ERI_tmp(:,:,:,:)
|
||||
double precision,allocatable :: Ca(:,:),Cb(:,:)
|
||||
|
||||
integer :: ixyz
|
||||
integer :: nBas2
|
||||
@ -128,7 +130,17 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
|
||||
|
||||
! 4-index transform
|
||||
|
||||
call AOtoMO_integral_transform(1,1,1,1,nBas2,cHF,ERI_AO,ERI_MO)
|
||||
allocate(Ca(nBas,nBas2),Cb(nBas,nBas2),ERI_tmp(nBas2,nBas2,nBas2,nBas2))
|
||||
Ca(:,:) = cHF(1:nBas,1:nBas2)
|
||||
Cb(:,:) = cHF(nBas+1:nBas2,1:nBas2)
|
||||
call AOtoMO_integral_transform_GHF(nBas,nBas2,Ca,Ca,Ca,Ca,ERI_AO,ERI_tmp)
|
||||
ERI_MO(:,:,:,:) = ERI_tmp(:,:,:,:)
|
||||
call AOtoMO_integral_transform_GHF(nBas,nBas2,Ca,Cb,Ca,Cb,ERI_AO,ERI_tmp)
|
||||
ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:)
|
||||
call AOtoMO_integral_transform_GHF(nBas,nBas2,Cb,Ca,Cb,Ca,ERI_AO,ERI_tmp)
|
||||
ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:)
|
||||
call AOtoMO_integral_transform_GHF(nBas,nBas2,Cb,Cb,Cb,Cb,ERI_AO,ERI_tmp)
|
||||
ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:)
|
||||
|
||||
call wall_time(end_AOtoMO)
|
||||
|
||||
@ -143,7 +155,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
|
||||
if(dostab) then
|
||||
|
||||
call wall_time(start_stab)
|
||||
! call GHF_stability(nBas,nC,nO,nV,nR,nS,epsHF,ERI_MO)
|
||||
call GHF_stability(nBas2,sum(nC),sum(nO),sum(nV),sum(nR),sum(nO)*sum(nV),epsHF,ERI_MO)
|
||||
call wall_time(end_stab)
|
||||
|
||||
t_stab = end_stab - start_stab
|
||||
@ -161,7 +173,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
|
||||
if(doMP) then
|
||||
|
||||
call wall_time(start_MP)
|
||||
! call GMP(doMP2,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
|
||||
call GMP(doMP2,reg_MP,nBas2,sum(nC),sum(nO),sum(nV),sum(nR),ERI_MO,ENuc,EHF,epsHF)
|
||||
call wall_time(end_MP)
|
||||
|
||||
t_MP = end_MP - start_MP
|
||||
@ -179,8 +191,8 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
|
||||
if(doRPA) then
|
||||
|
||||
call wall_time(start_RPA)
|
||||
! call GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip, &
|
||||
! nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S)
|
||||
call GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel, &
|
||||
nBas2,sum(nC),sum(nO),sum(nV),sum(nR),sum(nO)*sum(nV),ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S)
|
||||
call wall_time(end_RPA)
|
||||
|
||||
t_RPA = end_RPA - start_RPA
|
||||
@ -219,10 +231,10 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
|
||||
if(doGW) then
|
||||
|
||||
call wall_time(start_GW)
|
||||
! call GGW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, &
|
||||
! exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA, &
|
||||
! lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
|
||||
! ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
|
||||
call GGW(doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, &
|
||||
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA, &
|
||||
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas2,sum(nC),sum(nO),sum(nV),sum(nR),sum(nO)*sum(nV), &
|
||||
EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
|
||||
call wall_time(end_GW)
|
||||
|
||||
t_GW = end_GW - start_GW
|
||||
|
||||
82
src/RPA/GRPA.f90
Normal file
82
src/RPA/GRPA.f90
Normal file
@ -0,0 +1,82 @@
|
||||
subroutine GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
|
||||
|
||||
! Random-phase approximation module
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical :: dophRPA
|
||||
logical :: dophRPAx
|
||||
logical :: doppRPA
|
||||
|
||||
logical,intent(in) :: TDA
|
||||
logical,intent(in) :: doACFDT
|
||||
logical,intent(in) :: exchange_kernel
|
||||
integer,intent(in) :: nBas
|
||||
integer,intent(in) :: nC(nspin)
|
||||
integer,intent(in) :: nO(nspin)
|
||||
integer,intent(in) :: nV(nspin)
|
||||
integer,intent(in) :: nR(nspin)
|
||||
integer,intent(in) :: nS(nspin)
|
||||
double precision,intent(in) :: ENuc
|
||||
double precision,intent(in) :: EHF
|
||||
double precision,intent(in) :: epsHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: start_RPA ,end_RPA ,t_RPA
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute (direct) RPA excitations
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(dophRPA) then
|
||||
|
||||
call wall_time(start_RPA)
|
||||
call phGRPA(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
|
||||
call wall_time(end_RPA)
|
||||
|
||||
t_RPA = end_RPA - start_RPA
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute RPAx (RPA with exchange) excitations
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(dophRPAx) then
|
||||
|
||||
call wall_time(start_RPA)
|
||||
call phGRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
|
||||
|
||||
call wall_time(end_RPA)
|
||||
|
||||
t_RPA = end_RPA - start_RPA
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPA,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute ppRPA excitations
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doppRPA) then
|
||||
|
||||
call wall_time(start_RPA)
|
||||
! call ppGRPA(TDA,doACFDT,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,epsHF)
|
||||
call wall_time(end_RPA)
|
||||
|
||||
t_RPA = end_RPA - start_RPA
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
end subroutine
|
||||
78
src/RPA/phGRPA.f90
Normal file
78
src/RPA/phGRPA.f90
Normal file
@ -0,0 +1,78 @@
|
||||
subroutine phGRPA(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
|
||||
|
||||
! Perform a direct random phase approximation calculation
|
||||
|
||||
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) :: EHF
|
||||
double precision,intent(in) :: e(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: ispin
|
||||
logical :: dRPA
|
||||
double precision,allocatable :: Aph(:,:)
|
||||
double precision,allocatable :: Bph(:,:)
|
||||
double precision,allocatable :: Om(:)
|
||||
double precision,allocatable :: XpY(:,:)
|
||||
double precision,allocatable :: XmY(:,:)
|
||||
|
||||
double precision :: EcRPA
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'***********************************************'
|
||||
write(*,*)'| Random-phase approximation calculation |'
|
||||
write(*,*)'***********************************************'
|
||||
write(*,*)
|
||||
|
||||
! TDA
|
||||
|
||||
if(TDA) then
|
||||
write(*,*) 'Tamm-Dancoff approximation activated!'
|
||||
write(*,*)
|
||||
end if
|
||||
|
||||
! Initialization
|
||||
|
||||
dRPA = .true.
|
||||
|
||||
EcRPA = 0d0
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS))
|
||||
if(.not.TDA) allocate(Bph(nS,nS))
|
||||
|
||||
ispin = 3
|
||||
|
||||
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
|
||||
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||
|
||||
call phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||
call print_excitation_energies('phRPA@GHF',ispin,nS,Om)
|
||||
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@phRPA correlation energy =',EcRPA
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@phRPA total energy =',ENuc + EHF + EcRPA
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
end subroutine
|
||||
79
src/RPA/phGRPAx.f90
Normal file
79
src/RPA/phGRPAx.f90
Normal file
@ -0,0 +1,79 @@
|
||||
subroutine phGRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,e)
|
||||
|
||||
! Perform random phase approximation calculation with exchange (aka TDHF)
|
||||
|
||||
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) :: EHF
|
||||
double precision,intent(in) :: e(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: ispin
|
||||
logical :: dRPA
|
||||
double precision,allocatable :: Aph(:,:)
|
||||
double precision,allocatable :: Bph(:,:)
|
||||
double precision,allocatable :: Om(:)
|
||||
double precision,allocatable :: XpY(:,:)
|
||||
double precision,allocatable :: XmY(:,:)
|
||||
|
||||
double precision :: EcRPA
|
||||
|
||||
! 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
|
||||
|
||||
dRPA = .false.
|
||||
|
||||
EcRPA = 0d0
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS))
|
||||
if(.not.TDA) allocate(Bph(nS,nS))
|
||||
|
||||
ispin = 3
|
||||
|
||||
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,Aph)
|
||||
if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
|
||||
|
||||
call phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
|
||||
call print_excitation_energies('phRPAx@GHF',ispin,nS,Om)
|
||||
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@phRPAx correlation energy =',EcRPA
|
||||
write(*,'(2X,A50,F20.10)') 'Tr@phRPAx total energy =',ENuc + EHF + EcRPA
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
end subroutine
|
||||
Loading…
Reference in New Issue
Block a user