Radii derivatives

src/ddx_constants.f90

@@ -138,6 +138,9 @@ module ddx_constants
     !> Derivative of the characteristic function U at all grid points of all
     !! spheres. Dimension is (3, ngrid, nsph).
     real(dp), allocatable :: zi(:, :, :)
+    !> Derivative of the characteristic function U at all grid points of all
+    !! spheres wrt to radii. Dimension is (ngrid, nsph).
+    real(dp), allocatable :: zi_dr(:, :)
     !> Number of external Lebedev grid points on a molecular surface.
     integer :: ncav
     !> Number of external Lebedev grid points on each sphere.
@@ -950,6 +953,13 @@ subroutine constants_geometry_init(params, constants, ddx_error)
             return
         endif
         constants % zi = zero
+
+        allocate(constants % zi_dr(params % ngrid, params % nsph), stat=info)
+        if (info .ne. 0) then
+            call update_error(ddx_error, "`zi_dr` allocation failed")
+            return
+        endif
+        constants % zi_dr = zero
     end if
     ! Build arrays fi, ui, zi
     !$omp parallel do default(none) shared(params,constants,swthr) &
@@ -983,6 +993,9 @@ subroutine constants_geometry_init(params, constants, ddx_error)
                             & params % rsph(jsph) / vv
                         constants % zi(:, igrid, isph) = &
                             & constants % zi(:, igrid, isph) + vv*v
+                        constants % zi_dr(igrid, isph) = &
+                            & constants % zi_dr(igrid, isph) &
+                            & + vv*dot_product(v, constants % cgrid(:, igrid))
                     end if
                 end if
             enddo
@@ -2010,6 +2023,12 @@ subroutine constants_free(constants, ddx_error)
             call update_error(ddx_error, "`zi` deallocation failed!")
         end if
     end if
+    if (allocated(constants % zi_dr)) then
+        deallocate(constants % zi_dr, stat=istat)
+        if (istat .ne. 0) then
+            call update_error(ddx_error, "`zi_dr` deallocation failed!")
+        end if
+    end if
     if (allocated(constants % ncav_sph)) then
         deallocate(constants % ncav_sph, stat=istat)
         if (istat .ne. 0) then

src/ddx_cosmo.f90

@@ -239,6 +239,68 @@ subroutine cosmo_solvation_force_terms(params, constants, workspace, &
 
 end subroutine cosmo_solvation_force_terms
 
+!> Compute the solvation term of the forces and the radial derivatives
+!> (solute aspecific). This must be summed to the solute specific term
+!! to get the full forces.
+!!
+!> @ingroup Fortran_interface_ddcosmo
+!! @param[in] params: User specified parameters
+!! @param[in] constants: Precomputed constants
+!! @param[inout] workspace: Preallocated workspaces
+!! @param[inout] state: ddx state (contains solutions and RHSs)
+!! @param[in] e_cav: electric field, size (3, ncav)
+!! @param[inout] force: force term
+!! @param[inout] ddx_error: ddX error
+!!
+subroutine cosmo_solvation_force_dr_terms(params, constants, workspace, &
+        & state, e_cav, force, dr, ddx_error)
+    implicit none
+    type(ddx_params_type), intent(in) :: params
+    type(ddx_constants_type), intent(in) :: constants
+    type(ddx_workspace_type), intent(inout) :: workspace
+    type(ddx_state_type), intent(inout) :: state
+    type(ddx_error_type), intent(inout) :: ddx_error
+    real(dp), intent(in) :: e_cav(3, constants % ncav)
+    real(dp), intent(inout) :: force(3, params % nsph)
+    real(dp), intent(inout) :: dr(params % nsph)
+    ! local variables
+    real(dp) :: start_time, finish_time
+    integer :: isph
+
+    ! dummy operation on unused interface arguments
+    if (ddx_error % flag .eq. 0) continue
+
+    start_time = omp_get_wtime()
+
+    force = zero
+    do isph = 1, params % nsph
+        call contract_grad_l(params, constants, isph, state % xs, &
+            & state % sgrid, workspace % tmp_vylm(:, 1), &
+            & workspace % tmp_vdylm(:, :, 1), workspace % tmp_vplm(:, 1), &
+            & workspace % tmp_vcos(:, 1), workspace % tmp_vsin(:, 1), &
+            & force(:, isph), dr(isph))
+        call contract_grad_u(params, constants, isph, state % sgrid, &
+            & state % phi_grid, force(:, isph), dr(isph))
+    end do
+
+    force = -pt5*force
+
+    call zeta_grad(params, constants, state, e_cav, force)
+
+    finish_time = omp_get_wtime()
+    state % force_time = finish_time - start_time
+
+end subroutine cosmo_solvation_force_dr_terms
+
+!> This routines precomputes the intermediates to be used in the evaluation
+!! of ddCOSMO analytical derivatives.
+!!
+!! @param[in] params: ddx parameters
+!! @param[in] constant: ddx constants
+!! @param[inout] workspace: ddx workspaces
+!! @param[inout] state: ddx state
+!!
+
 !> This routines precomputes the intermediates to be used in the evaluation
 !! of ddCOSMO analytical derivatives.
 !!