4
1
mirror of https://github.com/pfloos/quack synced 2024-06-26 15:12:17 +02:00
quack/src/QuAcK/TDHF.f90
2020-01-14 16:36:11 +01:00

102 lines
3.1 KiB
Fortran

subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e)
! Perform random phase approximation calculation
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
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) :: e(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 :: rho
double precision :: EcRPA(nspin)
logical :: adiabatic_connection
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Time-dependent Hartree-Fock calculation |'
write(*,*)'************************************************'
write(*,*)
! Initialization
EcRPA(:) = 0d0
! Memory allocation
allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
! Singlet manifold
if(singlet_manifold) then
ispin = 1
call linear_response(ispin,.false.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
endif
! Triplet manifold
if(triplet_manifold) then
ispin = 2
call linear_response(ispin,.false.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
endif
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'Tr@TDHF correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A40,F15.6)') 'Tr@TDHF correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@TDHF correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@TDHF total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Compute the correlation energy via the adiabatic connection
adiabatic_connection = .true.
if(adiabatic_connection) then
write(*,*) '-------------------------------------------------------'
write(*,*) 'Adiabatic connection version of TDHF correlation energy'
write(*,*) '-------------------------------------------------------'
write(*,*)
call ACDFT(.false.,.false.,.false.,.false.,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho)
end if
end subroutine TDHF