Implementation of EEQ-BC and unit tests #45

Thomas3R · 2025-08-13T15:03:42Z

This pull request introduces significant improvements and refactoring to the charge model and coordination number logic in the codebase. The main changes include a redesign of the coordination number calculation, the addition of a new charge model (EEQBCModel), and enhanced flexibility for coordination number handling via new base and derived classes. The changes also streamline the interface for charge calculation and improve memory management.

Charge Model and Coordination Number Refactoring

The ChargeModel class now owns a NCoordErf member and delegates the calculation of coordination numbers to a new virtual method get_cn, allowing derived models to customize how coordination numbers are computed. This enables more flexible and extensible charge model implementations.

Addition of New Charge Model

A new EEQBCModel class is introduced, supporting bond capacitance-based charge calculation. This class includes extensive element-specific parameters and methods for local charge, capacitance matrix. It offers analytic partial charges and numerical gradients.

Coordination Number Class Hierarchy Enhancement

The coordination number logic is refactored into a base class NCoordBase and several derived classes (NCoordErf, NCoordErfEN, NCoordErfD4). This supports different counting functions, covalent radii handling, and electronegativity scaling, making the coordination number calculation more modular and adaptable to different charge models.

API and Interface Improvements

Method signatures in the charge and coordination number classes are updated for clarity, consistency, and to better separate concerns. For example, the counting function now takes explicit atomic distances and covalent radii, improving correctness and readability.

Memory Management Fix

A memory leak is fixed in the matrix inversion routine by ensuring that the pivot array is deleted in all error cases, improving code safety and reliability.

src/dftd_eeq.cpp

src/dftd_ncoord.cpp

Copilot

Pull Request Overview

This pull request implements the EEQ-BC (Electronegativity Equilibration with Bond Capacitance) charge model and significantly refactors the coordination number calculation system. The changes introduce a more flexible class hierarchy for coordination number computations and add comprehensive unit tests for the new functionality.

Key changes include:

Implementation of EEQBCModel class with bond capacitance-based charge calculation
Refactoring of coordination number classes into a base-derived hierarchy (NCoordBase, NCoordErf, NCoordErfEN, NCoordErfD4)
Enhanced coordination number calculation with configurable parameters and electronegativity scaling
Addition of comprehensive unit tests for the new EEQ-BC model and parameter validation

Reviewed Changes

Copilot reviewed 11 out of 11 changed files in this pull request and generated 7 comments.

Show a summary per file

File	Description
test/test_multi.h	Header file defining reference values and test function declaration for EEQ-BC model validation
test/test_multi.cpp	Unit test implementation for EEQ-BC parameters and charge calculation functionality
test/molecules.h	Addition of hydroxide ion test molecules with different atom orderings
test/meson.build	Integration of new multi test into the build system
test/main.cpp	Registration of multi test in the main test runner
src/dftd_ncoord.cpp	Major refactoring of coordination number calculation with new class hierarchy and enhanced functionality
src/dftd_eeq.cpp	Implementation of EEQ-BC model and refactoring of charge model interface
include/dftd_ncoord.h	Header definitions for the new coordination number class hierarchy

_{Tip: Customize your code reviews with copilot-instructions.md. Create the file or learn how to get started.}

test/test_multi.cpp

src/dftd_ncoord.cpp

src/dftd_eeq.cpp

marvinfriede · 2025-08-15T11:26:21Z

src/dftd_ncoord.cpp

  } else {
    info = ncoord_base(mol, realIdx, dist);
  }
  if (info != EXIT_SUCCESS) return info;

-  info = cut_coordination_number(cn_max, cn, dcndr, lgrad);
-  if (info != EXIT_SUCCESS) return info;
+  if (cn_max != -1.0) { // cn_max = -1.0 for EEQ-BC for cn and qloc; using NCoordErf and NCoordErfEN


Suggested change

if (cn_max != -1.0) { // cn_max = -1.0 for EEQ-BC for cn and qloc; using NCoordErf and NCoordErfEN

if (cn_max < 0.0) { // cn_max = -1.0 for EEQ-BC for cn and qloc; using NCoordErf and NCoordErfEN

Should also work instead of float comparison.

What is the advantage here? Why not use the exact value that cn_max is set to?

Equality comparison of floats is discouraged due to rounding errors.

marvinfriede · 2025-08-15T11:26:57Z

src/dftd_ncoord.cpp

  } else {
    info = ncoord_base(mol, realIdx, dist);
  }
  if (info != EXIT_SUCCESS) return info;

-  info = cut_coordination_number(cn_max, cn, dcndr, lgrad);
-  if (info != EXIT_SUCCESS) return info;
+  if (cn_max != -1.0) { // cn_max = -1.0 for EEQ-BC for cn and qloc; using NCoordErf and NCoordErfEN


Should also work instead of float comparison.

src/dftd_ncoord.cpp

marvinfriede · 2025-08-15T11:37:37Z

include/dftd_multicharge_param.h

+namespace eeqbc {
+
+// Maximum value of the atom number implemented in the EEQ_BC model
+constexpr int MAXELEMENT = 104;


Might be too new for ORCA's C++ standard. Not sure though.

include/dftd_eeq.h

include/dftd_ncoord.h

marvinfriede · 2025-08-15T11:53:37Z

src/dftd_ncoord.cpp


      f_en = get_en_factor(mol.ATNO(i), mol.ATNO(j));
-      countf = f_en * count_fct(rr);
+      countf = f_en * count_fct(r, rcovij);
      cn(ii) += countf;
      cn(jj) += countf;


No f_directed here?

Good catch, the JJ index should always get the directed factor. Apparently, this is not covered by the eeq-bc unit test. Maybe we can find a system where this makes a difference.

src/dftd_ncoord.cpp

marvinfriede · 2025-08-15T12:57:24Z

src/dftd_eeq.cpp

+int EEQBCModel::get_amat_0d(
+  const TMolecule &mol,
+  const TIVector &realIdx,
+  const TMatrix<double> &dist,
+  TMatrix<double> &Amat
+) const {


No ghost atom handling.

This does not work because the matrices use the compressed format, i.e., only contain real atoms, but mol.NAtoms contains all atoms.

…Adjust EEQ-BC functions for ghost atoms.

marvinfriede · 2025-08-19T07:22:51Z

include/dftd_ncoord.h

      bool lgrad);
-    /**
+          /**


marvinfriede · 2025-08-19T07:23:40Z

include/dftd_ncoord.h

+    //TVector<double> cn;  // coordination number
+    //TMatrix<double> dcndr;  // derivative of the coordination number


Remove if not needed anymore.

marvinfriede · 2025-08-19T07:37:06Z

src/dftd_eeq.cpp


 int EEQModel::get_vrhs(
  const TMolecule &mol,
  const TIVector &realIdx,
  const int &charge,
+  const TMatrix<double> &dist,


Not needed in the function

The dist matrix is needed in the EEQBCModel::get_vrhs

marvinfriede · 2025-08-19T07:37:28Z

src/dftd_eeq.cpp

  const TVector<double> &cn,
-  TVector<double> &Xvec,
-  TVector<double> &dXvec,
+  const TMatrix<double> &dcndr,


Derivative not used in method

marvinfriede · 2025-08-19T07:39:11Z

src/dftd_eeq.cpp

+  const TVector<double> &cn,
+  const TMatrix<double> &dcndr,


They are also not needed. Please double-check the signatures (compare with the Fortran implementation at https://github.com/grimme-lab/multicharge/blob/main/src/multicharge/model.F90)

marvinfriede · 2025-08-19T08:02:39Z

include/dftd_eeq.h

+    TVector<double> qloc;     // local charges
+    TMatrix<double> cmat;     // capacitance matrix
+
+    //


Empty comment

marvinfriede · 2025-08-19T08:03:32Z

include/dftd_eeq.h

      const TVector<double> &cn,
+      const TMatrix<double> &dcndr,


See also below, derivative should not be required.

marvinfriede · 2025-08-19T08:03:56Z

include/dftd_eeq.h

    EEQModel();

    int get_vrhs(
      const TMolecule &mol,
      const TIVector &realIdx,
      const int &charge,
+      const TMatrix<double> &dist,


marvinfriede · 2025-08-19T08:04:29Z

include/dftd_eeq.h

+      const TVector<double> &cn,
+      const TMatrix<double> &dcndr,


Required? Especially, the derivative should not be required.

include/dftd_eeq.h

Thomas3R added 2 commits August 13, 2025 16:51

Implementation of EEQ-BC and unit tests

cc834db

Code cleanup: numerical gradient dqdr

ccfa00c

Thomas3R commented Aug 14, 2025

View reviewed changes

src/dftd_eeq.cpp Outdated Show resolved Hide resolved

Thomas3R commented Aug 14, 2025

View reviewed changes

src/dftd_ncoord.cpp Show resolved Hide resolved

marvinfriede requested a review from Copilot August 15, 2025 11:20

Copilot AI reviewed Aug 15, 2025

View reviewed changes

marvinfriede reviewed Aug 15, 2025

View reviewed changes

src/dftd_ncoord.cpp Outdated Show resolved Hide resolved

marvinfriede reviewed Aug 15, 2025

View reviewed changes

Code cleanup and debugging. Remove cn and dcndr as member variables. …

9f53702

…Adjust EEQ-BC functions for ghost atoms.

marvinfriede requested changes Aug 19, 2025

View reviewed changes

	if (cn_max != -1.0) { // cn_max = -1.0 for EEQ-BC for cn and qloc; using NCoordErf and NCoordErfEN
	if (cn_max < 0.0) { // cn_max = -1.0 for EEQ-BC for cn and qloc; using NCoordErf and NCoordErfEN

		//TVector<double> cn; // coordination number
		//TMatrix<double> dcndr; // derivative of the coordination number

Implementation of EEQ-BC and unit tests #45

Are you sure you want to change the base?

Implementation of EEQ-BC and unit tests #45

Uh oh!

Conversation

Thomas3R commented Aug 13, 2025

Charge Model and Coordination Number Refactoring

Addition of New Charge Model

Coordination Number Class Hierarchy Enhancement

API and Interface Improvements

Memory Management Fix

Uh oh!

Uh oh!

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull Request Overview

Reviewed Changes

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Thomas3R Aug 18, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Thomas3R Aug 18, 2025 •

edited

Loading