EEQ-BC Implementation (#56)

This PR contains the implementation of the EEQ-BC charge model (see
[here](https://doi.org/10.1063/5.0268978)).

The implementation contains full support for
- electrostatic energy with gradients and stress tensor
- atomic partial charges with gradients and lattice vector derivatives
- both molecular and periodic calculations

EEQ-BC is implemented as an alternative model. The general solution of
the linear system remains unchanged. The new capacitance matrix is
handled in the cache and modifies only the final A-matrix and b-vector.

Additionally, we added explicit routines to calculate the EEQ and EEQ-BC
charges, avoiding the through DFT-D4.

---------

Co-authored-by: Polt <[email protected]>
Co-authored-by: Leopold Seidler <[email protected]>
Co-authored-by: Marvin Friede <[email protected]>

Author: 3 people

README.md

@@ -5,7 +5,7 @@
 [![CI](https://github.com/grimme-lab/multicharge/workflows/CI/badge.svg)](https://github.com/grimme-lab/multicharge/actions)
 [![codecov](https://codecov.io/gh/grimme-lab/multicharge/branch/main/graph/badge.svg)](https://codecov.io/gh/grimme-lab/multicharge)
 
-Electronegativity equilibration model for atomic partial charges.
+Electronegativity equilibration models for atomic partial charges.
 
 
 ## Installation
@@ -135,6 +135,33 @@ multicharge.git = "https://github.com/grimme-lab/multicharge"
 
 For an overview over all command line arguments use the ``--help`` argument or checkout the [``multicharge(1)``](man/multicharge.1.adoc) manpage.
 
+## Citation
+
+For the electronegativity equilibration model (EEQ): 
+
+Eike Caldeweyher, Sebastian Ehlert, Andreas Hansen, Hagen Neugebauer, Sebastian Spicher, Christoph Bannwarth and Stefan Grimme, *J. Chem Phys*, **2019**, 150, 154122.
+DOI: [10.1063/1.5090222](https://doi.org/10.1063/1.5090222)
+chemrxiv: [10.26434/chemrxiv.7430216](https://doi.org/10.26434/chemrxiv.7430216.v2)
+
+... the EEQ extention to Fr, Ra, and the full Actinide series: 
+
+Lukas Wittmann, Igor Gordiy, Marvin Friede, Benjamin Helmich-Paris, Stefan Grimme, Andreas Hansen and Markus Bursch, *Phys. Chem. Chem. Phys.*, **2024**, 26, 21379-21394.
+DOI: [10.1039/D4CP01514B](10.1039/D4CP01514B)
+
+... the periodic EEQ implementation: 
+
+Eike Caldeweyher, Jan-Michael Mewes, Sebastian Ehlert and Stefan Grimme, *Phys. Chem. Chem. Phys.*, **2020**, 22, 8499-8512.
+DOI: [10.1039/D0CP00502A](https://doi.org/10.1039/D0CP00502A)
+chemrxiv: [10.26434/chemrxiv.10299428](https://doi.org/10.26434/chemrxiv.10299428.v1)
+
+<br>
+
+For the bond capacity electronegativity equilibration charge model (EEQ<sub>BC</sub>): 
+
+Thomas Froitzheim, Marcel Müller, Andreas Hansen, and Stefan Grimme, *J. Chem. Phys.*, **2025**, 162, 214109.
+DOI: [10.1039/10.1063/5.0268978](https://doi.org/10.1063/5.0268978)
+chemrxiv: [10.26434/chemrxiv-2025-1nxwg](https://doi.org/10.26434/chemrxiv-2025-1nxwg)
+
 
 ## License
 

app/main.f90

@@ -14,38 +14,42 @@
 ! limitations under the License.
 
 program main
-   use, intrinsic :: iso_fortran_env, only : output_unit, error_unit, input_unit
-   use mctc_env, only : error_type, fatal_error, get_argument, wp
-   use mctc_io, only : structure_type, read_structure, filetype, get_filetype
-   use multicharge, only : mchrg_model_type, new_eeq2019_model, &
-      & write_ascii_model, write_ascii_properties, write_ascii_results, &
-      & get_multicharge_version
-   use multicharge_output, only : json_results
+   use, intrinsic :: iso_fortran_env, only: output_unit, error_unit, input_unit
+   use mctc_env, only: error_type, fatal_error, get_argument, wp
+   use mctc_io, only: structure_type, read_structure, filetype, get_filetype
+   use mctc_cutoff, only: get_lattice_points
+   use multicharge, only: mchrg_model_type, mchrg_model, new_eeq2019_model, &
+      & new_eeqbc2025_model, get_multicharge_version, &
+      & write_ascii_model, write_ascii_properties, write_ascii_results
+   use multicharge_output, only: json_results
    implicit none
    character(len=*), parameter :: prog_name = "multicharge"
    character(len=*), parameter :: json_output = "multicharge.json"
 
    character(len=:), allocatable :: input, chargeinput
    integer, allocatable :: input_format
-   integer :: stat, unit
+   integer :: stat, unit, model_id
    type(error_type), allocatable :: error
    type(structure_type) :: mol
-   type(mchrg_model_type) :: model
+   class(mchrg_model_type), allocatable :: model
    logical :: grad, json, exist
    real(wp), parameter :: cn_max = 8.0_wp, cutoff = 25.0_wp
-   real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :)
+   real(wp), allocatable :: cn(:), rcov(:), trans(:, :)
+   real(wp), allocatable :: qloc(:)
+   real(wp), allocatable :: dcndr(:, :, :), dcndL(:, :, :), dqlocdr(:, :, :), dqlocdL(:, :, :)
    real(wp), allocatable :: energy(:), gradient(:, :), sigma(:, :)
-   real(wp), allocatable :: qvec(:), dqdr(:, :, :), dqdL(:, :, :)
+   real(wp), allocatable :: qvec(:)
+   real(wp), allocatable :: dqdr(:, :, :), dqdL(:, :, :)
    real(wp), allocatable :: charge
 
-   call get_arguments(input, input_format, grad, charge, json, error)
+   call get_arguments(input, model_id, input_format, grad, charge, json, error)
    if (allocated(error)) then
       write(error_unit, '(a)') error%message
       error stop
    end if
 
    if (input == "-") then
-      if (.not.allocated(input_format)) input_format = filetype%xyz
+      if (.not. allocated(input_format)) input_format = filetype%xyz
       call read_structure(mol, input_unit, input_format, error)
    else
       call read_structure(mol, input, error, input_format)
@@ -76,36 +80,47 @@ program main
       end if
    end if
 
-   call new_eeq2019_model(mol, model, error)
-   if(allocated(error)) then
+   if (model_id == mchrg_model%eeq2019) then
+      call new_eeq2019_model(mol, model, error)
+   else if (model_id == mchrg_model%eeqbc2025) then
+      call new_eeqbc2025_model(mol, model, error)
+   else
+      call fatal_error(error, "Invalid model was choosen.")
+   end if
+   if (allocated(error)) then
       write(error_unit, '(a)') error%message
       error stop
    end if
 
    call write_ascii_model(output_unit, mol, model)
 
-   allocate(cn(mol%nat))
-   if (grad) then
-      allocate(dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat))
-   end if
-
-   call model%ncoord%get_cn(mol, cn, dcndr, dcndL)
-
    allocate(energy(mol%nat), qvec(mol%nat))
    energy(:) = 0.0_wp
+
+   allocate(cn(mol%nat), qloc(mol%nat))
    if (grad) then
       allocate(gradient(3, mol%nat), sigma(3, 3))
       gradient(:, :) = 0.0_wp
       sigma(:, :) = 0.0_wp
+
+      allocate(dqdr(3, mol%nat, mol%nat), dqdL(3, 3, mol%nat))
+      dqdr(:, :, :) = 0.0_wp
+      dqdL(:, :, :) = 0.0_wp
+
+      allocate(dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat))
+      allocate(dqlocdr(3, mol%nat, mol%nat), dqlocdL(3, 3, mol%nat))
    end if
 
-   call model%solve(mol, error, cn, dcndr, dcndL, energy, gradient, sigma, &
-      & qvec, dqdr, dqdL)
+   call get_lattice_points(mol%periodic, mol%lattice, model%ncoord%cutoff, trans)
+   call model%ncoord%get_coordination_number(mol, trans, cn, dcndr, dcndL)
+   call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL)
+   call model%solve(mol, error, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL, &
+      & energy, gradient, sigma, qvec, dqdr, dqdL)
+
    if (allocated(error)) then
       write(error_unit, '(a)') error%message
       error stop
    end if
-   
 
    call write_ascii_properties(output_unit, mol, model, cn, qvec)
    call write_ascii_results(output_unit, mol, energy, gradient, sigma)
@@ -115,12 +130,11 @@ program main
       call json_results(unit, "  ", energy=sum(energy), gradient=gradient, charges=qvec, cn=cn)
       close(unit)
       write(output_unit, '(a)') &
-         "[Info] JSON dump of results written to '"// json_output //"'"
+         "[Info] JSON dump of results written to '"//json_output//"'"
    end if
 
 contains
 
-
 subroutine help(unit)
    integer, intent(in) :: unit
 
@@ -133,7 +147,8 @@ subroutine help(unit)
       "higher multipole moments", &
       ""
 
-   write(unit, '(2x, a, t25, a)') &
+   write(unit, '(2x, a, t35, a)') &
+      "-m, -model, --model <model>", "Choose the charge model", &
       "-i, -input, --input <format>", "Hint for the format of the input file", &
       "-c, -charge, --charge <value>", "Set the molecular charge", &
       "-g, -grad, --grad", "Evaluate molecular gradient and virial", &
@@ -145,7 +160,6 @@ subroutine help(unit)
 
 end subroutine help
 
-
 subroutine version(unit)
    integer, intent(in) :: unit
    character(len=:), allocatable :: version_string
@@ -156,12 +170,15 @@ subroutine version(unit)
 
 end subroutine version
 
-
-subroutine get_arguments(input, input_format, grad, charge, json, error)
+subroutine get_arguments(input, model_id, input_format, grad, charge, &
+   & json, error)
 
    !> Input file name
    character(len=:), allocatable :: input
 
+   !> ID of choosen model type
+   integer, intent(out) :: model_id
+
    !> Input file format
    integer, allocatable, intent(out) :: input_format
 
@@ -180,6 +197,7 @@ subroutine get_arguments(input, input_format, grad, charge, json, error)
    integer :: iarg, narg, iostat
    character(len=:), allocatable :: arg
 
+   model_id = mchrg_model%eeq2019
    grad = .false.
    json = .false.
    iarg = 0
@@ -195,24 +213,39 @@ subroutine get_arguments(input, input_format, grad, charge, json, error)
          call version(output_unit)
          stop
       case default
-         if (.not.allocated(input)) then
+         if (.not. allocated(input)) then
             call move_alloc(arg, input)
             cycle
          end if
          call fatal_error(error, "Too many positional arguments present")
          exit
+      case("-m", "-model", "--model")
+         iarg = iarg + 1
+         call get_argument(iarg, arg)
+         if (.not. allocated(arg)) then
+            call fatal_error(error, "Missing argument for model")
+            exit
+         end if
+         if (arg == "eeq2019" .or. arg == "eeq") then
+            model_id = mchrg_model%eeq2019
+         else if (arg == "eeqbc2025" .or. arg == "eeqbc") then
+            model_id = mchrg_model%eeqbc2025
+         else
+            call fatal_error(error, "Invalid model")
+            exit
+         end if
       case("-i", "-input", "--input")
          iarg = iarg + 1
          call get_argument(iarg, arg)
-         if (.not.allocated(arg)) then
+         if (.not. allocated(arg)) then
             call fatal_error(error, "Missing argument for input format")
             exit
          end if
          input_format = get_filetype("."//arg)
       case("-c", "-charge", "--charge")
          iarg = iarg + 1
          call get_argument(iarg, arg)
-         if (.not.allocated(arg)) then
+         if (.not. allocated(arg)) then
             call fatal_error(error, "Missing argument for charge")
             exit
          end if
@@ -229,8 +262,8 @@ subroutine get_arguments(input, input_format, grad, charge, json, error)
       end select
    end do
 
-   if (.not.allocated(input)) then
-      if (.not.allocated(error)) then
+   if (.not. allocated(input)) then
+      if (.not. allocated(error)) then
          call help(output_unit)
          error stop
       end if

man/multicharge.1.adoc

@@ -14,18 +14,21 @@ Electronegativity equilibration model for atomic partial charges.
 
 == Options
 
+*-m, -model, --model* _model_::
+Choose the charge model (eeq or eeqbc)
+
 *-i, -input, --input* _format_::
 Hint for the format of the input file
 
 *-c, -charge, --charge* _value_::
 Provide the molecular charge
 
-*-j, -json, --json*::
-Provide output in JSON format to the file 'multicharge.json'
-
 *-g, -grad, --grad*::
 Evaluate molecular gradient and virial
 
+*-j, -json, --json*::
+Provide output in JSON format to the file 'multicharge.json'
+
 *-v, -version, --version*::
 Print program version and exit
 

src/multicharge.f90

@@ -14,11 +14,11 @@
 ! limitations under the License.
 
 module multicharge
-   use multicharge_cutoff, only : get_lattice_points
+   use multicharge_charge, only : get_charges, get_eeq_charges, get_eeqbc_charges
    use multicharge_model, only : mchrg_model_type
    use multicharge_output, only : write_ascii_model, write_ascii_properties, &
       & write_ascii_results
-   use multicharge_param, only : new_eeq2019_model
+   use multicharge_param, only : new_eeq2019_model, new_eeqbc2025_model, mchrg_model
    use multicharge_version, only : get_multicharge_version
    implicit none
    public

src/multicharge/CMakeLists.txt

@@ -13,21 +13,23 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+add_subdirectory("model")
 add_subdirectory("param")
 
 set(dir "${CMAKE_CURRENT_SOURCE_DIR}")
 
 list(
   APPEND srcs
   "${dir}/blas.F90"
-  "${dir}/cutoff.f90"
+  "${dir}/charge.f90"
   "${dir}/ewald.f90"
   "${dir}/lapack.F90"
-  "${dir}/model.F90"
+  "${dir}/model.f90"
   "${dir}/output.f90"
   "${dir}/param.f90"
   "${dir}/version.f90"
   "${dir}/wignerseitz.f90"
+  "${dir}/cache.f90"
 )
 
 set(srcs "${srcs}" PARENT_SCOPE)

src/multicharge/cache.f90

@@ -0,0 +1,46 @@
+! This file is part of multicharge.
+! SPDX-Identifier: Apache-2.0
+!
+! Licensed under the Apache License, Version 2.0 (the "License");
+! you may not use this file except in compliance with the License.
+! You may obtain a copy of the License at
+!
+!     http://www.apache.org/licenses/LICENSE-2.0
+!
+! Unless required by applicable law or agreed to in writing, software
+! distributed under the License is distributed on an "AS IS" BASIS,
+! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+! See the License for the specific language governing permissions and
+! limitations under the License.
+
+!> @file multicharge/cache.f90
+!> Contains the cache baseclass for the charge models and a container for mutable cache data
+
+!> Cache for charge models
+module multicharge_model_cache
+   use mctc_env, only: wp
+   use mctc_io, only: structure_type
+   use multicharge_wignerseitz, only: wignerseitz_cell_type
+   implicit none
+   private
+
+   type, public :: cache_container
+      !> Mutable data attribute
+      class(*), allocatable :: raw
+   end type cache_container
+
+   !> Cache for the charge model
+   type, abstract, public :: model_cache
+      !> Coordination number array
+      real(wp), allocatable :: cn(:)
+      !> Coordination number gradient w.r.t the positions
+      real(wp), allocatable :: dcndr(:, :, :)
+      !> Coordination number gradient w.r.t the lattice vectors
+      real(wp), allocatable :: dcndL(:, :, :)
+      !> Ewald separation parameter
+      real(wp) :: alpha
+      !> Wigner-Seitz cell
+      type(wignerseitz_cell_type) :: wsc
+   end type model_cache
+
+end module multicharge_model_cache