subroutine oqp_tdhf_z_vector(infos)
use oqp_tagarray_driver
use strings, only: Cstring, fstring
use messages, only: show_message, with_abort
use util, only: measure_time
use tdhf_lib, only: iatogen, mntoia, esum, &
tdhf_unrelaxed_density
use tdhf_sf_lib, only: sfrorhs, &
sfromcal, sfrogen, sfrolhs, &
pcgb, sfropcal, sfrowcal
use dft, only: dft_initialize, dftclean
use mod_dft_gridint_fxc, only: tddft_fxc
use mod_dft_gridint_gxc, only: tddft_gxc
use mathlib, only: symmetrize_matrix, orthogonal_transform
use mod_dft_molgrid, only: dft_grid_t
use mathlib, only: pack_matrix, unpack_matrix
use mathlib, only: triangular_to_full
use tdhf_lib, only: int2_rpagrd_data_t
use pcg_mod
use printing, only: print_module_info
implicit none
character(len=*), parameter :: subroutine_name = "oqp_tdhf_z_vector"
type(basis_set), pointer :: basis
type(information), target, intent(inout) :: infos
type(dft_grid_t), target :: molGrid
type(int2_compute_t), target :: int2_driver
class(int2_fock_data_t), allocatable, target :: int2_data
type(tdhf_cg_data) :: cgdata
type(pcg_t) :: pcg
real(kind=dp), allocatable :: hpp(:,:,:), hpt(:,:,:), hmm(:,:,:), gxp(:,:,:)
real(kind=dp), allocatable, target :: wrk1(:,:), &
rhs(:), xm(:), zvec(:), xminv(:), pa(:,:,:)
real(kind=dp), contiguous, pointer :: ppa(:,:,:,:), &
pxm(:,:), prhs(:,:), wmo(:,:)
logical :: dft, tda
integer :: scf_type, mol_mult
integer :: i, j, iter, ok
integer :: nocc, nvir, lexc, nbf, nbf2
real(kind=dp) :: cnvtol, scale_exch
! tagarray
real(kind=dp), contiguous, pointer :: &
mo_a(:,:), mo_energy_a(:), wao(:), td_p(:,:), td_t(:,:), &
ta(:), xpy(:,:), xmy(:,:), td_energies(:)
character(len=*), parameter :: &
tags_required(*) = [ character(len=80) :: &
OQP_E_MO_A, OQP_VEC_MO_A, OQP_TD_T, OQP_TD_XPY, OQP_TD_XMY, &
OQP_td_energies &
]
! Log file
open (unit=IW, file=infos%log_filename, position="append")
!
call print_module_info('TDHF_Z_Vector','Solving Z-Vector for TDDFT')
mol_mult = infos%mol_prop%mult
scf_type = infos%control%scftype
dft = infos%control%hamilton == 20
tda = infos%tddft%tda
! Load basis set
basis => infos%basis
basis%atoms => infos%atoms
nbf = basis%nbf
nbf2 = nbf*(nbf+1)/2
if (dft) call dft_initialize(infos, basis, molGrid)
! Parameter it should be inputed later
! convergence tolerance in the iterative TD-DFT step.
cnvtol = infos%tddft%zvconv
nocc = infos%mol_prop%nocc
nvir = nbf-nocc
lexc = nocc*nvir
allocate(&
rhs(lexc), &
pa(nbf,nbf,1), &
wrk1(nbf,nbf), &
stat=ok, &
source=0.0_dp)
if( ok/=0 ) call show_message('Cannot allocate memory', with_abort)
call infos%dat%remove_records([character(80) :: OQP_WAO, OQP_TD_P])
call data_has_tags(infos%dat, tags_required, module_name, subroutine_name, WITH_ABORT)
call tagarray_get_data(infos%dat, OQP_E_MO_A, mo_energy_a)
call tagarray_get_data(infos%dat, OQP_VEC_MO_A, mo_a)
call tagarray_get_data(infos%dat, OQP_TD_T, td_t)
call tagarray_get_data(infos%dat, OQP_TD_XPY, xpy)
call tagarray_get_data(infos%dat, OQP_TD_XMY, xmy)
call tagarray_get_data(infos%dat, OQP_td_energies, td_energies)
ta => td_t(:,1)
write(iw,'(/1x,71("-")&
&/19x,"TD-DFT ENERGY GRADIENT CALCULATION"&
&/1x,71("-")/)')
write(iw,fmt='(5x,a/&
&5x,16("-")/&
&5x,a,x,i0,x,f17.10,x,"Hartree"/&
&5x,a,x,i0/&
&5x,a,x,e10.4/&
&5x,a,x,i0)') &
'Z-vector options' &
, 'Target state is', infos%tddft%target_state, infos%mol_energy%energy+td_energies(infos%tddft%target_state) &
, 'Multiplicity is', infos%tddft%mult &
, 'Convergence is', infos%tddft%zvconv &
, 'Maximum iterations is', infos%control%maxit_zv
call flush(iw)
! 1. Compute right-hand side of Z-vector equation
allocate(hpp(nbf,nbf,1), & ! H+[X+Y]
hpt(nbf,nbf,1), & ! H+[T]
hmm(nbf,nbf,1), & ! H-[X-Y]
gxp(nbf,nbf,1), & ! g_xc[X+Y][X+Y]
source=0.0d0)
call tdhf_unrelaxed_density(xmy(:,infos%tddft%target_state), xpy(:,infos%tddft%target_state), mo_a, td_t(:,1), nocc, tda)
! Initialize ERI calculations
call int2_driver%init(basis, infos)
call int2_driver%set_screening()
! Compute H+[X+Y], H+[T], H-[X-Y], and G_xc
call compute_r_terms(infos, basis, int2_driver, molGrid, mo_a, xpy(:,infos%tddft%target_state), &
xmy(:,infos%tddft%target_state), ta, hpp, hpt, hmm, gxp)
! Transform Gxc to MO basis
call orthogonal_transform('n', nbf, mo_a, gxp(:,:,1))
! Assemble RHS in MO basis
prhs(1:nocc,1:nvir) => rhs(1:)
call compute_r_mo(prhs, mo_a, xpy(:,infos%tddft%target_state), xmy(:,infos%tddft%target_state), hpt, hpp, hmm, gxp(:,:,1))
rhs = -rhs
! 2. Initialize CG for Z-vector solution
write(iw,'(/3x,25("-")&
&/6x,"START Z-VECTOR LOOP"&
&/3x,25("-")/)')
call flush(iw)
allocate(xminv(lexc), &
xm(lexc), &
stat=ok, &
source=0.0_dp)
if( ok/=0 ) call show_message('Cannot allocate memory', with_abort)
scale_exch = 1.0_dp
if (dft) scale_exch = infos%dft%HFscale
int2_data = int2_td_data_t(d2=pa, &
int_apb=.true., &
int_amb=.false., &
tamm_dancoff=tda, &
scale_exchange=scale_exch)
pxm(1:nocc, 1:nvir) => xm(1:)
do i = 1, nvir
do j = 1, nocc
pxm(j,i) = mo_energy_a(nocc+i) - mo_energy_a(j)
end do
end do
xminv = 1.0d0/xm
cgdata = tdhf_cg_data( &
infos=infos, int2_driver=int2_driver, &
int2_data=int2_data, molgrid=molgrid, &
wrk=wrk1, mo=mo_a, pa=pa, xm=xm, xminv=xminv, &
nbf=nbf, nocc = nocc, dft = dft &
)
call pcg%init(b=rhs, &
update=compute_apbx, precond=precond, &
dat=cgdata, tol=sqrt(abs(cnvtol)))
write(iw,'(" INITIAL ERROR =",3X,1P,E10.3,1X,"/",1P,E10.3)') pcg%error**2, cnvtol
! Begin CG iterations
do iter = 1, infos%control%maxit_zv
if (pcg%errcode /= PCG_OK) exit
call pcg%step()
write(iw,'(" ITER#",I2," ERROR =",3X,1P,E10.3,1X,"/",1P,E10.3)') &
iter, pcg%error**2, cnvtol
call flush(iw)
end do
select case (pcg%errcode)
case (PCG_CONVERGED)
write(iw,'(/3x,24("-")&
&/6x,"Z-Vector converged"&
&/3x,24("-")/)')
infos%mol_energy%Z_Vector_converged=.true.
case default
write(iw,'(/3x,24("-")&
&/6x,"Z-Vector not converged"&
&/3x,24("-")/)')
infos%mol_energy%Z_Vector_converged=.false.
end select
call flush(iw)
! Save Z-vector and clean up memory
deallocate(xminv, rhs, xm)
call int2_data%clean()
deallocate(int2_data)
allocate(zvec, source=pcg%x)
call pcg%clean()
! 3. Now, compute relaxed energy-weighted difference density matrix W
! Convert Z-vector from MO to AO and assemble P = T+Z
pa = 0
call iatogen(zvec, pa(:,:,1), nocc, nocc)
call symmetrize_matrix(pa(:,:,1), nbf)
pa(:,:,1) = 0.5*pa(:,:,1)
call orthogonal_transform('t', nbf, mo_a, pa(:,:,1))
call unpack_matrix(ta, wrk1) ! T
pa(:,:,1) = pa(:,:,1) + wrk1 ! T+Z
! Store relaxed difference density matrix P to global memory
call infos%dat%reserve_data(OQP_td_p, TA_TYPE_REAL64, nbf2, [nbf2, 1], comment=OQP_td_p_comment)
call tagarray_get_data(infos%dat, OQP_td_p, td_p)
call pack_matrix(pa(:,:,1), td_p(:,1))
! Compute H+[P]
ppa(1:nbf,1:nbf, 1:1, 1:1) => pa
int2_data = int2_rpagrd_data_t(&
xpy=null(), &
xmy=null(), &
t=ppa, &
nspin = 1, &
tamm_dancoff=tda, &
scale_exchange=scale_exch)
call int2_driver%run(int2_data, &
cam=dft.and.infos%dft%cam_flag, &
alpha=infos%dft%cam_alpha, &
beta=infos%dft%cam_beta,&
mu=infos%dft%cam_mu)
select type (int2_data)
type is (int2_rpagrd_data_t)
hpt = int2_data%hpt(:,:,:,1,1)
end select
call int2_data%clean()
deallocate(int2_data)
if (dft) then
pa = pa*2
call tddft_fxc(basis=basis, &
molGrid=molGrid, &
isVecs=.true., &
wf=MO_A, &
fx=hpt(:,:,1:1), &
dx=pa(:,:,1:1), &
nmtx=1, &
!threshold=1.0d-15, &
threshold=0.0d0, &
infos=infos)
end if
! Transform H+[P], H+[X+Y], H-[X-Y] from AO to MO basis
call orthogonal_transform('n', nbf, mo_a, hpt(:,:,1), wrk=wrk1) ! P
call orthogonal_transform('n', nbf, mo_a, hpp(:,:,1), wrk=wrk1) ! X+Y
call orthogonal_transform('n', nbf, mo_a, hmm(:,:,1), wrk=wrk1) ! X-Y
! Calculate W in MO basis
wmo(1:nbf,1:nbf) => wrk1
call compute_w_mo(wmo, xpy(:,infos%tddft%target_state), xmy(:,infos%tddft%target_state), zvec, &
hpt(:,:,1), hpp(:,:,1), hmm(:,:,1), &
mo_energy_a, td_energies(infos%tddft%target_state), nocc, gxp(:,:,1))
! Transform W from MO to AO basis
call orthogonal_transform('t', nbf, mo_a, wmo)
! Store W to global memory
call infos%dat%reserve_data(OQP_WAO, TA_TYPE_REAL64, nbf2, comment=OQP_WAO_comment)
call tagarray_get_data(infos%dat, OQP_WAO, wao)
call pack_matrix(wmo, wao)
wao = wao*0.5_dp
! Cleanup
call int2_driver%clean()
if (dft) call dftclean(infos)
call measure_time(print_total=1, log_unit=iw)
close(iw)
end subroutine oqp_tdhf_z_vector