subroutine tdhf_gradient(infos)
use io_constants, only: iw
use oqp_tagarray_driver
use types, only: information
use strings, only: Cstring, fstring
use basis_tools, only: basis_set
use messages, only: show_message, with_abort
use grd1, only: eijden, print_gradient
use util, only: measure_time
use tdhf_lib, only: &
iatogen, mntoia
use dft, only: dft_initialize, dftclean
use mathlib, only: symmetrize_matrix, orthogonal_transform
use mod_dft_molgrid, only: dft_grid_t
use mod_dft_gridint_tdxc_grad, only: tddft_xc_gradient
use mathlib, only: unpack_matrix
use printing, only: print_module_info
implicit none
character(len=*), parameter :: subroutine_name = "tdhf_gradient"
type(basis_set), pointer :: basis
type(information), target, intent(inout) :: infos
integer :: nbf, nocc
type(dft_grid_t) :: molGrid
! General data
logical :: dft
integer :: scf_type, mol_mult
real(kind=dp), allocatable :: p(:,:,:), d(:,:,:), xpy2(:,:,:), xmy2(:,:,:), wrk1(:,:), wrk2(:,:)
! tagarray
real(kind=dp), contiguous, pointer :: &
dmat_a(:), mo_a(:,:), td_p(:,:), xpy(:), xmy(:)
character(len=*), parameter :: tags_general(*) = [ character(len=80) :: &
OQP_TD_XPY, OQP_TD_XMY, OQP_DM_A, OQP_VEC_MO_A, OQP_TD_P ]
mol_mult = infos%mol_prop%mult
if (mol_mult/=1) call show_message(&
'RPA-TDDFT are only available for RHF reference', with_abort)
scf_type = infos%control%scftype
if (scf_type/=1) error stop
dft = infos%control%hamilton == 20
! Files open
open (unit=IW, file=infos%log_filename, position="append")
!
call print_module_info('TDHF_Grad','Computing Grdient of TDDFT')
!
write(iw,'(/5X,"Gradient options"/&
&5X,18("-")/&
&5X,"Target State: ",I8/)') infos%tddft%target_state
! Load basis set
basis => infos%basis
basis%atoms => infos%atoms
call data_has_tags(infos%dat, tags_general, module_name, subroutine_name, WITH_ABORT)
call tagarray_get_data(infos%dat, OQP_DM_A, dmat_a)
call tagarray_get_data(infos%dat, OQP_VEC_MO_A, mo_a)
call tagarray_get_data(infos%dat, OQP_td_p, td_p)
call tagarray_get_data(infos%dat, OQP_td_xpy, xpy)
call tagarray_get_data(infos%dat, OQP_td_xmy, xmy)
! Allocate H, S ,T and D matrices
nbf = basis%nbf
nocc = infos%mol_prop%nocc
! Compute 1e gradient
call flush(iw)
call tdhf_1e_grad(infos, basis)
!call print_gradient(infos)
write(iw,"(' ..... End Of 1-Eelectron Gradient ......')")
call measure_time(print_total=1, log_unit=iw)
call flush(iw)
allocate(d(nbf,nbf,2), &
p(nbf,nbf,2), &
xpy2(nbf,nbf,2), &
xmy2(nbf,nbf,2), &
wrk1(nbf,nbf), &
wrk2(nbf,nbf), &
source=0.0d0)
call iatogen(xpy, wrk1, nocc, nocc)
call symmetrize_matrix(wrk1, nbf)
wrk1 = wrk1*0.5
call orthogonal_transform('t', nbf, mo_a, wrk1, xpy2(:,:,1), wrk2)
call iatogen(xmy, wrk1, nocc, nocc)
call orthogonal_transform('t', nbf, mo_a, wrk1, xmy2(:,:,1), wrk2)
call unpack_matrix(td_p(:,1), p(:,:,1))
call unpack_matrix(dmat_a, d(:,:,1))
! Compute xc gradient
if (dft) then
call dft_initialize(infos, basis, molGrid)
d(:,:,2) = d(:,:,1)
p(:,:,2) = p(:,:,1)
xpy2(:,:,2) = xpy2(:,:,1)
call tddft_xc_gradient(basis=basis, &
molGrid=molGrid, &
dedft=infos%atoms%grad, &
da=d(:,:,1), &
pa=p(:,:,1:1), &
xa=xpy2(:,:,1:1), &
nmtx=1, &
threshold=1.0d-14, &
infos=infos)
call dftclean(infos)
call measure_time(print_total=1, log_unit=iw)
call flush(iw)
end if
! Compute 2e gradient
call tdhf_2e_grad(basis, infos, d(:,:,1:1), p(:,:,1:1), xpy2(:,:,1:1), xmy2(:,:,1:1))
call print_gradient(infos)
! Print timings
call measure_time(print_total=1, log_unit=iw)
close(iw)
end subroutine tdhf_gradient