From da345dc869fe49be43df5add03a3dd80501316b5 Mon Sep 17 00:00:00 2001
From: Michele Nottoli <michele.nottoli@gmail.com>
Date: Wed, 5 Nov 2025 13:35:47 +0100
Subject: [PATCH 1/6] Workflow.

---
 src/ddx_cosmo.f90     | 62 +++++++++++++++++++++++++++++++++++++++++++
 src/ddx_gradients.f90 | 48 ++++++++++++++++++++++++---------
 2 files changed, 98 insertions(+), 12 deletions(-)

diff --git a/src/ddx_cosmo.f90 b/src/ddx_cosmo.f90
index 1bb3b19a..5a597131 100644
--- a/src/ddx_cosmo.f90
+++ b/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.
 !!
diff --git a/src/ddx_gradients.f90 b/src/ddx_gradients.f90
index cd98a4f2..ce01d666 100644
--- a/src/ddx_gradients.f90
+++ b/src/ddx_gradients.f90
@@ -15,24 +15,36 @@ module ddx_gradients
 contains
 
 !> Compute the gradients of the ddCOSMO matrix
-subroutine contract_grad_L(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx)
+subroutine contract_grad_L(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx, dr)
 type(ddx_params_type), intent(in) :: params
     type(ddx_constants_type), intent(in) :: constants
       integer,                         intent(in)    :: isph
-      real(dp),  dimension(constants % nbasis, params % nsph), intent(in)    :: sigma
-      real(dp),  dimension(params % ngrid, params % nsph),       intent(in)    :: xi
-      real(dp),  dimension(constants % nbasis),      intent(inout) :: basloc, vplm
-      real(dp),  dimension(3, constants % nbasis),    intent(inout) :: dbsloc
-      real(dp),  dimension(params % lmax+1),      intent(inout) :: vcos, vsin
-      real(dp),  dimension(3),           intent(inout) :: fx
+      real(dp), dimension(constants % nbasis, params % nsph), intent(in)    :: sigma
+      real(dp), dimension(params % ngrid, params % nsph),       intent(in)    :: xi
+      real(dp), dimension(constants % nbasis), intent(inout) :: basloc, vplm
+      real(dp), dimension(3, constants % nbasis), intent(inout) :: dbsloc
+      real(dp), dimension(params % lmax+1), intent(inout) :: vcos, vsin
+      real(dp), dimension(3), intent(inout) :: fx
+      real(dp), optional, intent(inout) :: dr
 
-      call contract_gradi_Lik(params, constants, isph, sigma, xi(:, isph), basloc, dbsloc, vplm, vcos, vsin, fx )
-      call contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx )
+      logical :: do_dr
+      real(dp) :: dr_local = 0.0d0
+
+      if (present(dr)) then
+          do_dr = .true.
+      else
+          do_dr = .false.
+      end if
+
+      call contract_gradi_Lik(params, constants, isph, sigma, xi(:, isph), basloc, dbsloc, vplm, vcos, vsin, fx, dr_local)
+      call contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx, dr_local)
+
+      if (do_dr) dr = dr_local
 
 end subroutine contract_grad_L
 
 !> Contribution to the gradients of the ddCOSMO matrix
-subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx )
+subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx, dr)
     type(ddx_params_type), intent(in) :: params
     type(ddx_constants_type), intent(in) :: constants
       integer,                         intent(in)    :: isph
@@ -46,6 +58,7 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
       real(dp)  :: vvij, tij, xij, oij, t, fac, fl, f1, f2, f3, beta, tlow, thigh
       real(dp)  :: vij(3), sij(3), alp(3), va(3)
       real(dp), external :: dnrm2
+      real(dp), intent(inout) :: dr
       tlow  = one - pt5*(one - params % se)*params % eta
       thigh = one + pt5*(one + params % se)*params % eta
 
@@ -100,7 +113,7 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
 end subroutine contract_gradi_Lik
 
 !> Contribution to the gradients of the ddCOSMO matrix
-subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx )
+subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx, dr)
     type(ddx_params_type), intent(in) :: params
     type(ddx_constants_type), intent(in) :: constants
       integer,                         intent(in)    :: isph
@@ -118,6 +131,7 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
       real(dp)  :: vji(3), sji(3), alp(3), vb(3), vjk(3), sjk(3), vc(3)
       real(dp) :: rho, ctheta, stheta, cphi, sphi
       real(dp), external :: dnrm2
+      real(dp), intent(inout) :: dr
 
       tlow  = one - pt5*(one - params % se)*params % eta
       thigh = one + pt5*(one + params % se)*params % eta
@@ -207,16 +221,26 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
 end subroutine contract_gradi_Lji
 
 !> Gradient of the characteristic function U
-subroutine contract_grad_U(params, constants, isph, xi, phi, fx )
+subroutine contract_grad_U(params, constants, isph, xi, phi, fx, dr)
     type(ddx_params_type), intent(in) :: params
     type(ddx_constants_type), intent(in) :: constants
       integer,                        intent(in)    :: isph
       real(dp),  dimension(params % ngrid, params % nsph), intent(in)    :: xi, phi
       real(dp),  dimension(3),          intent(inout) :: fx
+      real(dp),  optional, intent(inout) :: dr
       integer :: ig, ji, jsph
       real(dp)  :: vvji, tji, fac, swthr
       real(dp)  :: alp(3), vji(3), sji(3)
       real(dp), external :: dnrm2
+      logical :: do_dr
+      real(dp) :: dr_local
+
+      if (present(dr)) then
+          do_dr = .true.
+      else
+          do_dr = .false.
+      end if
+
       do ig = 1, params % ngrid
         alp = zero
         if (constants % ui(ig,isph) .gt. zero .and. constants % ui(ig,isph).lt.one) then

From bf00b63453f08de4582c8b3a0388cf2b7b488be7 Mon Sep 17 00:00:00 2001
From: Michele Nottoli <michele.nottoli@gmail.com>
Date: Wed, 5 Nov 2025 15:21:45 +0100
Subject: [PATCH 2/6] Testing the gradients on the individual pieces.

---
 src/CMakeLists.txt     |   2 +
 src/test_gradients.f90 | 263 +++++++++++++++++++++++++++++++++++++++++
 2 files changed, 265 insertions(+)
 create mode 100644 src/test_gradients.f90

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 4ad1b06b..ee9e23de 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -75,7 +75,9 @@ if (DDX_LIBRARY)
     endif()
 
     add_executable(ddx_driver ddx_driver.f90)
+    add_executable(test_gradients test_gradients.f90)
     target_link_libraries(ddx_driver PUBLIC ddx)
+    target_link_libraries(test_gradients PUBLIC ddx)
 
     #
     # Installation rules
diff --git a/src/test_gradients.f90 b/src/test_gradients.f90
new file mode 100644
index 00000000..6e075a22
--- /dev/null
+++ b/src/test_gradients.f90
@@ -0,0 +1,263 @@
+!> @copyright (c) 2020-2021 RWTH Aachen. All rights reserved.
+!!
+!! ddX software
+!!
+!! @file src/ddx_driver.f90
+!! Main ddx driver.
+!!
+!! @version 1.0.0
+!! @author Aleksandr Mikhalev
+!! @date 2021-02-25
+
+!> Standalone application of ddX: this program can compute the solvation
+!! energy and forces for a solute made of point charges. All the
+!! relevant steps are here outlined.
+program test_gradients
+use ddx
+use ddx_multipolar_solutes
+use omp_lib
+implicit none
+
+real(dp), parameter :: step = 1d-6
+character(len=255) :: dummy_file_name = ''
+character(len=255) :: fname
+character(len=2047) :: banner
+type(ddx_type) :: ddx_data
+type(ddx_state_type) :: state
+type(ddx_error_type) :: ddx_error
+type(ddx_electrostatics_type) :: electrostatics
+real(dp), allocatable :: psi(:, :), force(:, :), charges(:), &
+    & multipoles(:, :), dr(:), numforce(:,:), numdr(:)
+real(dp) :: tol, esolv, start_time, finish_time
+integer :: i, isph, info, j
+real(dp), allocatable :: x(:,:), s(:,:), tmp_lx(:,:)
+real(dp) :: xpsi, slx, sphi
+real(dp), allocatable :: grad_xpsi(:,:), grad_slx(:,:), &
+    & grad_sphi(:,:), grad_xpsi_num(:,:), grad_slx_num(:,:), &
+    & grad_sphi_num(:,:), dr_xpsi(:), dr_slx(:), dr_sphi(:), &
+    & dr_xpsi_num(:), dr_slx_num(:), dr_sphi_num(:)
+real(dp), external :: ddot
+
+call get_command_argument(1, fname)
+call ddfromfile(fname, ddx_data, tol, charges, ddx_error)
+call check_error(ddx_error)
+
+call allocate_state(ddx_data % params, ddx_data % constants, state, &
+    & ddx_error)
+call check_error(ddx_error)
+
+allocate( &
+    & grad_xpsi(3,ddx_data%params%nsph), &
+    & grad_slx(3,ddx_data%params%nsph), &
+    & grad_sphi(3,ddx_data%params%nsph), &
+    & grad_xpsi_num(3,ddx_data%params%nsph), &
+    & grad_slx_num(3,ddx_data%params%nsph), &
+    & grad_sphi_num(3,ddx_data%params%nsph), &
+    & dr_xpsi(ddx_data%params%nsph), &
+    & dr_slx(ddx_data%params%nsph), &
+    & dr_sphi(ddx_data%params%nsph), &
+    & dr_xpsi_num(ddx_data%params%nsph), &
+    & dr_slx_num(ddx_data%params%nsph), &
+    & dr_sphi_num(ddx_data%params%nsph), &
+    & x(ddx_data%constants%nbasis,ddx_data%params%nsph), &
+    & s(ddx_data%constants%nbasis,ddx_data%params%nsph), &
+    & tmp_lx(ddx_data%constants%nbasis,ddx_data%params%nsph), &
+    & force(3, ddx_data % params % nsph), &
+    & dr(ddx_data % params % nsph), &
+    & numforce(3, ddx_data % params % nsph), &
+    & numdr(ddx_data % params % nsph), &
+    & multipoles(1, ddx_data % params % nsph), stat=info)
+if (info .ne. 0) then
+    write(6, *) "Allocation failed in ddx_driver"
+    stop 1
+end if
+multipoles(1, :) = charges/sqrt4pi
+call multipole_electrostatics(ddx_data % params, ddx_data % constants, &
+    & ddx_data % workspace, multipoles, 0, electrostatics, ddx_error)
+
+allocate(psi(ddx_data % constants % nbasis, ddx_data % params % nsph), &
+    & stat=info)
+if (info .ne. 0) then
+    write(6, *) "Allocation failed in ddx_driver"
+    stop 1
+end if
+call multipole_psi(ddx_data % params, multipoles, 0, psi)
+
+! Reference calculation
+call ddrun(ddx_data,state,electrostatics,psi,tol,esolv,ddx_error,force)
+
+x(:,:) = state%xs(:,:)
+s(:,:) = state%s(:,:)
+! Compute the first Lagrangian term 0.5 X^T Psi
+xpsi = pt5*ddot(ddx_data%constants%n,x,1,psi,1)
+
+! Compute the second Lagrangian term -0.5 S^T L X
+ddx_data%constants%dodiag = .true.
+call lx(ddx_data%params,ddx_data%constants,ddx_data%workspace, &
+    & x,tmp_lx,ddx_error)
+slx = -pt5*ddot(ddx_data%constants%n,s,1,tmp_lx,1)
+
+! Compute the third Lagrangian term - 0.5 S^T Phi
+! Note that Phi is saved as -Phi, so no minus required
+sphi = pt5*ddot(ddx_data%constants%n,s,1,state%phi,1)
+
+write(6,*) xpsi, slx, sphi
+! Check if the Lagrangian formulation returns the energy
+if (abs(esolv - (xpsi + slx + sphi)).gt.1e-10) then
+    write(6, *) "Inconsistency"
+    stop 1
+end if
+
+call sgradlx(ddx_data%params,ddx_data%constants,ddx_data%workspace, &
+    & state,grad_slx,ddx_error)
+
+call sgradphi(ddx_data%params,ddx_data%constants,ddx_data%workspace, &
+    & state,grad_sphi,ddx_error,electrostatics%e_cav,multipoles)
+
+do isph = 1, ddx_data%params%nsph
+    do j = 1, 3
+        ddx_data%params%csph(j,isph) = ddx_data%params%csph(j,isph) - step
+        call displaced_run(ddx_data,multipoles,tol,xpsi,slx,sphi,psi, &
+            & esolv,force,x,s)
+
+        grad_xpsi_num(j, isph) = - xpsi
+        grad_slx_num(j, isph) = - slx
+        grad_sphi_num(j, isph) = - sphi
+
+        ddx_data%params%csph(j,isph) = ddx_data%params%csph(j,isph) + 2.0d0*step
+        call displaced_run(ddx_data,multipoles,tol,xpsi,slx,sphi,psi, &
+            & esolv,force,x,s)
+
+        grad_xpsi_num(j, isph) = grad_xpsi_num(j, isph) + xpsi
+        grad_slx_num(j, isph) = grad_slx_num(j, isph) + slx
+        grad_sphi_num(j, isph) = grad_sphi_num(j, isph) + sphi
+        grad_xpsi_num(j, isph) = grad_xpsi_num(j, isph)/(2.0d0*step)
+        grad_slx_num(j, isph) = grad_slx_num(j, isph)/(2.0d0*step)
+        grad_sphi_num(j, isph) = grad_sphi_num(j, isph)/(2.0d0*step)
+
+        ddx_data%params%csph(j,isph) = ddx_data%params%csph(j,isph) - step
+    end do
+end do
+
+write(6,*) "Difference S^T grad L X", minval(abs(grad_slx_num - grad_slx))
+write(6,*) "Difference S^T grad Phi", maxval(abs(grad_sphi_num - grad_sphi))
+
+deallocate(psi, multipoles, charges, force, numforce, dr, numdr,&
+    & x, s, tmp_lx, &
+    & grad_xpsi, grad_slx, grad_sphi, &
+    & grad_xpsi_num, grad_slx_num, grad_sphi_num, &
+    & dr_xpsi, dr_slx, dr_sphi, &
+    & dr_xpsi_num, dr_slx_num, dr_sphi_num, &
+    & stat=info)
+if (info .ne. 0) then
+    write(6, *) "Deallocation failed in ddx_driver"
+    stop 1
+end if
+call deallocate_electrostatics(electrostatics, ddx_error)
+call deallocate_state(state, ddx_error)
+call deallocate_model(ddx_data, ddx_error)
+
+contains
+
+subroutine displaced_run(ddx_data,multipoles,tol,xpsi,slx,sphi, &
+        & psi,esolv,force,x,s)
+    type(ddx_type), intent(inout) :: ddx_data
+    real(dp), intent(in) :: tol
+    real(dp), intent(out) :: xpsi, slx, sphi
+    type(ddx_type) :: ddx_data2
+    type(ddx_error_type) :: error2
+    type(ddx_state_type) :: state2
+    type(ddx_electrostatics_type) :: electrostatics2
+    real(dp), intent(in) :: multipoles(1,ddx_data%params%nsph)
+    real(dp), intent(out) :: esolv
+    real(dp), intent(out) :: force(3,ddx_data%params%nsph)
+    real(dp), intent(out) :: psi(ddx_data%constants%nbasis, &
+        & ddx_data%params%nsph)
+    real(dp), intent(in) :: &
+        & x(ddx_data%constants%nbasis,ddx_data%params%nsph), &
+        & s(ddx_data%constants%nbasis,ddx_data%params%nsph)
+
+    ! Make a copy of the model
+    call allocate_model(ddx_data%params%nsph,ddx_data%params%csph(1,:), &
+        & ddx_data%params%csph(2,:),ddx_data%params%csph(3,:), &
+        & ddx_data%params%rsph,ddx_data%params%model, &
+        & ddx_data%params%lmax,ddx_data%params%ngrid,ddx_data%params%force, &
+        & ddx_data%params%fmm,ddx_data%params%pm,ddx_data%params%pl, &
+        & ddx_data%params%se,ddx_data%params%eta,ddx_data%params%eps, &
+        & ddx_data%params%kappa,ddx_data%params%matvecmem, &
+        & ddx_data%params%maxiter,ddx_data%params%jacobi_ndiis, &
+        & ddx_data%params%nproc,dummy_file_name,ddx_data2,error2)
+    call check_error(error2)
+
+    call allocate_state(ddx_data2%params,ddx_data2%constants,state2, &
+        & error2)
+    call check_error(error2)
+
+    call multipole_electrostatics(ddx_data2%params,ddx_data2%constants, &
+        & ddx_data2%workspace,multipoles,0,electrostatics2,error2)
+    call check_error(error2)
+
+    call multipole_psi(ddx_data2%params,multipoles,0,psi)
+    call ddrun(ddx_data2,state2,electrostatics2,psi,tol,esolv,ddx_error,force)
+
+    ! Compute the first Lagrangian term 0.5 X^T Psi
+    xpsi = pt5*ddot(ddx_data2%constants%n,x,1,psi,1)
+    ! Compute the second Lagrangian term -0.5 S^T L X
+    ddx_data2%constants%dodiag = .true.
+    call lx(ddx_data2%params,ddx_data2%constants,ddx_data2%workspace, &
+        & x,tmp_lx,error2)
+    slx = -pt5*ddot(ddx_data2%constants%n,s,1,tmp_lx,1)
+    ! Compute the third Lagrangian term - 0.5 S^T Phi
+    ! Note that Phi is saved as -Phi, so no minus required
+    sphi = pt5*ddot(ddx_data2%constants%n,s,1,state2%phi,1)
+
+    call deallocate_electrostatics(electrostatics2,error2)
+    call deallocate_state(state2,error2)
+    call deallocate_model(ddx_data2,error2)
+end subroutine displaced_run
+
+subroutine sgradlx(params, constants, workspace, &
+        & state, force, 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(inout) :: force(3, params % nsph)
+    integer :: isph
+    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))
+    end do
+    force = -pt5*force
+end subroutine sgradlx
+
+subroutine sgradphi(params, constants, workspace, &
+        & state,force,ddx_error,e_cav,multipoles)
+    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(inout) :: force(3, params % nsph)
+    real(dp), intent(in) :: e_cav(3, constants % ncav)
+    real(dp), intent(in) :: multipoles(1, ddx_data % params % nsph)
+    integer :: isph
+    force = zero
+    do isph = 1, params % nsph
+        call contract_grad_u(params, constants, isph, state % sgrid, &
+            & state % phi_grid, force(:, isph))
+    end do
+    force = -pt5*force
+    call zeta_grad(params, constants, state, e_cav, force)
+    call multipole_force_terms(params,constants, &
+        & workspace,state,0,multipoles,force,ddx_error)
+end subroutine sgradphi
+
+end program test_gradients

From e06f5ef4e7b37ac1880b6cc709cc872995643e3b Mon Sep 17 00:00:00 2001
From: robin-dahl <robin.dahl@rwth-aachen.de>
Date: Fri, 7 Nov 2025 14:37:02 +0100
Subject: [PATCH 3/6] Work in progress radii derivative

---
 Input_cosmo.txt        |  57 ++++++++++++++++
 src/ddx_constants.f90  |  18 +++++
 src/ddx_gradients.f90  | 137 ++++++++++++++++++++++++++++++++-----
 src/test_gradients.f90 | 152 ++++++++++++++++++++++++++++++++++++++++-
 4 files changed, 344 insertions(+), 20 deletions(-)
 create mode 100644 Input_cosmo.txt

diff --git a/Input_cosmo.txt b/Input_cosmo.txt
new file mode 100644
index 00000000..129de4ca
--- /dev/null
+++ b/Input_cosmo.txt
@@ -0,0 +1,57 @@
+''          ! Name of log-file for additional information
+1           ! Number of OpenMP cores to be used
+1           ! Model. 1 for COSMO, 2 for PCM and 3 for LPB
+7           ! Maximal degree of modeling spherical harmonics
+302         ! Approximate number of Lebedev grid points
+78.3553     ! Dielectric permittivity constant
+0.0         ! Shift of the regularized characteristic function
+0.1         ! Regularization parameter
+0.0         ! Debye H\"{u}ckel parameter
+0           ! Whether to compute and store sparse matrices (1) or not (0)
+1d-14       ! Relative threshold for the iterative solver
+200         ! Maximum number of iterations
+25          ! Number of DIIS extrapolation points
+1           ! Whether to compute (1) or not (0) forces
+0           ! Whether to use (1) or not (0) the FMM
+20          ! Max degree of multipole spherical harmonics for the FMM
+20          ! Max degree of local spherical harmonics for the FMM
+3                           # nsph
+ 1.09581   0.64052   0.26811   0.00079   0.529177
+ 1.17352   0.53040  -0.07460  -0.00051   0.529177
+-1.25694  -0.93056  -1.33500  -0.00018   0.529177
+
+ 0.33039  -3.37146  -0.33148   0.00013   4.00253
+
+-0.27154  -5.88699  -0.46410  -0.00009   3.63772
+-0.25299  -2.46950   2.04845   0.00073   3.80402
+ 0.21982   0.08761   2.49865  -0.00032   4.00253
+-0.28184   1.30432   4.81168  -0.00073   3.63772
+-0.22781   4.22705  -0.09253   0.00089   3.80402
+-0.20253   2.46931  -3.89829   0.00003   3.80402
+ 0.12910   4.70148  -2.75142  -0.00017   4.00253
+-0.09227   6.17186   1.91355  -0.00056   4.00253
+-0.08846  -3.06677  -4.83502  -0.00048   4.00253
+-0.08265  -3.86502   4.25693  -0.00006   4.00253
+ 0.04393   6.41587  -3.85765  -0.00002   2.99956
+ 0.09224   5.95218   3.08325   1.68217   2.99956
+ 0.09228   5.95188   3.08333  -1.68106   2.99956
+ 0.05299   8.04286   1.05689  -0.00040   2.99956
+ 0.07429  -2.46586  -5.86225   1.68268   2.99956
+ 0.07434  -2.46688  -5.86083  -1.68167   2.99956
+ 0.10225  -5.11800  -4.65665  -0.00014   2.99956
+ 0.08476  -5.05274   4.31652   1.68185   2.99956
+ 0.08479  -5.05304   4.31660  -1.68138   2.99956
+ 0.10651  -2.59079   5.87351   0.00130   2.99956
+12          ! number of spheres
+-0.04192   0.00000   2.29035   1.32281   4.00253
+-0.04192   0.00000   2.29035  -1.32281   4.00253
+-0.04198   0.00000   0.00000  -2.64562   4.00253
+-0.04192   0.00000  -2.29035  -1.32281   4.00253
+-0.04192   0.00000  -2.29035   1.32281   4.00253
+-0.04198   0.00000   0.00000   2.64562   4.00253
+ 0.04193   0.00103   4.05914   2.34326   2.99956
+ 0.04193   0.00103   4.05914  -2.34326   2.99956
+ 0.04197   0.00000   0.00000  -4.68652   2.99956
+ 0.04193  -0.00103  -4.05914  -2.34326   2.99956
+ 0.04193  -0.00103  -4.05914   2.34326   2.99956
+ 0.04197   0.00000   0.00000   4.68652   2.99956
diff --git a/src/ddx_constants.f90 b/src/ddx_constants.f90
index 66577adc..d9f300e6 100644
--- a/src/ddx_constants.f90
+++ b/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,8 @@ 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 +2022,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
diff --git a/src/ddx_gradients.f90 b/src/ddx_gradients.f90
index ce01d666..6cef9fa7 100644
--- a/src/ddx_gradients.f90
+++ b/src/ddx_gradients.f90
@@ -55,8 +55,9 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
       real(dp),  dimension(params % lmax+1),      intent(inout) :: vcos, vsin
       real(dp),  dimension(3),           intent(inout) :: fx
       integer :: ig, ij, jsph, l, ind, m
-      real(dp)  :: vvij, tij, xij, oij, t, fac, fl, f1, f2, f3, beta, tlow, thigh
+      real(dp)  :: vvij, tij, xij, oij, t, fac, fl, f1, fr1, f2, fr2, f3, fr3, beta, tlow, thigh
       real(dp)  :: vij(3), sij(3), alp(3), va(3)
+      real(dp) :: dsij, dtij(3), alp1(3), alp2(3), alp_rad, dsij_rad(3), dtij_rad, alp1_rad, alp2_rad, va_rad
       real(dp), external :: dnrm2
       real(dp), intent(inout) :: dr
       tlow  = one - pt5*(one - params % se)*params % eta
@@ -64,6 +65,7 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
 
       do ig = 1, params % ngrid
         va = zero
+        va_rad = zero
         do ij = constants % inl(isph), constants % inl(isph+1) - 1
           jsph = constants % nl(ij)
           vij  = params % csph(:,isph) + &
@@ -77,8 +79,23 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
           else
               sij = one
           end if
+
+          dsij = 1.0_dp / (tij*params%rsph(jsph)) !scalar
+          dtij(:) = sij(:)/params%rsph(jsph)      !vector
+
+          dsij_rad =  zero !constants%cgrid(:,ig)/vvij &
+            ! & - vij * dot_product(vij, constants%cgrid(:,ig)) / (vvij**3) !vector
+          dtij_rad = dot_product(vij, constants%cgrid(:,ig)) / (vvij* params%rsph(jsph)) !scalar
+
           call dbasis(params, constants, sij, basloc, dbsloc, vplm, vcos, vsin)
           alp  = zero
+          alp1 = zero
+          alp2 = zero
+
+          alp_rad  = zero
+          alp1_rad = zero
+          alp2_rad = zero
+      
           t    = one
           do l = 1, params % lmax
             ind = l*l + l + 1
@@ -87,7 +104,15 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
             do m = -l, l
               f2 = fac*sigma(ind+m,jsph)
               f1 = f2*fl*basloc(ind+m)
-              alp(:) = alp(:) + f1*sij(:) + f2*dbsloc(:,ind+m)
+
+              alp1(:) = f1*dtij
+              alp2(:) = f2*tij* dsij * dbsloc(:,ind+m)
+              alp(:) = alp(:) + alp1(:) + alp2(:)
+
+              alp1_rad = f1*dtij_rad 
+              alp2_rad = f2*tij * dot_product(dsij_rad, dbsloc(:,ind+m))
+              alp_rad = alp_rad + alp1_rad + alp2_rad
+
             end do
             t = t*tij
           end do
@@ -100,15 +125,18 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
             oij = xij
             f2  = zero
           end if
-          f1 = oij/params % rsph(jsph)
+          f1 = oij 
           va(:) = va(:) + f1*alp(:) + beta*f2*constants % zi(:,ig,isph)
+          va_rad = va_rad + f1*alp_rad + beta*f2*constants % zi_dr(ig,isph)
           if (tij .gt. tlow) then
-            f3 = beta*dfsw(tij,params % se,params % eta)/params % rsph(jsph)
+            f3 = beta*dfsw(tij,params % se,params % eta) 
             if (constants % fi(ig,isph).gt.one) f3 = f3/constants % fi(ig,isph)
-            va(:) = va(:) + f3*sij(:)
+            va(:) = va(:) + f3*dtij 
+            va_rad = va_rad + f3*dtij_rad
           end if
         end do
         fx = fx - constants % wgrid(ig)*xi(ig)*va(:)
+        dr = dr - constants % wgrid(ig)*xi(ig)*va_rad
       end do
 end subroutine contract_gradi_Lik
 
@@ -129,6 +157,7 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
       real(dp)  :: vvji, tji, xji, oji, t, fac, fl, f1, f2, beta, di, tlow, thigh
       real(dp)  :: b, g1, g2, vvjk, tjk, xjk
       real(dp)  :: vji(3), sji(3), alp(3), vb(3), vjk(3), sjk(3), vc(3)
+      real(dp) :: dsji, dtji(3), alp1(3), alp2(3), alp_rad, dsji_rad(3), dtji_rad, alp1_rad, alp2_rad, vb_rad, vc_rad
       real(dp) :: rho, ctheta, stheta, cphi, sphi
       real(dp), external :: dnrm2
       real(dp), intent(inout) :: dr
@@ -139,6 +168,10 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
       do ig = 1, params % ngrid
         vb = zero
         vc = zero
+
+        vb_rad = zero
+        vc_rad = zero
+
         do ji = constants % inl(isph), constants % inl(isph+1) - 1
           jsph = constants % nl(ji)
           vji  = params % csph(:,jsph) + &
@@ -152,9 +185,25 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
           else
               sji = one
           end if
+
+          dsji = 1.0_dp/(tji*params % rsph(isph))
+          dtji(:) = -sji/params % rsph(isph)
+
+          dsji_rad = constants%cgrid(:,ig)/vvji &
+            & - vji * dot_product(vji, constants%cgrid(:,ig)) / (vvji**3)
+          dtji_rad = -vvji/(params%rsph(isph)**2) 
+
           call dbasis(params, constants, sji, basloc, dbsloc, vplm, vcos, vsin)
-          alp = zero
-          t   = one
+          alp  = zero
+          alp1 = zero
+          alp2 = zero
+
+          alp_rad  = zero
+          alp1_rad = zero
+          alp2_rad = zero
+
+
+          t    = one
           do l = 1, params % lmax
             ind = l*l + l + 1
             fl  = dble(l)
@@ -162,18 +211,28 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
             do m = -l, l
               f2 = fac*sigma(ind+m,isph)
               f1 = f2*fl*basloc(ind+m)
-              alp = alp + f1*sji + f2*dbsloc(:,ind+m)
+
+              alp1(:) = f1*dtji 
+              alp2(:) = f2*tji* dsji * dbsloc(:,ind+m)
+              alp(:) = alp(:) - alp1(:) + alp2(:)
+
+              alp1_rad = f1*dtji_rad 
+              alp2_rad = f2*tji * dot_product(dsji_rad, dbsloc(:,ind+m))
+              alp_rad = alp_rad - alp1_rad + alp2_rad
+
             end do
             t = t*tji
-          end do
+          end do 
+
           xji = fsw(tji, params % se, params % eta)
           if (constants % fi(ig,jsph).gt.one) then
             oji = xji/constants % fi(ig,jsph)
           else
             oji = xji
           end if
-          f1 = oji/params % rsph(isph)
+          f1 = oji 
           vb = vb + f1*alp*xi(ig,jsph)
+          vb_rad = vb_rad + f1*alp_rad*xi(ig,jsph)
           if (tji .gt. tlow) then
             beta = intmlp(params, constants, tji, sigma(:,isph), basloc)
             if (constants % fi(ig,jsph) .gt. one) then
@@ -204,22 +263,27 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
                 end if
               end do
               if (proc) then
-                g1 = di*di*dfsw(tji, params % se, params % eta)/params % rsph(isph)
+                g1 = di*di*dfsw(tji, params % se, params % eta) 
                 g2 = g1*xi(ig,jsph)*b
-                vc = vc + g2*sji
+                vc = vc - g2*dtji 
+                vc_rad = vc_rad - g2*dtji_rad
               end if
             else
               di  = one
               fac = zero
             end if
-            f2 = (one-fac)*di*dfsw(tji, params % se, params % eta)/params % rsph(isph)
-            vb = vb + f2*xi(ig,jsph)*beta*sji
+            f2 = (one-fac)*di*dfsw(tji, params % se, params % eta) 
+            vb = vb - f2*xi(ig,jsph)*beta*dtji 
+            vb_rad = vb_rad - f2*xi(ig,jsph)*beta*dtji_rad
           end if
         end do
         fx = fx + constants % wgrid(ig)*(vb - vc)
+        dr = dr + constants % wgrid(ig)*(vb_rad - vc_rad)
       end do
 end subroutine contract_gradi_Lji
 
+
+
 !> Gradient of the characteristic function U
 subroutine contract_grad_U(params, constants, isph, xi, phi, fx, dr)
     type(ddx_params_type), intent(in) :: params
@@ -230,7 +294,8 @@ subroutine contract_grad_U(params, constants, isph, xi, phi, fx, dr)
       real(dp),  optional, intent(inout) :: dr
       integer :: ig, ji, jsph
       real(dp)  :: vvji, tji, fac, swthr
-      real(dp)  :: alp(3), vji(3), sji(3)
+      real(dp)  :: alp(3), vji(3), sji(3), dtji(3)
+      real(dp) :: dtji_rad, alp_rad
       real(dp), external :: dnrm2
       logical :: do_dr
       real(dp) :: dr_local
@@ -243,8 +308,10 @@ subroutine contract_grad_U(params, constants, isph, xi, phi, fx, dr)
 
       do ig = 1, params % ngrid
         alp = zero
+        alp_rad = zero
         if (constants % ui(ig,isph) .gt. zero .and. constants % ui(ig,isph).lt.one) then
           alp = alp + phi(ig,isph)*xi(ig,isph)*constants % zi(:,ig,isph)
+          alp_rad = alp_rad + phi(ig,isph)*xi(ig,isph)*constants % zi_dr(ig,isph)
         end if
         do ji = constants % inl(isph), constants % inl(isph+1) - 1
           jsph  = constants % nl(ji)
@@ -257,11 +324,21 @@ subroutine contract_grad_U(params, constants, isph, xi, phi, fx, dr)
           swthr = one + (params % se + 1.d0)*params % eta / 2.d0
           if (tji.lt.swthr .and. tji.gt.swthr-params % eta .and. constants % ui(ig,jsph).gt.zero) then
             sji = vji/vvji
-            fac = - dfsw(tji, params % se, params % eta)/params % rsph(isph)
-            alp = alp + fac*phi(ig,jsph)*xi(ig,jsph)*sji
+            dtji = - sji / params % rsph(isph)
+            dtji_rad = - vvji/(params%rsph(isph)**2) 
+            ! dtji_rad = dot_product(vji, constants%cgrid(:,ig)) / (vvji* params%rsph(isph)) 
+            fac = dfsw(tji, params % se, params % eta) 
+            alp = alp + fac*phi(ig,jsph)*xi(ig,jsph) * dtji 
+            alp_rad = alp_rad + fac*phi(ig,jsph)*xi(ig,jsph) * dtji_rad
           end if
         end do
-        fx = fx - constants % wgrid(ig)*alp
+      !   fx = fx - constants % wgrid(ig)*alp
+        if (present(dr)) then
+         dr = dr - constants % wgrid(ig)*alp_rad
+        else 
+         fx = fx - constants % wgrid(ig)*alp
+      end if
+
       end do
 end subroutine contract_grad_U
 
@@ -2261,4 +2338,28 @@ subroutine zeta_grad(params, constants, state, e_cav, forces)
     end do
 end subroutine zeta_grad
 
+!> Force term: interaction of the external electric field with the zeta
+!! intermediate. This routine is called by the gradients of ddCOSMO
+!! ddPCM and ddLPB.
+subroutine zeta_grad_dr(params, constants, state, e_cav, dr)
+    implicit none
+    type(ddx_params_type), intent(in) :: params
+    type(ddx_constants_type), intent(in) :: constants
+    type(ddx_state_type), intent(inout) :: state
+    real(dp), intent(inout) :: dr(params % nsph)
+    real(dp), intent(in) :: e_cav(3, constants % ncav)
+    ! local variables
+    integer :: icav, isph, igrid
+
+    icav = 0
+    do isph = 1, params % nsph
+        do igrid = 1, params % ngrid
+            if (constants % ui(igrid, isph) .eq. zero) cycle
+            icav = icav + 1
+            dr(isph) = dr(isph) + pt5 &
+                & * dot_product(constants % cgrid(:,igrid), e_cav(:, icav))*state % zeta(icav)
+        end do
+    end do
+end subroutine zeta_grad_dr
+
 end module ddx_gradients
diff --git a/src/test_gradients.f90 b/src/test_gradients.f90
index 6e075a22..d99d0b12 100644
--- a/src/test_gradients.f90
+++ b/src/test_gradients.f90
@@ -38,6 +38,7 @@ program test_gradients
     & dr_xpsi_num(:), dr_slx_num(:), dr_sphi_num(:)
 real(dp), external :: ddot
 
+
 call get_command_argument(1, fname)
 call ddfromfile(fname, ddx_data, tol, charges, ddx_error)
 call check_error(ddx_error)
@@ -71,6 +72,7 @@ program test_gradients
     write(6, *) "Allocation failed in ddx_driver"
     stop 1
 end if
+
 multipoles(1, :) = charges/sqrt4pi
 call multipole_electrostatics(ddx_data % params, ddx_data % constants, &
     & ddx_data % workspace, multipoles, 0, electrostatics, ddx_error)
@@ -83,9 +85,13 @@ program test_gradients
 end if
 call multipole_psi(ddx_data % params, multipoles, 0, psi)
 
+
 ! Reference calculation
+! call draco(ddx_data%params%csph, ddx_data%params%nsph, &
+   !  & ddx_data%params%rsph)
 call ddrun(ddx_data,state,electrostatics,psi,tol,esolv,ddx_error,force)
 
+
 x(:,:) = state%xs(:,:)
 s(:,:) = state%s(:,:)
 ! Compute the first Lagrangian term 0.5 X^T Psi
@@ -108,17 +114,24 @@ program test_gradients
     stop 1
 end if
 
+
 call sgradlx(ddx_data%params,ddx_data%constants,ddx_data%workspace, &
     & state,grad_slx,ddx_error)
 
 call sgradphi(ddx_data%params,ddx_data%constants,ddx_data%workspace, &
     & state,grad_sphi,ddx_error,electrostatics%e_cav,multipoles)
 
+call sdrlx(ddx_data%params,ddx_data%constants,ddx_data%workspace, &
+    & state,grad_slx,dr_slx,ddx_error)
+
+call sdrphi(ddx_data%params,ddx_data%constants,ddx_data%workspace, &
+    & state,grad_sphi,dr_sphi,ddx_error,electrostatics%e_cav,multipoles)
+
 do isph = 1, ddx_data%params%nsph
     do j = 1, 3
         ddx_data%params%csph(j,isph) = ddx_data%params%csph(j,isph) - step
         call displaced_run(ddx_data,multipoles,tol,xpsi,slx,sphi,psi, &
-            & esolv,force,x,s)
+             & esolv,force,x,s)
 
         grad_xpsi_num(j, isph) = - xpsi
         grad_slx_num(j, isph) = - slx
@@ -139,9 +152,37 @@ program test_gradients
     end do
 end do
 
-write(6,*) "Difference S^T grad L X", minval(abs(grad_slx_num - grad_slx))
+write(6,*) "Difference S^T grad L X", maxval(abs(grad_slx_num - grad_slx))
 write(6,*) "Difference S^T grad Phi", maxval(abs(grad_sphi_num - grad_sphi))
 
+
+do isph = 1, ddx_data%params%nsph
+        ddx_data%params%rsph(isph) = ddx_data%params%rsph(isph) - step
+        call displaced_run(ddx_data,multipoles,tol,xpsi,slx,sphi,psi, &
+             & esolv,force,x,s)
+
+        dr_xpsi_num( isph) = - xpsi
+        dr_slx_num(isph) = - slx
+        dr_sphi_num( isph) = - sphi
+
+        ddx_data%params%rsph(isph) = ddx_data%params%rsph(isph) + 2.0d0*step
+        call displaced_run(ddx_data,multipoles,tol,xpsi,slx,sphi,psi, &
+            & esolv,force,x,s)
+
+        dr_xpsi_num( isph) = dr_xpsi_num( isph) + xpsi
+        dr_slx_num(isph) = dr_slx_num(isph) + slx
+        dr_sphi_num( isph) = dr_sphi_num(isph) + sphi
+        dr_xpsi_num( isph) = dr_xpsi_num(isph)/(2.0d0*step)
+        dr_slx_num(isph) = dr_slx_num( isph)/(2.0d0*step)
+        dr_sphi_num( isph) = dr_sphi_num(isph)/(2.0d0*step)
+
+        ddx_data%params%rsph(isph) = ddx_data%params%rsph(isph) - step
+end do
+
+write(6,*) "-----" 
+write(6,*) "Difference S^T grad L X", maxval(abs(dr_slx_num - dr_slx))
+write(6,*) "Difference S^T grad Phi", maxval(abs(dr_sphi_num - dr_sphi))
+
 deallocate(psi, multipoles, charges, force, numforce, dr, numdr,&
     & x, s, tmp_lx, &
     & grad_xpsi, grad_slx, grad_sphi, &
@@ -193,6 +234,7 @@ subroutine displaced_run(ddx_data,multipoles,tol,xpsi,slx,sphi, &
         & error2)
     call check_error(error2)
 
+
     call multipole_electrostatics(ddx_data2%params,ddx_data2%constants, &
         & ddx_data2%workspace,multipoles,0,electrostatics2,error2)
     call check_error(error2)
@@ -237,6 +279,30 @@ subroutine sgradlx(params, constants, workspace, &
     force = -pt5*force
 end subroutine sgradlx
 
+subroutine sdrlx(params, constants, workspace, &
+        & state, 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(inout) :: force(3, params % nsph)
+    real(dp), intent(inout) :: dr(params % nsph)
+    integer :: isph
+    force = zero
+    dr = 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=dr(isph))
+    end do
+    dr = -pt5*dr
+    force = -pt5*force
+end subroutine sdrlx
+
 subroutine sgradphi(params, constants, workspace, &
         & state,force,ddx_error,e_cav,multipoles)
     implicit none
@@ -260,4 +326,86 @@ subroutine sgradphi(params, constants, workspace, &
         & workspace,state,0,multipoles,force,ddx_error)
 end subroutine sgradphi
 
+
+subroutine sdrphi(params, constants, workspace, &
+        & state,force,dr,ddx_error,e_cav,multipoles)
+    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(inout) :: force(3, params % nsph)
+    real(dp), intent(inout) :: dr(params % nsph)
+    real(dp), intent(in) :: e_cav(3, constants % ncav)
+    real(dp), intent(in) :: multipoles(1, ddx_data % params % nsph)
+    integer :: isph
+    dr = zero
+    do isph = 1, params % nsph
+        call contract_grad_u(params, constants, isph, state % sgrid, &
+            & state % phi_grid, force(:, isph), dr=dr(isph))
+    end do
+    dr = -pt5*dr
+    call zeta_grad_dr(params, constants, state, e_cav, dr)
+end subroutine sdrphi
+
+
+
+
+! subroutine draco(xyz, nat, rvdw)
+!    real(dp), intent(in) :: xyz(:, :)
+!    real(dp), intent(inout) :: rvdw(:)
+!    integer, intent(in) :: nat
+!                                        !    x/y/z,  I,  A
+!    real(dp), allocatable :: cn(:), dcndr(:,:,:)
+
+!    allocate(cn(nat))
+
+!    call ncoord_dexp(nat, xyz, cn)
+
+!    rvdw = 7.0_dp !+ cn / 2.0_dp
+
+! end subroutine draco
+
+
+
+! subroutine ncoord_dexp(nat, xyz, cn)
+
+!    !> Molecular structure data
+!    integer, intent(in) :: nat   
+!    real(dp), intent(in), dimension(3, nat) :: xyz
+
+!    !> Error function coordination number.
+!    real(dp), intent(out) :: cn(:)
+
+
+!    integer :: iat, jat, itr
+!    real(dp) :: r2, r1, rc, rij(3)
+
+!    real(dp) :: kcn
+
+!    kcn = 0.5_dp
+
+!    cn(:) = 0.0_dp
+
+!    do iat = 1, nat
+!       do jat = 1, nat
+
+!          if (iat /= jat) then
+
+!             rij = xyz(:, iat) - xyz(:, jat)
+!             r2 = sum(rij**2)
+!             r1 = sqrt(r2)
+
+!             cn(iat) = cn(iat) + exp(-kcn * (1.0_dp - r1))
+
+
+!          end if 
+
+!       end do
+!    end do
+
+! end subroutine ncoord_dexp
+
+
 end program test_gradients

From 4695d6543a6162c6eb0af349ed444122b93ae4db Mon Sep 17 00:00:00 2001
From: Michele Nottoli <michele.nottoli@gmail.com>
Date: Fri, 7 Nov 2025 14:59:23 +0100
Subject: [PATCH 4/6] Introduced jacobian in ddCOSMO derivatives.

---
 src/ddx_gradients.f90 | 67 ++++++++++++++++++++++++++++++++++---------
 1 file changed, 53 insertions(+), 14 deletions(-)

diff --git a/src/ddx_gradients.f90 b/src/ddx_gradients.f90
index 6cef9fa7..a256b67d 100644
--- a/src/ddx_gradients.f90
+++ b/src/ddx_gradients.f90
@@ -60,6 +60,7 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
       real(dp) :: dsij, dtij(3), alp1(3), alp2(3), alp_rad, dsij_rad(3), dtij_rad, alp1_rad, alp2_rad, va_rad
       real(dp), external :: dnrm2
       real(dp), intent(inout) :: dr
+      real(dp) :: sjac(3,3), qij
       tlow  = one - pt5*(one - params % se)*params % eta
       thigh = one + pt5*(one + params % se)*params % eta
 
@@ -83,6 +84,19 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
           dsij = 1.0_dp / (tij*params%rsph(jsph)) !scalar
           dtij(:) = sij(:)/params%rsph(jsph)      !vector
 
+          ! build the jacobian of sik
+          sjac = zero
+          sjac(1,1) = one
+          sjac(2,2) = one
+          sjac(3,3) = one
+          qij = one/vvij
+          do icomp = 1, 3
+              do jcomp = 1, 3
+              sjac(icomp,jcomp) = qij*(sjac(icomp,jcomp) &
+                  & - sij(icomp)*sij(jcomp))
+              end do
+          end do
+
           dsij_rad =  zero !constants%cgrid(:,ig)/vvij &
             ! & - vij * dot_product(vij, constants%cgrid(:,ig)) / (vvij**3) !vector
           dtij_rad = dot_product(vij, constants%cgrid(:,ig)) / (vvij* params%rsph(jsph)) !scalar
@@ -95,7 +109,7 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
           alp_rad  = zero
           alp1_rad = zero
           alp2_rad = zero
-      
+
           t    = one
           do l = 1, params % lmax
             ind = l*l + l + 1
@@ -106,12 +120,18 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
               f1 = f2*fl*basloc(ind+m)
 
               alp1(:) = f1*dtij
-              alp2(:) = f2*tij* dsij * dbsloc(:,ind+m)
+              alp2(1) = sjac(1,1)*dbsloc(1,ind+m) + &
+                  & sjac(1,2)*dbsloc(2,ind+m) + sjac(1,3)*dbsloc(3,ind+m)
+              alp2(2) = sjac(2,1)*dbsloc(1,ind+m) + &
+                  & sjac(2,2)*dbsloc(2,ind+m) + sjac(2,3)*dbsloc(3,ind+m)
+              alp2(3) = sjac(3,1)*dbsloc(1,ind+m) + &
+                  & sjac(3,2)*dbsloc(2,ind+m) + sjac(3,3)*dbsloc(3,ind+m)
+              alp2(:) = alp2(:)*tij*f2
               alp(:) = alp(:) + alp1(:) + alp2(:)
 
               alp1_rad = f1*dtij_rad 
               alp2_rad = f2*tij * dot_product(dsij_rad, dbsloc(:,ind+m))
-              alp_rad = alp_rad + alp1_rad + alp2_rad
+              alp_rad = alp_rad - alp1_rad - alp2_rad
 
             end do
             t = t*tij
@@ -161,6 +181,7 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
       real(dp) :: rho, ctheta, stheta, cphi, sphi
       real(dp), external :: dnrm2
       real(dp), intent(inout) :: dr
+      real(dp) :: sjac(3,3), qji
 
       tlow  = one - pt5*(one - params % se)*params % eta
       thigh = one + pt5*(one + params % se)*params % eta
@@ -186,6 +207,19 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
               sji = one
           end if
 
+          ! build the jacobian of sji
+          sjac = zero
+          sjac(1,1) = - one
+          sjac(2,2) = - one
+          sjac(3,3) = - one
+          qji = one/vvji
+          do icomp = 1, 3
+              do jcomp = 1, 3
+                  sjac(icomp,jcomp) = qji*(sjac(icomp,jcomp) &
+                      & + sji(icomp)*sji(jcomp))
+              end do
+          end do
+
           dsji = 1.0_dp/(tji*params % rsph(isph))
           dtji(:) = -sji/params % rsph(isph)
 
@@ -202,7 +236,6 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
           alp1_rad = zero
           alp2_rad = zero
 
-
           t    = one
           do l = 1, params % lmax
             ind = l*l + l + 1
@@ -212,17 +245,24 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
               f2 = fac*sigma(ind+m,isph)
               f1 = f2*fl*basloc(ind+m)
 
-              alp1(:) = f1*dtji 
-              alp2(:) = f2*tji* dsji * dbsloc(:,ind+m)
-              alp(:) = alp(:) - alp1(:) + alp2(:)
+              alp1(:) = f1*dtji
+              alp2(1) = sjac(1,1)*dbsloc(1,ind+m) + &
+                  & sjac(1,2)*dbsloc(2,ind+m) + sjac(1,3)*dbsloc(3,ind+m)
+              alp2(2) = sjac(2,1)*dbsloc(1,ind+m) + &
+                  & sjac(2,2)*dbsloc(2,ind+m) + sjac(2,3)*dbsloc(3,ind+m)
+              alp2(3) = sjac(3,1)*dbsloc(1,ind+m) + &
+                  & sjac(3,2)*dbsloc(2,ind+m) + sjac(3,3)*dbsloc(3,ind+m)
+              alp2(:) = alp2(:)*tji*f2
 
-              alp1_rad = f1*dtji_rad 
+              alp(:) = alp(:) - alp1(:) - alp2(:)
+
+              alp1_rad = f1*dtji_rad
               alp2_rad = f2*tji * dot_product(dsji_rad, dbsloc(:,ind+m))
               alp_rad = alp_rad - alp1_rad + alp2_rad
 
             end do
             t = t*tji
-          end do 
+          end do
 
           xji = fsw(tji, params % se, params % eta)
           if (constants % fi(ig,jsph).gt.one) then
@@ -230,9 +270,8 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
           else
             oji = xji
           end if
-          f1 = oji 
-          vb = vb + f1*alp*xi(ig,jsph)
-          vb_rad = vb_rad + f1*alp_rad*xi(ig,jsph)
+          vb = vb + oji*alp*xi(ig,jsph)
+          vb_rad = vb_rad + oji*alp_rad*xi(ig,jsph)
           if (tji .gt. tlow) then
             beta = intmlp(params, constants, tji, sigma(:,isph), basloc)
             if (constants % fi(ig,jsph) .gt. one) then
@@ -265,7 +304,7 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
               if (proc) then
                 g1 = di*di*dfsw(tji, params % se, params % eta) 
                 g2 = g1*xi(ig,jsph)*b
-                vc = vc - g2*dtji 
+                vc = vc - g2*dtji
                 vc_rad = vc_rad - g2*dtji_rad
               end if
             else
@@ -273,7 +312,7 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
               fac = zero
             end if
             f2 = (one-fac)*di*dfsw(tji, params % se, params % eta) 
-            vb = vb - f2*xi(ig,jsph)*beta*dtji 
+            vb = vb - f2*xi(ig,jsph)*beta*dtji
             vb_rad = vb_rad - f2*xi(ig,jsph)*beta*dtji_rad
           end if
         end do

From 7367a003d7ec30883b2f4bd649545f99148711a0 Mon Sep 17 00:00:00 2001
From: Michele Nottoli <michele.nottoli@gmail.com>
Date: Mon, 10 Nov 2025 14:04:19 +0100
Subject: [PATCH 5/6] Works.

---
 Input_cosmo.txt        | 57 ----------------------------------
 src/ddx_constants.f90  |  3 +-
 src/ddx_gradients.f90  | 69 ++++++++++++++++++++++--------------------
 src/test_gradients.f90 |  8 ++---
 4 files changed, 43 insertions(+), 94 deletions(-)
 delete mode 100644 Input_cosmo.txt

diff --git a/Input_cosmo.txt b/Input_cosmo.txt
deleted file mode 100644
index 129de4ca..00000000
--- a/Input_cosmo.txt
+++ /dev/null
@@ -1,57 +0,0 @@
-''          ! Name of log-file for additional information
-1           ! Number of OpenMP cores to be used
-1           ! Model. 1 for COSMO, 2 for PCM and 3 for LPB
-7           ! Maximal degree of modeling spherical harmonics
-302         ! Approximate number of Lebedev grid points
-78.3553     ! Dielectric permittivity constant
-0.0         ! Shift of the regularized characteristic function
-0.1         ! Regularization parameter
-0.0         ! Debye H\"{u}ckel parameter
-0           ! Whether to compute and store sparse matrices (1) or not (0)
-1d-14       ! Relative threshold for the iterative solver
-200         ! Maximum number of iterations
-25          ! Number of DIIS extrapolation points
-1           ! Whether to compute (1) or not (0) forces
-0           ! Whether to use (1) or not (0) the FMM
-20          ! Max degree of multipole spherical harmonics for the FMM
-20          ! Max degree of local spherical harmonics for the FMM
-3                           # nsph
- 1.09581   0.64052   0.26811   0.00079   0.529177
- 1.17352   0.53040  -0.07460  -0.00051   0.529177
--1.25694  -0.93056  -1.33500  -0.00018   0.529177
-
- 0.33039  -3.37146  -0.33148   0.00013   4.00253
-
--0.27154  -5.88699  -0.46410  -0.00009   3.63772
--0.25299  -2.46950   2.04845   0.00073   3.80402
- 0.21982   0.08761   2.49865  -0.00032   4.00253
--0.28184   1.30432   4.81168  -0.00073   3.63772
--0.22781   4.22705  -0.09253   0.00089   3.80402
--0.20253   2.46931  -3.89829   0.00003   3.80402
- 0.12910   4.70148  -2.75142  -0.00017   4.00253
--0.09227   6.17186   1.91355  -0.00056   4.00253
--0.08846  -3.06677  -4.83502  -0.00048   4.00253
--0.08265  -3.86502   4.25693  -0.00006   4.00253
- 0.04393   6.41587  -3.85765  -0.00002   2.99956
- 0.09224   5.95218   3.08325   1.68217   2.99956
- 0.09228   5.95188   3.08333  -1.68106   2.99956
- 0.05299   8.04286   1.05689  -0.00040   2.99956
- 0.07429  -2.46586  -5.86225   1.68268   2.99956
- 0.07434  -2.46688  -5.86083  -1.68167   2.99956
- 0.10225  -5.11800  -4.65665  -0.00014   2.99956
- 0.08476  -5.05274   4.31652   1.68185   2.99956
- 0.08479  -5.05304   4.31660  -1.68138   2.99956
- 0.10651  -2.59079   5.87351   0.00130   2.99956
-12          ! number of spheres
--0.04192   0.00000   2.29035   1.32281   4.00253
--0.04192   0.00000   2.29035  -1.32281   4.00253
--0.04198   0.00000   0.00000  -2.64562   4.00253
--0.04192   0.00000  -2.29035  -1.32281   4.00253
--0.04192   0.00000  -2.29035   1.32281   4.00253
--0.04198   0.00000   0.00000   2.64562   4.00253
- 0.04193   0.00103   4.05914   2.34326   2.99956
- 0.04193   0.00103   4.05914  -2.34326   2.99956
- 0.04197   0.00000   0.00000  -4.68652   2.99956
- 0.04193  -0.00103  -4.05914  -2.34326   2.99956
- 0.04193  -0.00103  -4.05914   2.34326   2.99956
- 0.04197   0.00000   0.00000   4.68652   2.99956
diff --git a/src/ddx_constants.f90 b/src/ddx_constants.f90
index d9f300e6..d4503ab7 100644
--- a/src/ddx_constants.f90
+++ b/src/ddx_constants.f90
@@ -994,7 +994,8 @@ subroutine constants_geometry_init(params, constants, ddx_error)
                         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))
+                            & constants % zi_dr(igrid, isph) &
+                            & + vv*dot_product(v, constants % cgrid(:, igrid))
                     end if
                 end if
             enddo
diff --git a/src/ddx_gradients.f90 b/src/ddx_gradients.f90
index a256b67d..ff0f0f07 100644
--- a/src/ddx_gradients.f90
+++ b/src/ddx_gradients.f90
@@ -28,7 +28,7 @@ subroutine contract_grad_L(params, constants, isph, sigma, xi, basloc, dbsloc, v
       real(dp), optional, intent(inout) :: dr
 
       logical :: do_dr
-      real(dp) :: dr_local = 0.0d0
+      real(dp) :: dr_local
 
       if (present(dr)) then
           do_dr = .true.
@@ -36,6 +36,8 @@ subroutine contract_grad_L(params, constants, isph, sigma, xi, basloc, dbsloc, v
           do_dr = .false.
       end if
 
+      dr_local = 0.0d0
+
       call contract_gradi_Lik(params, constants, isph, sigma, xi(:, isph), basloc, dbsloc, vplm, vcos, vsin, fx, dr_local)
       call contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx, dr_local)
 
@@ -67,7 +69,7 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
       do ig = 1, params % ngrid
         va = zero
         va_rad = zero
-        do ij = constants % inl(isph), constants % inl(isph+1) - 1
+       do ij = constants % inl(isph), constants % inl(isph+1) - 1
           jsph = constants % nl(ij)
           vij  = params % csph(:,isph) + &
               & params % rsph(isph)*constants % cgrid(:,ig) - &
@@ -97,9 +99,10 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
               end do
           end do
 
-          dsij_rad =  zero !constants%cgrid(:,ig)/vvij &
-            ! & - vij * dot_product(vij, constants%cgrid(:,ig)) / (vvij**3) !vector
-          dtij_rad = dot_product(vij, constants%cgrid(:,ig)) / (vvij* params%rsph(jsph)) !scalar
+          dsij_rad = constants%cgrid(:,ig)/vvij &
+             & - vij * dot_product(vij, constants%cgrid(:,ig)) / (vvij**3)
+          dtij_rad = dot_product(vij, constants%cgrid(:,ig)) &
+              & / (vvij*params%rsph(jsph))
 
           call dbasis(params, constants, sij, basloc, dbsloc, vplm, vcos, vsin)
           alp  = zero
@@ -120,18 +123,20 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
               f1 = f2*fl*basloc(ind+m)
 
               alp1(:) = f1*dtij
+              alp1_rad = f1*dtij_rad
+
               alp2(1) = sjac(1,1)*dbsloc(1,ind+m) + &
-                  & sjac(1,2)*dbsloc(2,ind+m) + sjac(1,3)*dbsloc(3,ind+m)
+                  & sjac(1,2)*dbsloc(2,ind+m)+sjac(1,3)*dbsloc(3,ind+m)
               alp2(2) = sjac(2,1)*dbsloc(1,ind+m) + &
-                  & sjac(2,2)*dbsloc(2,ind+m) + sjac(2,3)*dbsloc(3,ind+m)
+                  & sjac(2,2)*dbsloc(2,ind+m)+sjac(2,3)*dbsloc(3,ind+m)
               alp2(3) = sjac(3,1)*dbsloc(1,ind+m) + &
-                  & sjac(3,2)*dbsloc(2,ind+m) + sjac(3,3)*dbsloc(3,ind+m)
+                  & sjac(3,2)*dbsloc(2,ind+m)+sjac(3,3)*dbsloc(3,ind+m)
               alp2(:) = alp2(:)*tij*f2
-              alp(:) = alp(:) + alp1(:) + alp2(:)
 
-              alp1_rad = f1*dtij_rad 
-              alp2_rad = f2*tij * dot_product(dsij_rad, dbsloc(:,ind+m))
-              alp_rad = alp_rad - alp1_rad - alp2_rad
+              alp2_rad = f2*tij*dot_product(dsij_rad, dbsloc(:,ind+m))
+
+              alp(:) = alp(:) + alp1(:) + alp2(:)
+              alp_rad = alp_rad + alp1_rad + alp2_rad
 
             end do
             t = t*tij
@@ -145,13 +150,13 @@ subroutine contract_gradi_Lik(params, constants, isph, sigma, xi, basloc, dbsloc
             oij = xij
             f2  = zero
           end if
-          f1 = oij 
+          f1 = oij
           va(:) = va(:) + f1*alp(:) + beta*f2*constants % zi(:,ig,isph)
           va_rad = va_rad + f1*alp_rad + beta*f2*constants % zi_dr(ig,isph)
           if (tij .gt. tlow) then
-            f3 = beta*dfsw(tij,params % se,params % eta) 
+            f3 = beta*dfsw(tij,params % se,params % eta)
             if (constants % fi(ig,isph).gt.one) f3 = f3/constants % fi(ig,isph)
-            va(:) = va(:) + f3*dtij 
+            va(:) = va(:) + f3*dtij
             va_rad = va_rad + f3*dtij_rad
           end if
         end do
@@ -223,9 +228,8 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
           dsji = 1.0_dp/(tji*params % rsph(isph))
           dtji(:) = -sji/params % rsph(isph)
 
-          dsji_rad = constants%cgrid(:,ig)/vvji &
-            & - vji * dot_product(vji, constants%cgrid(:,ig)) / (vvji**3)
-          dtji_rad = -vvji/(params%rsph(isph)**2) 
+          dsji_rad = zero
+          dtji_rad = -tji/params%rsph(isph)
 
           call dbasis(params, constants, sji, basloc, dbsloc, vplm, vcos, vsin)
           alp  = zero
@@ -246,19 +250,20 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
               f1 = f2*fl*basloc(ind+m)
 
               alp1(:) = f1*dtji
+              alp1_rad = f1*dtji_rad
+
               alp2(1) = sjac(1,1)*dbsloc(1,ind+m) + &
-                  & sjac(1,2)*dbsloc(2,ind+m) + sjac(1,3)*dbsloc(3,ind+m)
+                  & sjac(1,2)*dbsloc(2,ind+m)+sjac(1,3)*dbsloc(3,ind+m)
               alp2(2) = sjac(2,1)*dbsloc(1,ind+m) + &
-                  & sjac(2,2)*dbsloc(2,ind+m) + sjac(2,3)*dbsloc(3,ind+m)
+                  & sjac(2,2)*dbsloc(2,ind+m)+sjac(2,3)*dbsloc(3,ind+m)
               alp2(3) = sjac(3,1)*dbsloc(1,ind+m) + &
-                  & sjac(3,2)*dbsloc(2,ind+m) + sjac(3,3)*dbsloc(3,ind+m)
+                  & sjac(3,2)*dbsloc(2,ind+m)+sjac(3,3)*dbsloc(3,ind+m)
               alp2(:) = alp2(:)*tji*f2
 
-              alp(:) = alp(:) - alp1(:) - alp2(:)
+              alp2_rad = f2*tji*dot_product(dsji_rad,dbsloc(:,ind+m))
 
-              alp1_rad = f1*dtji_rad
-              alp2_rad = f2*tji * dot_product(dsji_rad, dbsloc(:,ind+m))
-              alp_rad = alp_rad - alp1_rad + alp2_rad
+              alp(:) = alp(:) + alp1(:) + alp2(:)
+              alp_rad = alp_rad + alp1_rad + alp2_rad
 
             end do
             t = t*tji
@@ -302,7 +307,7 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
                 end if
               end do
               if (proc) then
-                g1 = di*di*dfsw(tji, params % se, params % eta) 
+                g1 = di*di*dfsw(tji, params % se, params % eta)
                 g2 = g1*xi(ig,jsph)*b
                 vc = vc - g2*dtji
                 vc_rad = vc_rad - g2*dtji_rad
@@ -311,13 +316,13 @@ subroutine contract_gradi_Lji(params, constants, isph, sigma, xi, basloc, dbsloc
               di  = one
               fac = zero
             end if
-            f2 = (one-fac)*di*dfsw(tji, params % se, params % eta) 
-            vb = vb - f2*xi(ig,jsph)*beta*dtji
-            vb_rad = vb_rad - f2*xi(ig,jsph)*beta*dtji_rad
+            f2 = (one-fac)*di*dfsw(tji, params % se, params % eta)
+            vb = vb + f2*xi(ig,jsph)*beta*dtji
+            vb_rad = vb_rad + f2*xi(ig,jsph)*beta*dtji_rad
           end if
         end do
-        fx = fx + constants % wgrid(ig)*(vb - vc)
-        dr = dr + constants % wgrid(ig)*(vb_rad - vc_rad)
+        fx = fx - constants % wgrid(ig)*(vb + vc)
+        dr = dr - constants % wgrid(ig)*(vb_rad + vc_rad)
       end do
 end subroutine contract_gradi_Lji
 
@@ -364,7 +369,7 @@ subroutine contract_grad_U(params, constants, isph, xi, phi, fx, dr)
           if (tji.lt.swthr .and. tji.gt.swthr-params % eta .and. constants % ui(ig,jsph).gt.zero) then
             sji = vji/vvji
             dtji = - sji / params % rsph(isph)
-            dtji_rad = - vvji/(params%rsph(isph)**2) 
+            dtji_rad = - vvji/(params%rsph(isph)**2)
             ! dtji_rad = dot_product(vji, constants%cgrid(:,ig)) / (vvji* params%rsph(isph)) 
             fac = dfsw(tji, params % se, params % eta) 
             alp = alp + fac*phi(ig,jsph)*xi(ig,jsph) * dtji 
diff --git a/src/test_gradients.f90 b/src/test_gradients.f90
index d99d0b12..dd7ba393 100644
--- a/src/test_gradients.f90
+++ b/src/test_gradients.f90
@@ -173,15 +173,15 @@ program test_gradients
         dr_slx_num(isph) = dr_slx_num(isph) + slx
         dr_sphi_num( isph) = dr_sphi_num(isph) + sphi
         dr_xpsi_num( isph) = dr_xpsi_num(isph)/(2.0d0*step)
-        dr_slx_num(isph) = dr_slx_num( isph)/(2.0d0*step)
+        dr_slx_num(isph) = dr_slx_num(isph)/(2.0d0*step)
         dr_sphi_num( isph) = dr_sphi_num(isph)/(2.0d0*step)
 
         ddx_data%params%rsph(isph) = ddx_data%params%rsph(isph) - step
 end do
 
-write(6,*) "-----" 
-write(6,*) "Difference S^T grad L X", maxval(abs(dr_slx_num - dr_slx))
-write(6,*) "Difference S^T grad Phi", maxval(abs(dr_sphi_num - dr_sphi))
+write(6,*) "-----"
+write(6,*) "Difference S^T dr L X", maxval(abs(dr_slx_num - dr_slx))
+write(6,*) "Difference S^T dr Phi", maxval(abs(dr_sphi_num - dr_sphi))
 
 deallocate(psi, multipoles, charges, force, numforce, dr, numdr,&
     & x, s, tmp_lx, &

From b1819f25bc642e2e51eca40d8a41340f9bccf920 Mon Sep 17 00:00:00 2001
From: Michele Nottoli <michele.nottoli@gmail.com>
Date: Mon, 30 Mar 2026 16:46:29 +0200
Subject: [PATCH 6/6] Moved test.

---
 src/CMakeLists.txt                |  2 -
 tests/CMakeLists.txt              |  1 +
 tests/Input_cosmo_small.txt       | 21 +++++++
 {src => tests}/test_gradients.f90 | 94 ++++++++++---------------------
 4 files changed, 51 insertions(+), 67 deletions(-)
 create mode 100644 tests/Input_cosmo_small.txt
 rename {src => tests}/test_gradients.f90 (89%)

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ee9e23de..4ad1b06b 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -75,9 +75,7 @@ if (DDX_LIBRARY)
     endif()
 
     add_executable(ddx_driver ddx_driver.f90)
-    add_executable(test_gradients test_gradients.f90)
     target_link_libraries(ddx_driver PUBLIC ddx)
-    target_link_libraries(test_gradients PUBLIC ddx)
 
     #
     # Installation rules
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 49209055..ac0a7c4a 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -92,6 +92,7 @@ ddx_add_test("ddlpb_esolv.f90" "${CMAKE_CURRENT_SOURCE_DIR}/data/ddlpb_force.txt
 ddx_add_test("matrix_derivatives.f90" "${CMAKE_CURRENT_SOURCE_DIR}/data/ddlpb_force.txt")
 ddx_add_test("matrix_adjoint.f90" "${CMAKE_CURRENT_SOURCE_DIR}/data/ddlpb_force.txt")
 ddx_add_test("matrix_solvers.f90" "${CMAKE_CURRENT_SOURCE_DIR}/data/ddlpb_force.txt")
+ddx_add_test("test_gradients.f90" "${CMAKE_CURRENT_SOURCE_DIR}/Input_cosmo_small.txt")
 ddx_add_test("m2l.f90")
 ddx_add_test("multipolar_solutes.f90")
 ddx_add_test("error.f90")
diff --git a/tests/Input_cosmo_small.txt b/tests/Input_cosmo_small.txt
new file mode 100644
index 00000000..284a3d05
--- /dev/null
+++ b/tests/Input_cosmo_small.txt
@@ -0,0 +1,21 @@
+''          ! Name of log-file for additional information
+1           ! Number of OpenMP cores to be used
+1           ! Model. 1 for COSMO, 2 for PCM and 3 for LPB
+7           ! Maximal degree of modeling spherical harmonics
+302         ! Approximate number of Lebedev grid points
+78.3553     ! Dielectric permittivity constant
+0.0         ! Shift of the regularized characteristic function
+0.1         ! Regularization parameter
+0.0         ! Debye H\"{u}ckel parameter
+0           ! Whether to compute and store sparse matrices (1) or not (0)
+1d-14       ! Relative threshold for the iterative solver
+200         ! Maximum number of iterations
+25          ! Number of DIIS extrapolation points
+1           ! Whether to compute (1) or not (0) forces
+0           ! Whether to use (1) or not (0) the FMM
+20          ! Max degree of multipole spherical harmonics for the FMM
+20          ! Max degree of local spherical harmonics for the FMM
+3                           # nsph
+ 1.09581   0.0 0.0 0.0 2.0
+ 1.17352   0.99 0.0 0.0 3.0
+ 0.33039  -3.37146  -0.33148   0.00013   4.00253
diff --git a/src/test_gradients.f90 b/tests/test_gradients.f90
similarity index 89%
rename from src/test_gradients.f90
rename to tests/test_gradients.f90
index dd7ba393..95c99b3e 100644
--- a/src/test_gradients.f90
+++ b/tests/test_gradients.f90
@@ -31,12 +31,13 @@ program test_gradients
 real(dp) :: tol, esolv, start_time, finish_time
 integer :: i, isph, info, j
 real(dp), allocatable :: x(:,:), s(:,:), tmp_lx(:,:)
-real(dp) :: xpsi, slx, sphi
+real(dp) :: xpsi, slx, sphi, diff
 real(dp), allocatable :: grad_xpsi(:,:), grad_slx(:,:), &
     & grad_sphi(:,:), grad_xpsi_num(:,:), grad_slx_num(:,:), &
     & grad_sphi_num(:,:), dr_xpsi(:), dr_slx(:), dr_sphi(:), &
     & dr_xpsi_num(:), dr_slx_num(:), dr_sphi_num(:)
 real(dp), external :: ddot
+real(dp), parameter :: threshold = 1e-8
 
 
 call get_command_argument(1, fname)
@@ -109,7 +110,8 @@ program test_gradients
 
 write(6,*) xpsi, slx, sphi
 ! Check if the Lagrangian formulation returns the energy
-if (abs(esolv - (xpsi + slx + sphi)).gt.1e-10) then
+diff = abs(esolv - (xpsi + slx + sphi))
+if (diff.gt.threshold) then
     write(6, *) "Inconsistency"
     stop 1
 end if
@@ -152,8 +154,19 @@ program test_gradients
     end do
 end do
 
-write(6,*) "Difference S^T grad L X", maxval(abs(grad_slx_num - grad_slx))
-write(6,*) "Difference S^T grad Phi", maxval(abs(grad_sphi_num - grad_sphi))
+diff = maxval(abs(grad_slx_num - grad_slx))
+write(6,*) "Difference S^T grad L X", diff
+if (diff.gt.threshold) then
+    write(6, *) "Inconsistency"
+    stop 1
+end if
+
+diff = maxval(abs(grad_sphi_num - grad_sphi))
+write(6,*) "Difference S^T grad Phi", diff
+if (diff.gt.threshold) then
+    write(6, *) "Inconsistency"
+    stop 1
+end if
 
 
 do isph = 1, ddx_data%params%nsph
@@ -180,8 +193,18 @@ program test_gradients
 end do
 
 write(6,*) "-----"
-write(6,*) "Difference S^T dr L X", maxval(abs(dr_slx_num - dr_slx))
-write(6,*) "Difference S^T dr Phi", maxval(abs(dr_sphi_num - dr_sphi))
+diff = maxval(abs(dr_slx_num - dr_slx))
+write(6,*) "Difference S^T dr L X", diff
+if (diff.gt.threshold) then
+    write(6, *) "Inconsistency"
+    stop 1
+end if
+diff = maxval(abs(dr_sphi_num - dr_sphi))
+write(6,*) "Difference S^T dr Phi", diff
+if (diff.gt.threshold) then
+    write(6, *) "Inconsistency"
+    stop 1
+end if
 
 deallocate(psi, multipoles, charges, force, numforce, dr, numdr,&
     & x, s, tmp_lx, &
@@ -349,63 +372,4 @@ subroutine sdrphi(params, constants, workspace, &
     call zeta_grad_dr(params, constants, state, e_cav, dr)
 end subroutine sdrphi
 
-
-
-
-! subroutine draco(xyz, nat, rvdw)
-!    real(dp), intent(in) :: xyz(:, :)
-!    real(dp), intent(inout) :: rvdw(:)
-!    integer, intent(in) :: nat
-!                                        !    x/y/z,  I,  A
-!    real(dp), allocatable :: cn(:), dcndr(:,:,:)
-
-!    allocate(cn(nat))
-
-!    call ncoord_dexp(nat, xyz, cn)
-
-!    rvdw = 7.0_dp !+ cn / 2.0_dp
-
-! end subroutine draco
-
-
-
-! subroutine ncoord_dexp(nat, xyz, cn)
-
-!    !> Molecular structure data
-!    integer, intent(in) :: nat   
-!    real(dp), intent(in), dimension(3, nat) :: xyz
-
-!    !> Error function coordination number.
-!    real(dp), intent(out) :: cn(:)
-
-
-!    integer :: iat, jat, itr
-!    real(dp) :: r2, r1, rc, rij(3)
-
-!    real(dp) :: kcn
-
-!    kcn = 0.5_dp
-
-!    cn(:) = 0.0_dp
-
-!    do iat = 1, nat
-!       do jat = 1, nat
-
-!          if (iat /= jat) then
-
-!             rij = xyz(:, iat) - xyz(:, jat)
-!             r2 = sum(rij**2)
-!             r1 = sqrt(r2)
-
-!             cn(iat) = cn(iat) + exp(-kcn * (1.0_dp - r1))
-
-
-!          end if 
-
-!       end do
-!    end do
-
-! end subroutine ncoord_dexp
-
-
 end program test_gradients