Analytical charge Hessian contractions

config/cmake/Findmctc-lib.cmake

@@ -16,7 +16,7 @@
 set(_lib "mctc-lib")
 set(_pkg "MCTCLIB")
 set(_url "https://github.com/grimme-lab/mctc-lib")
-set(_rev "v0.6.1")
+set(_rev "v0.7.0")
 
 if(NOT DEFINED "${_pkg}_FIND_METHOD")
   if(DEFINED "${PROJECT_NAME}-dependency-method")

config/meson.build

@@ -153,7 +153,7 @@ mctc_dep = dependency('mctc-lib', required: false)
 if not mctc_dep.found()
   mctc_prj = subproject(
     'mctc-lib',
-    version: '>=0.5.1',
+    version: '>=0.7.0',
     default_options: [
       'default_library=static',
     ],

fpm.toml

@@ -10,7 +10,7 @@ link = ["lapack", "blas"]
 
 [dependencies]
 mctc-lib.git = "https://github.com/grimme-lab/mctc-lib.git"
-mctc-lib.tag = "v0.6.1"
+mctc-lib.tag = "v0.7.0"
 
 [dev-dependencies]
 mstore.git = "https://github.com/grimme-lab/mstore.git"

src/multicharge.f90

@@ -14,8 +14,9 @@
 ! limitations under the License.
 
 module multicharge
-   use multicharge_charge, only : get_charges, get_eeq_charges, get_eeqbc_charges
+   use multicharge_charge, only : get_charges, get_eeq_charges, get_eeqbc_charges, add_charge_hessian
    use multicharge_model, only : mchrg_model_type
+   use multicharge_model_cache, only : cache_container
    use multicharge_output, only : write_ascii_model, write_ascii_properties, &
       & write_ascii_results
    use multicharge_param, only : new_eeq2019_model, new_eeqbc2025_model, mchrg_model

src/multicharge/cache.f90

@@ -27,6 +27,8 @@ module multicharge_model_cache
    type, public :: cache_container
       !> Mutable data attribute
       class(*), allocatable :: raw
+      !> Optional charge-response data, refreshed by get_charges/solve.
+      real(wp), allocatable :: qvec(:), ainv(:, :), dqdr(:, :, :)
    end type cache_container
 
    !> Cache for the charge model

src/multicharge/charge.f90

@@ -23,18 +23,19 @@ module multicharge_charge
    use mctc_env, only : error_type, wp
    use mctc_io, only : structure_type
    use multicharge_model, only : mchrg_model_type
+   use multicharge_model_cache, only : cache_container
    use multicharge_param, only : new_eeq2019_model, new_eeqbc2025_model
    implicit none
    private
 
-   public :: get_charges, get_eeq_charges, get_eeqbc_charges
+   public :: get_charges, get_eeq_charges, get_eeqbc_charges, add_charge_hessian
 
 
 contains
 
 
 !> Classical electronegativity equilibration charges
-subroutine get_charges(mchrg_model, mol, error, qvec, dqdr, dqdL)
+subroutine get_charges(mchrg_model, mol, error, qvec, dqdr, dqdL, cache)
 
    !> Multicharge model
    class(mchrg_model_type), intent(in) :: mchrg_model
@@ -54,12 +55,16 @@ subroutine get_charges(mchrg_model, mol, error, qvec, dqdr, dqdL)
    !> Derivative of the partial charges w.r.t. strain deformations
    real(wp), intent(out), contiguous, optional :: dqdL(:, :, :)
 
+   !> Retain the solution and Cartesian response for add_charge_hessian.
+   !> Recompute after any change to the geometry, charge, or model.
+   type(cache_container), intent(inout), optional :: cache
+
    logical :: grad
    real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :)
    real(wp), allocatable :: qloc(:), dqlocdr(:, :, :), dqlocdL(:, :, :)
    real(wp), allocatable :: trans(:, :)
 
-   grad = present(dqdr) .and. present(dqdL)
+   grad = (present(dqdr) .and. present(dqdL)) .or. present(cache)
 
    allocate(cn(mol%nat), qloc(mol%nat))
    if (grad) then
@@ -72,11 +77,25 @@ subroutine get_charges(mchrg_model, mol, error, qvec, dqdr, dqdL)
    call mchrg_model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL)
 
    call mchrg_model%solve(mol, error, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL, &
-      & qvec=qvec, dqdr=dqdr, dqdL=dqdL)
+      & qvec=qvec, dqdr=dqdr, dqdL=dqdL, response=cache)
 
 end subroutine get_charges
 
 
+!> Add sum_i weights(i) d2q_i/dR2, holding weights and lattice fixed.
+!> The cache must come from get_charges with the same structure and model.
+subroutine add_charge_hessian(mchrg_model, mol, error, weights, hessian, cache)
+   class(mchrg_model_type), intent(in) :: mchrg_model
+   type(structure_type), intent(in) :: mol
+   type(error_type), allocatable, intent(out) :: error
+   real(wp), intent(in) :: weights(:)
+   real(wp), intent(inout) :: hessian(:, :)
+   type(cache_container), intent(inout) :: cache
+
+   call mchrg_model%add_charge_hessian(mol, error, weights, hessian, cache)
+end subroutine add_charge_hessian
+
+
 !> Obtain charges from electronegativity equilibration model
 subroutine get_eeq_charges(mol, error, qvec, dqdr, dqdL)
 

src/multicharge/model/CMakeLists.txt

@@ -18,6 +18,7 @@ set(dir "${CMAKE_CURRENT_SOURCE_DIR}")
 
 list(
   APPEND srcs
+  "${dir}/hessian.f90"
   "${dir}/eeq.f90"
   "${dir}/eeqbc.f90"
   "${dir}/type.F90"