module tdhf_gradient_mod use precision, only: dp use grd2, only: grd2_driver, grd2_compute_data_t use basis_tools, only: basis_set, bas_norm_matrix, build_cart_density use constants, only: HARMONIC_ACTIVE, NUM_CART_BF use types, only: information use io_constants, only: iw !############################################################################### implicit none !############################################################################### character(len=*), parameter :: module_name = "tdhf_gradient_mod" !############################################################################### public tdhf_gradient !############################################################################### type, extends(grd2_compute_data_t) :: grd2_tdhf_compute_data_t real(kind=dp), pointer :: d2(:,:,:) => null() real(kind=dp), pointer :: p2(:,:,:) => null() real(kind=dp), pointer :: xpy2(:,:,:) => null() real(kind=dp), pointer :: xmy2(:,:,:) => null() ! Cartesian-effective (bfnrm-folded) copies + offsets for HARMONIC_ACTIVE. real(kind=dp), allocatable :: d_cart(:,:), p_cart(:,:), xpy_cart(:,:), xmy_cart(:,:) integer, allocatable :: cart_off(:) integer :: nbf = 0 contains procedure :: init => grd2_tdhf_compute_data_t_init procedure :: clean => grd2_tdhf_compute_data_t_clean procedure :: get_density => grd2_tdhf_compute_data_t_get_density procedure :: build_cart => grd2_tdhf_build_cart end type contains !############################################################################### subroutine tdhf_gradient_C(c_handle) bind(C, name="tdhf_gradient") use c_interop, only: oqp_handle_t, oqp_handle_get_info use strings, only: Cstring implicit none type(oqp_handle_t) :: c_handle type(information), pointer :: inf inf => oqp_handle_get_info(c_handle) call tdhf_gradient(inf) end subroutine tdhf_gradient_c !############################################################################### 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(:,infos%tddft%target_state), 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(:,infos%tddft%target_state), 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 !############################################################################### subroutine tdhf_1e_grad(infos, basis) use oqp_tagarray_driver use types, only: information use basis_tools, only: basis_set use util, only: measure_time use messages, only: show_message, WITH_ABORT use precision, only: dp use constants, only: tol_int use grd1, only: eijden, print_gradient, & grad_nn, grad_ee_overlap, & grad_ee_kinetic, grad_en_hellman_feynman, grad_en_pulay, grad_1e_ecp use mathlib, only: symmetrize_matrix implicit none character(len=*), parameter :: subroutine_name = "tdhf_1e_grad" type(information), intent(inout) :: infos type(basis_set), intent(inout) :: basis real(kind=dp), allocatable :: dens(:) real(kind=dp) :: tol integer :: nbf, nbf2, ok ! tagarray real(kind=dp), pointer :: dmat_a(:), wao(:), td_p(:,:) character(len=*), parameter :: tags_general(*) = (/ character(len=80) :: & OQP_DM_A, OQP_WAO, OQP_TD_P /) nbf = basis%nbf nbf2 = nbf*(nbf+1)/2 tol = tol_int*log(10.0_dp) ! initial memory allocation allocate(dens(nbf2), source=0.0_dp, stat=ok) if(ok/=0) call show_message('Cannot allocate memory', WITH_ABORT) 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_WAO, wao) call tagarray_get_data(infos%dat, OQP_td_p, td_p) associate( grad => infos%atoms%grad & , xyz => infos%atoms%xyz & , zn => infos%atoms%zn - infos%basis%ecp_zn_num & , p => td_p & , w => wao & , urohf => infos%control%scftype>=2 & ) ! Zero out gradient grad = 0.0d0 ! Nuclear repulsion force call grad_nn(infos%atoms, infos%basis%ecp_zn_num) ! Obtain Lagrangian matrix (`dens`) call eijden(dens, nbf, infos) ! Add W matrix: dens = dens - 2*w ! Overlap gradient call grad_ee_overlap(basis, dens, grad, logtol=tol) ! Compute total density matrix, discard Lagrangian dens = dmat_a + 2*p(:,1) ! Hellmann-Feynman force call grad_en_hellman_feynman(basis, xyz, zn, dens, grad, logtol=tol) ! KE gradient call grad_ee_kinetic(basis, dens, grad, logtol=tol) ! Pulay force call grad_en_pulay(basis, xyz, zn, dens, grad, logtol=tol) ! Effective core potential gradient call grad_1e_ecp(infos, basis, xyz, dens, grad, logtol=tol) end associate end subroutine !############################################################################### !> @brief The driver for the two electron gradient subroutine tdhf_2e_grad(basis, infos, d, p, xpy, xmy) use basis_tools, only: basis_set use precision, only: dp use messages, only: show_message, WITH_ABORT use types, only: information implicit none type(information), target, intent(inout) :: infos type(basis_set) :: basis real(kind=dp), contiguous, target :: p(:,:,:), d(:,:,:), xpy(:,:,:), xmy(:,:,:) logical :: urohf real(kind=dp) :: hfscale integer :: ok real(kind=dp), allocatable :: de(:,:) class(grd2_compute_data_t), allocatable :: gcomp urohf = infos%control%scftype >= 2 hfscale = 1.0d0 if(infos%control%hamilton.ge.20) hfscale = infos%dft%hfscale allocate(de(3,ubound(infos%atoms%zn,1)), & source=0.0d0, & stat=ok) if(ok/=0) call show_message('cannot allocate memory', WITH_ABORT) gcomp = grd2_tdhf_compute_data_t( d2 = d & , p2 = p & , xpy2 = xpy & , xmy2 = xmy & , hfscale = hfscale & , nbf = basis%nbf & ) call gcomp%init() select type (gcomp) class is (grd2_tdhf_compute_data_t) call gcomp%build_cart(basis) end select call grd2_driver(infos, basis, de, gcomp) infos%atoms%grad = infos%atoms%grad + de call gcomp%clean() end subroutine !############################################################################### subroutine grd2_tdhf_compute_data_t_init(this) !use messages, only: show_message, WITH_ABORT implicit none class(grd2_tdhf_compute_data_t), target, intent(inout) :: this call this%clean() end subroutine !############################################################################### !> @brief Build Cartesian-effective (bfnrm-folded) copies of the four TDDFT !> gradient densities (ground d, relaxed p, transition X+Y and X-Y) + !> Cartesian offsets, for contraction with Cartesian derivative ERIs under !> HARMONIC_ACTIVE. Transition densities may be non-symmetric; the per-block !> expansion handles that. subroutine grd2_tdhf_build_cart(this, basis) class(grd2_tdhf_compute_data_t), intent(inout) :: this type(basis_set), intent(in) :: basis integer, allocatable :: od(:) integer :: nc if (.not. HARMONIC_ACTIVE) return call tdhf_cart_one(basis, this%d2(:,:,1), this%d_cart, this%cart_off, nc) call tdhf_cart_one(basis, this%p2(:,:,1), this%p_cart, od, nc) call tdhf_cart_one(basis, this%xpy2(:,:,1), this%xpy_cart, od, nc) call tdhf_cart_one(basis, this%xmy2(:,:,1), this%xmy_cart, od, nc) end subroutine grd2_tdhf_build_cart subroutine tdhf_cart_one(basis, m, m_cart, off, nc) type(basis_set), intent(in) :: basis real(kind=dp), intent(in) :: m(:,:) real(kind=dp), allocatable, intent(out) :: m_cart(:,:) integer, allocatable, intent(out) :: off(:) integer, intent(out) :: nc real(kind=dp), allocatable :: tmp(:,:) tmp = m call bas_norm_matrix(tmp, basis%bfnrm, basis%nbf) call build_cart_density(basis, tmp, m_cart, off, nc) end subroutine tdhf_cart_one !############################################################################### subroutine grd2_tdhf_compute_data_t_clean(this) implicit none class(grd2_tdhf_compute_data_t), target, intent(inout) :: this end subroutine !############################################################################### !> @brief Compute density factors for \Gamma term of 2-electron contribution !> to TD-DFT gradients subroutine grd2_tdhf_compute_data_t_get_density(this, basis, id, dab, dabmax) implicit none class(grd2_tdhf_compute_data_t), target, intent(inout) :: this type(basis_set), intent(in) :: basis integer, intent(in) :: id(4) real(kind=dp), target, intent(out) :: dab(*) real(kind=dp), intent(out) :: dabmax real(kind=dp) :: coul, coul2, coul3, exch, exch2, exch3, exch4 real(kind=dp) :: df1, bfn real(kind=dp) :: coulfact, xcfact integer :: i_, j_, k_, l_ integer :: i, j, k, l integer :: loc(4) integer :: nbf(4) real(kind=dp), pointer :: ab(:,:,:,:) real(kind=dp), pointer :: p(:,:), d(:,:), xpy(:,:), xmy(:,:) logical :: usecart coulfact = 4*this%coulscale xcfact = this%hfscale dabmax = 0 usecart = HARMONIC_ACTIVE if (usecart) then p => this%p_cart; d => this%d_cart; xpy => this%xpy_cart; xmy => this%xmy_cart loc = this%cart_off(id) - 1 nbf = NUM_CART_BF(basis%am(id)) else p => this%p2(:,:,1); d => this%d2(:,:,1) xpy => this%xpy2(:,:,1); xmy => this%xmy2(:,:,1) loc = basis%ao_offset(id) - 1 nbf = basis%naos(id) end if ab(1:nbf(4),1:nbf(3),1:nbf(2),1:nbf(1)) => dab(1:product(nbf)) do i_ = 1, nbf(1) i = loc(1) + i_ do j_ = 1, nbf(2) j = loc(2) + j_ do k_ = 1, nbf(3) k = loc(3) + k_ do l_ = 1, nbf(4) l = loc(4) + l_ coul = & + p(j,i) * d(l,k) & + p(l,k) * d(j,i) coul2 = & + xpy(i,j) * xpy(k,l) coul3 = d(i,j) * d(k,l) coul = 2*coul + 8*coul2 + coul3 df1 = coulfact*coul if (xcfact/=0.0_dp) then exch = & + p(i,k) * d(l,j) & + p(j,k) * d(l,i) & + p(l,i) * d(k,j) & + p(l,j) * d(k,i) exch2 = & + xpy(k,i) * xpy(l,j) & + xpy(l,i) * xpy(k,j) exch3 = & + (xmy(i,l)-xmy(l,i)) * (xmy(j,k)-xmy(k,j)) & + (xmy(j,l)-xmy(l,j)) * (xmy(i,k)-xmy(k,i)) exch4 = d(i,k) * d(j,l) & + d(i,l) * d(j,k) exch = 2*exch + 8*exch2 + 2*exch3 + exch4 df1 = df1 - xcfact*exch end if dabmax = max(dabmax, abs(df1)) bfn = 1.0_dp if (.not. usecart) bfn = product(basis%bfnrm([i,j,k,l])) ab(l_,k_,j_,i_) = df1*bfn end do end do end do end do end subroutine grd2_tdhf_compute_data_t_get_density !############################################################################### end module tdhf_gradient_mod