subroutine process_intersection_derivatives(inter_idx, k, sign_factor, intersection_results, &
intersection_sets, sys_geom, total_gradient, &
total_hessian, bonds, compute_grad, compute_hess, hess_dim, world_comm)
!! Process derivatives for a single intersection fragment
use mqc_result_types, only: calculation_result_t
use mqc_physical_fragment, only: build_fragment_from_atom_list, &
redistribute_cap_gradients, redistribute_cap_hessian
use mqc_gmbe_utils, only: find_fragment_intersection
use mqc_config_parser, only: bond_t
use mqc_error, only: error_t
integer, intent(in) :: inter_idx, k
real(dp), intent(in) :: sign_factor
type(calculation_result_t), intent(in) :: intersection_results(:)
integer, intent(in) :: intersection_sets(:, :)
type(system_geometry_t), intent(in) :: sys_geom
real(dp), intent(inout) :: total_gradient(:, :)
real(dp), intent(inout), optional :: total_hessian(:, :)
type(bond_t), intent(in), optional :: bonds(:)
logical, intent(in) :: compute_grad, compute_hess
integer, intent(in) :: hess_dim
type(comm_t), intent(in), optional :: world_comm
integer :: j, current_n, next_n, n_intersect
integer, allocatable :: current_atoms(:), next_atoms(:), intersect_atoms(:)
real(dp), allocatable :: term_gradient(:, :), term_hessian(:, :)
logical :: has_intersection
type(physical_fragment_t) :: inter_fragment
type(error_t) :: inter_error
! Build the intersection atom list
call get_monomer_atom_list(sys_geom, intersection_sets(1, inter_idx), current_atoms, current_n)
if (current_n > 0) then
do j = 2, k
call get_monomer_atom_list(sys_geom, intersection_sets(j, inter_idx), next_atoms, next_n)
has_intersection = find_fragment_intersection(current_atoms, current_n, &
next_atoms, next_n, &
intersect_atoms, n_intersect)
deallocate (current_atoms, next_atoms)
if (.not. has_intersection) then
current_n = 0
exit
end if
current_n = n_intersect
allocate (current_atoms(current_n))
current_atoms = intersect_atoms
deallocate (intersect_atoms)
end do
end if
if (current_n > 0) then
call build_fragment_from_atom_list(sys_geom, current_atoms, current_n, &
inter_fragment, inter_error, bonds)
if (inter_error%has_error()) then
call logger%error(inter_error%get_full_trace())
if (present(world_comm)) then
call abort_comm(world_comm, 1)
else
error stop "Failed to build intersection fragment in GMBE"
end if
end if
if (compute_grad .and. intersection_results(inter_idx)%has_gradient) then
allocate (term_gradient(3, sys_geom%total_atoms))
term_gradient = 0.0_dp
call redistribute_cap_gradients(inter_fragment, intersection_results(inter_idx)%gradient, &
term_gradient)
total_gradient = total_gradient + sign_factor*term_gradient
deallocate (term_gradient)
end if
if (compute_hess .and. intersection_results(inter_idx)%has_hessian) then
allocate (term_hessian(hess_dim, hess_dim))
term_hessian = 0.0_dp
call redistribute_cap_hessian(inter_fragment, intersection_results(inter_idx)%hessian, &
term_hessian)
total_hessian = total_hessian + sign_factor*term_hessian
deallocate (term_hessian)
end if
call inter_fragment%destroy()
else
call logger%warning("GMBE intersection has no atoms; skipping derivatives")
end if
if (allocated(current_atoms)) deallocate (current_atoms)
end subroutine process_intersection_derivatives