diff --git a/include/dftd_cblas.h b/include/dftd_cblas.h
index 7dadc78..75f6fb9 100644
--- a/include/dftd_cblas.h
+++ b/include/dftd_cblas.h
@@ -74,7 +74,7 @@ inline int BLAS_Add_Mat_x_Vec(
   };
 
   return EXIT_FAILURE;
-};
+}
 
 /**
  * @brief General matrix-matrix multiplication (`C = alpha * A * B + C`).
@@ -97,8 +97,9 @@ inline int BLAS_Add_Mat_x_Mat(
 ) {
   // check for size 0 matrices
   if (A.cols == 0 || A.rows == 0 || B.cols == 0 || B.rows == 0 || C.cols == 0 ||
-      C.rows == 0)
-    exit(EXIT_FAILURE);
+      C.rows == 0) {
+      exit(EXIT_FAILURE);
+  };
 
   // check for transpositions
   if (!TransposeA) {
@@ -197,7 +198,7 @@ inline int BLAS_Add_Mat_x_Mat(
     }; // B transposed
   };
   return EXIT_SUCCESS;
-};
+}
 
 /**
  * @brief Compute inverse of a matrix using LU decomposition.
@@ -234,6 +235,51 @@ inline int BLAS_InvertMatrix(TMatrix<double> &a) {
   if (info != 0) { return EXIT_FAILURE; }
 
   return EXIT_SUCCESS;
-};
+}
+
+/**
+ * @brief Solve a symmetric linear system A * X = B for X.
+ *
+ * This routine factorizes a symmetric matrix A using Bunch-Kaufman factorization
+ * and solves for the right-hand side vector B. The matrix A is overwritten
+ * by its factorization. The solution overwrites B.
+ *
+ * @param A Symmetric matrix of size (m x m). Overwritten by the factorization.
+ * @param B Right-hand side vector of size m. Overwritten by the solution.
+ * @return int Returns EXIT_SUCCESS (0) on success, EXIT_FAILURE (1) on error.
+ */
+inline int BLAS_SolveSymmetric(
+  TMatrix<double> &A,   // symmetric matrix
+  TVector<double> &B    // RHS vector (becomes solution)
+) {
+  const lapack_int m    = A.rows;
+  const lapack_int nrhs = 1;
+
+  if (A.cols != m || B.N != m) {
+    fprintf(stderr, "BLAS_SolveSymmetric error: dimension mismatch\n");
+    return EXIT_FAILURE;
+  }
+
+  lapack_int info;
+  lapack_int *ipiv = new lapack_int[A.rows];
+
+  // Factorization
+  info = LAPACKE_dsytrf(LAPACK_ROW_MAJOR, 'L', m, A.p, m, ipiv);
+  if (info != 0) {
+    delete[] ipiv;
+    fprintf(stderr, "dsytrf failed: info=%d\n", (int)info);
+    return EXIT_FAILURE;
+  }
+
+  // Solve for all RHS columns
+  info = LAPACKE_dsytrs(LAPACK_ROW_MAJOR, 'L', m, nrhs, A.p, m, ipiv, B.p, nrhs);
+  delete[] ipiv;
+
+  if (info != 0) {
+    fprintf(stderr, "dsytrs failed: info=%d\n", (int)info);
+    return EXIT_FAILURE;
+  }
+  return EXIT_SUCCESS;
+}
 
 } // namespace dftd4
diff --git a/include/dftd_eeq.h b/include/dftd_eeq.h
index b730501..b88dd93 100644
--- a/include/dftd_eeq.h
+++ b/include/dftd_eeq.h
@@ -134,10 +134,15 @@ class ChargeModel
   * @param charge     Target total charge of the molecule.
   * @param dist       Interatomic distance matrix.
   * @param cn         Coordination numbers for each atom.
-  * @param Xvec       Output vector representing the effective driving term for the
+  * @param dcndr      Derivative of the coordination number w.r.t. atom positions
+  * @param xvec       Output vector representing the effective driving term for the
   *                   EEQ system (size: number of atoms).
-  * @param dXvec      Output vector of derivatives of Xvec w.r.t. atomic positions.
+  * @param dxvecdr    Output vector of derivatives of xvec w.r.t. atom positions.
   *                   Only filled if `lgrad` is true.
+  * @param qloc       Local charge (EEQ-BC specific).
+  * @param dqlocdr    Derivative of the local charge (EEQ-BC specific).
+  * @param cmat       Capacitance matrix.
+  * @param dcmatdr    Derivative of the capacitance matrix.
   * @param lgrad      If true, compute derivatives of Xvec.
   *
   * @returns EXIT_SUCCESS (0) on success, or a nonzero error code if the computation fails.
@@ -148,8 +153,13 @@ class ChargeModel
     const int &charge,
     const TMatrix<double> &dist,
     const TVector<double> &cn,
-    TVector<double> &Xvec,
-    TVector<double> &dXvec,
+    const TMatrix<double> &dcndr,
+    TVector<double> &xvec,
+    TMatrix<double> &dxvecdr,
+    TVector<double> &qloc,
+    TMatrix<double> &dqlocdr,
+    TMatrix<double> &cmat,
+    TMatrix<double> &dcmatdr,
     bool lgrad
   ) = 0;
 
@@ -163,6 +173,8 @@ class ChargeModel
   * @param realIdx Mapping of real atom indices (excludes dummy atoms).
   * @param dist    Interatomic distance matrix.
   * @param cn      Coordination numbers for each atom.
+  * @param qloc    Local charge (EEQ-BC specific)
+  * @param cmat    Capacitance matrix (EEQ-BC specific)
   * @param Amat    Output Coulomb matrix (size: number of atoms + 1 for charge constraint).
   * 
   * @returns EXIT_SUCCESS (0) on success, or a nonzero error code if the computation fails.
@@ -172,6 +184,8 @@ class ChargeModel
     const TIVector &realIdx,
     const TMatrix<double> &dist,
     const TVector<double> &cn,
+    const TVector<double> &qloc,
+    const TMatrix<double> &cmat,
     TMatrix<double> &Amat
   ) const = 0;
   
@@ -185,8 +199,14 @@ class ChargeModel
   * @param dist    Interatomic distance matrix.
   * @param q       Atomic charges computed from the EEQ model.
   * @param Amat    Output coulomb matrix.
-  * @param dAmat   Output derivative matrix with respect to Cartesian coordinates.
+  * @param dAmatdr   Output derivative matrix with respect to Cartesian coordinates.
   * @param atrace  Output trace contributions for the derivative.
+  * @param cn         Coordination numbers for each atom.
+  * @param dcndr      Derivative of the coordination number w.r.t. atom positions
+  * @param qloc       Local charge (EEQ-BC specific).
+  * @param dqlocdr    Derivative of the local charge (EEQ-BC specific).
+  * @param cmat       Capacitance matrix.
+  * @param dcmatdr    Derivative of the capacitance matrix.
   * 
   * @returns EXIT_SUCCESS (0) on success, or a nonzero error code if the computation fails.
   */
@@ -196,8 +216,14 @@ class ChargeModel
     const TMatrix<double> &dist,
     const TVector<double> &q,
     const TMatrix<double> &Amat,
-    TMatrix<double> &dAmat,
-    TMatrix<double> &atrace
+    TMatrix<double> &dAmatdr,
+    TMatrix<double> &atrace,
+    const TVector<double> &cn,
+    const TMatrix<double> &dcndr,
+    const TVector<double> &qloc,
+    const TMatrix<double> &dqlocdr,
+    const TMatrix<double> &cmat,
+    const TMatrix<double> &dcmatdr
   ) const = 0;
 
   /**
@@ -240,7 +266,7 @@ class EEQModel : public ChargeModel {
     // initializes its coordination number function.
     EEQModel();
   
-    // Computes the right-hand side vector (Xvec) for the EEQ charge equations.
+    // Computes the right-hand side vector (xvec) for the EEQ charge equations.
     // Optionally computes its derivatives with respect to atomic positions.
     int get_vrhs(
       const TMolecule &mol,
@@ -248,8 +274,13 @@ class EEQModel : public ChargeModel {
       const int &charge,
       const TMatrix<double> &dist,
       const TVector<double> &cn,
-      TVector<double> &Xvec,
-      TVector<double> &dXvec,
+      const TMatrix<double> &dcndr,
+      TVector<double> &xvec,
+      TMatrix<double> &dxvecdr,
+      TVector<double> &qloc,
+      TMatrix<double> &dqlocdr,
+      TMatrix<double> &cmat,
+      TMatrix<double> &dcmatdr,
       bool lgrad
     ) override;
   
@@ -259,6 +290,8 @@ class EEQModel : public ChargeModel {
       const TIVector &realIdx,
       const TMatrix<double> &dist,
       const TVector<double> &cn,
+      const TVector<double> &qloc,
+      const TMatrix<double> &cmat,
       TMatrix<double> &Amat
     ) const override;
     
@@ -269,8 +302,14 @@ class EEQModel : public ChargeModel {
       const TMatrix<double> &dist,
       const TVector<double> &q,
       const TMatrix<double> &Amat,
-      TMatrix<double> &dAmat,
-      TMatrix<double> &atrace
+      TMatrix<double> &dAmatdr,
+      TMatrix<double> &atrace,
+      const TVector<double> &cn,
+      const TMatrix<double> &dcndr,
+      const TVector<double> &qloc,
+      const TMatrix<double> &dqlocdr,
+      const TMatrix<double> &cmat,
+      const TMatrix<double> &dcmatdr
     ) const override;
 
     // Calculate the coordination number, forwarding to get_ncoord
@@ -324,8 +363,13 @@ class EEQBCModel : public ChargeModel {
       const int &charge,
       const TMatrix<double> &dist,
       const TVector<double> &cn,
-      TVector<double> &Xvec,
-      TVector<double> &dXvec,
+      const TMatrix<double> &dcndr,
+      TVector<double> &xvec,
+      TMatrix<double> &dxvecdr,
+      TVector<double> &qloc,
+      TMatrix<double> &dqlocdr,
+      TMatrix<double> &cmat,
+      TMatrix<double> &dcmatdr,
       bool lgrad
     ) override;
   
@@ -335,6 +379,8 @@ class EEQBCModel : public ChargeModel {
       const TIVector &realIdx,
       const TMatrix<double> &dist,
       const TVector<double> &cn,
+      const TVector<double> &qloc,
+      const TMatrix<double> &cmat,
       TMatrix<double> &Amat
     ) const override;
     
@@ -345,8 +391,14 @@ class EEQBCModel : public ChargeModel {
       const TMatrix<double> &dist,
       const TVector<double> &q,
       const TMatrix<double> &Amat,
-      TMatrix<double> &dAmat,
-      TMatrix<double> &atrace
+      TMatrix<double> &dAmatdr,
+      TMatrix<double> &atrace,
+      const TVector<double> &cn,
+      const TMatrix<double> &dcndr,
+      const TVector<double> &qloc,
+      const TMatrix<double> &dqlocdr,
+      const TMatrix<double> &cmat,
+      const TMatrix<double> &dcmatdr
     ) const override;
 
     // Calculate the coordination number, forwarding to get_ncoord
@@ -361,11 +413,13 @@ class EEQBCModel : public ChargeModel {
 
     // Get purely geometry-dependent local charges
     int get_qloc(
-      const TMolecule&,
-      const TIVector &,
-      const TMatrix<double>&,
-      const double,
-      TVector<double>&
+      const TMolecule &mol,
+      const TIVector &realIdx,
+      const TMatrix<double>& dist,
+      const double q_tot,
+      TVector<double> &qloc,
+      TMatrix<double> &dqlocdr,
+      bool lgrad
     );
 
     // Get the capacitance for bond between atoms i and j
@@ -376,31 +430,57 @@ class EEQBCModel : public ChargeModel {
       double &c_ij
     ) const;
 
+    // Get the capacitance for bond between atoms i and j
+    int get_dcpair(
+      double dist_ij,  // distance between atom j to atom i
+      double rvdw_ijat,  // pairwise van der Waals radii for ij atom types
+      double cap_ij,  // product of bond capacitances for atom types of i and j
+      TVector<double> &vec, // Vector from j to i
+      TVector<double> &dcdr_ij  // Out: Capacitance for bond ij
+    ) const;
+
     // Get the capacitance matrix
     int get_cmat(
-      const TMolecule&,
+      const TMolecule& mol,
       const TIVector &realIdx,
-      const TMatrix<double>&,
-      TMatrix<double>&
+      const TMatrix<double>& dist,
+      TMatrix<double>& cmat
+    );
+
+    // Get the capacitance matrix
+    int get_dcmatdr(
+      const TMolecule &mol,
+      const TIVector &realIdx,
+      const TMatrix<double> &dist,
+      TMatrix<double> &cmat
     );
 
     // Get the right-hand side (electronegativity) of the set of linear equations
     int get_xvec(
-      const TMolecule&,
-      const TIVector&,
-      const TMatrix<double>&,
+      const TMolecule &mol,
+      const TIVector &realIdx,
+      const TMatrix<double> &dist,
       const TVector<double> &cn,
-      TMatrix<double>&,
-      int,
-      TVector<double>&
+      TMatrix<double> &cmat,
+      int charge,
+      TVector<double> &qloc,
+      TVector<double> &xvec
     );
 
-    // numerical gradient of partial charges w.r.t. atom positions
-    int num_grad_dqdr(
-      TMolecule&,  // molecular geometry
-      const TIVector&,
-      int,
-      TMatrix<double>& // numerical gradient
+    // Get the derivative of the right-hand side of the set of linear equations
+    int get_xvec_derivs(
+      const TMolecule &mol,
+      const TIVector &realIdx,
+      const TMatrix<double> &dist,
+      const TVector<double> &cn,
+      const TMatrix<double> &dcndr,
+      TMatrix<double> &cmat,
+      int charge,
+      TVector<double> &xvec,
+      TMatrix<double> &dxvecdr,
+      TVector<double> &qloc,
+      TMatrix<double> &dqlocdr,
+      TMatrix<double> &dcmatdr
     );
 
 };
diff --git a/include/dftd_parameters.h b/include/dftd_parameters.h
index 112dc2b..ce93fda 100644
--- a/include/dftd_parameters.h
+++ b/include/dftd_parameters.h
@@ -28,7 +28,8 @@
 
 namespace dftd4 {
 
-constexpr int MAXELEMENT = 104; // 103 + dummy
+// TODO keep in sync with qcelems.h???
+constexpr int MAXELEMENT = 118+2; //dummy + H..Og + point charge
 
 // clang-format off
 /**
@@ -266,6 +267,7 @@ static const int refn[MAXELEMENT]{
   3, 4, 4, 2, 2, 3, 5, 4, 3, 2, 1, 3, 4, 3, 4, 4, 4, 3, 3, 4, 3, 2, 2, 4, 5, 4,
   3, 2, 1, 3, 4, 3, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 4, 3, 3, 3, 5,
   3, 2, 2, 4, 5, 4, 3, 2, 1, 2, 3, 7, 5, 6, 7, 7, 7, 7, 7, 5, 7, 7, 7, 5, 7, 7,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
 };
 
 // 1
@@ -6146,7 +6148,9 @@ static const double *refcovcn[MAXELEMENT]{
   refcovcn84,  refcovcn85,  refcovcn86, refcovcn87, refcovcn88,  refcovcn89,
   refcovcn90,  refcovcn91,  refcovcn92, refcovcn93, refcovcn94,  refcovcn95,
   refcovcn96,  refcovcn97,  refcovcn98, refcovcn99, refcovcn100, refcovcn101,
-  refcovcn102, refcovcn103,
+  refcovcn102, refcovcn103, nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr
 };
 static const double *refq_eeq[MAXELEMENT]{
   nullptr,     refq_eeq1,   refq_eeq2,  refq_eeq3,  refq_eeq4,   refq_eeq5,
@@ -6166,7 +6170,9 @@ static const double *refq_eeq[MAXELEMENT]{
   refq_eeq84,  refq_eeq85,  refq_eeq86, refq_eeq87, refq_eeq88,  refq_eeq89,
   refq_eeq90,  refq_eeq91,  refq_eeq92, refq_eeq93, refq_eeq94,  refq_eeq95,
   refq_eeq96,  refq_eeq97,  refq_eeq98, refq_eeq99, refq_eeq100, refq_eeq101,
-  refq_eeq102, refq_eeq103,
+  refq_eeq102, refq_eeq103, nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr
 };
 static const double *refsq[MAXELEMENT]{
   nullptr, refsq1,  refsq2,  refsq3,  refsq4,   refsq5,   refsq6,   refsq7,
@@ -6182,7 +6188,11 @@ static const double *refsq[MAXELEMENT]{
   refsq80, refsq81, refsq82, refsq83, refsq84,  refsq85,  refsq86,  refsq87,
   refsq88, refsq89, refsq90, refsq91, refsq92,  refsq93,  refsq94,  refsq95,
   refsq96, refsq97, refsq98, refsq99, refsq100, refsq101, refsq102, refsq103,
+                            nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr
 };
+
 static const double *refalpha[MAXELEMENT]{
   nullptr,     refalpha1,   refalpha2,  refalpha3,  refalpha4,   refalpha5,
   refalpha6,   refalpha7,   refalpha8,  refalpha9,  refalpha10,  refalpha11,
@@ -6201,7 +6211,9 @@ static const double *refalpha[MAXELEMENT]{
   refalpha84,  refalpha85,  refalpha86, refalpha87, refalpha88,  refalpha89,
   refalpha90,  refalpha91,  refalpha92, refalpha93, refalpha94,  refalpha95,
   refalpha96,  refalpha97,  refalpha98, refalpha99, refalpha100, refalpha101,
-  refalpha102, refalpha103,
+  refalpha102, refalpha103, nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr
 };
 static const double *refascale[MAXELEMENT]{
   nullptr,      refascale1,   refascale2,   refascale3,   refascale4,
@@ -6225,6 +6237,9 @@ static const double *refascale[MAXELEMENT]{
   refascale90,  refascale91,  refascale92,  refascale93,  refascale94,
   refascale95,  refascale96,  refascale97,  refascale98,  refascale99,
   refascale100, refascale101, refascale102, refascale103,
+                            nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr
 };
 static const double *refscount[MAXELEMENT]{
   nullptr,      refscount1,   refscount2,   refscount3,   refscount4,
@@ -6248,6 +6263,9 @@ static const double *refscount[MAXELEMENT]{
   refscount90,  refscount91,  refscount92,  refscount93,  refscount94,
   refscount95,  refscount96,  refscount97,  refscount98,  refscount99,
   refscount100, refscount101, refscount102, refscount103,
+                            nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr
 };
 static const int *refsys[MAXELEMENT]{
   nullptr,  refsys1,  refsys2,   refsys3,   refsys4,   refsys5,   refsys6,
@@ -6265,6 +6283,9 @@ static const int *refsys[MAXELEMENT]{
   refsys84, refsys85, refsys86,  refsys87,  refsys88,  refsys89,  refsys90,
   refsys91, refsys92, refsys93,  refsys94,  refsys95,  refsys96,  refsys97,
   refsys98, refsys99, refsys100, refsys101, refsys102, refsys103,
+                            nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr
 };
 static const int *refc[MAXELEMENT]{
   nullptr, refc1,   refc2,   refc3,   refc4,   refc5,  refc6,  refc7,  refc8,
@@ -6279,6 +6300,9 @@ static const int *refc[MAXELEMENT]{
   refc81,  refc82,  refc83,  refc84,  refc85,  refc86, refc87, refc88, refc89,
   refc90,  refc91,  refc92,  refc93,  refc94,  refc95, refc96, refc97, refc98,
   refc99,  refc100, refc101, refc102, refc103,
+                            nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr,
+  nullptr,     nullptr,     nullptr,    nullptr,    nullptr,     nullptr
 };
 
 // sec. 1
diff --git a/scripts/make-qcdftd4.sh b/scripts/make-qcdftd4.sh
index 23db57b..0bb145c 100755
--- a/scripts/make-qcdftd4.sh
+++ b/scripts/make-qcdftd4.sh
@@ -151,7 +151,24 @@ cat >> qcdftd4.h << 'EOT'
 // https://doi.org/10.1063/1.5090222
 // Bernardo de Souza, 14/09/2023
  -------------------------------------------------------------------------------------- */
-int CalcEEQCharges(TRMatrix &XYZ, TIVector &ATNO, int NAtoms, int totalcharge, TRVector &q, bool printerror=true);
+int CalcEEQCharges(TRMatrix &XYZ, TIVector &ATNO, int NAtoms, int totalcharge, TRVector &q,
+                   bool printerror=true, bool allowallatoms=false);
+
+/* --------------------------------------------------------------------------------------
+// Calculates the EEQ-BC charges according to the paper
+// https://doi.org/10.1063/5.0268978
+// Thomas Rose, 14/10/2025
+ -------------------------------------------------------------------------------------- */
+int CalcEEQBCCharges(TRMatrix &XYZ, TIVector &ATNO, int NAtoms, int totalcharge, TRVector &q
+                   , bool printerror=true, bool allowalaltoms=false);
+
+/* --------------------------------------------------------------------------------------
+// Calculates the coordination number and the EEQ charges as done in the D4 paper
+// https://doi.org/10.1063/1.5090222
+// Miquel Garcia-Ratés, 03/2024
+ -------------------------------------------------------------------------------------- */
+int CalcCNEEQ(TOrcaInfo &DAT, TGeomInput &GEO, TRVector &CN, TRVector &Q, TRMatrix &GRADQ,
+              bool lgrad, bool printerror=true);
 
 #endif // __QCDFTD4_H
 
@@ -196,18 +213,55 @@ cat > qcdftd4.cpp << 'EOT'
 // https://doi.org/10.1063/1.5090222
 // Bernardo de Souza, 14/09/2023
  -------------------------------------------------------------------------------------- */
-int CalcEEQCharges(TRMatrix &XYZ, TIVector &ATNO, int NAtoms, int totalcharge, TRVector &q, bool printerror){
+int CalcEEQCharges(TRMatrix &XYZ, TIVector &ATNO, int NAtoms,
+                   int totalcharge, TRVector &q,
+                   bool printerror, bool allowallatoms){
 
   // initialize the charges to zero
   q.Init();
 
+  TGeomInput molecule;
+  molecule.NAtoms=NAtoms;
+  molecule.CC.CopyMat(XYZ);
+  molecule.ATNO.CopyVec(ATNO);
+  TIVector Index(NAtoms);
+  for (int i=0;i<NAtoms;i++) Index(i)=i;
+  TRMatrix dist;
+  dist.NewMatrix(NAtoms, NAtoms);
+  multicharge::EEQModel chrg_model;
   // check atomic numbers to guarantee we have all parameters
   // ghost atoms (ATNO=0) will have no charge
-  for (int i=0;i<NAtoms;i++)
-    if (ATNO(i)>86){
-      if (printerror) printMessage("Atomic number %d detected. EEQ charges can not be calculated.\n",ATNO(i));
-      return 1;
+  for (int i=0;i<NAtoms;i++){
+    if (molecule.ATNO(i)>86){
+      if (!allowallatoms){
+        if (printerror) printMessage("Atomic number %d detected. EEQ charges can not be calculated.\n",ATNO(i));
+        return 1;
+      }
+      // reduce to the previous row
+      else{
+        molecule.ATNO(i) -= 32;
+      }
     }
+  }
+
+  dftd4::calc_distances(molecule, Index, dist);
+
+  TRMatrix dqdr;
+
+  return chrg_model.get_charges(molecule, Index, dist, totalcharge, 25, q, dqdr, false);
+}
+
+/* --------------------------------------------------------------------------------------
+// Calculates the EEQ-BC charges according to the paper
+// https://doi.org/10.1063/5.0268978
+// Thomas Rose, 14/10/2025
+ -------------------------------------------------------------------------------------- */
+int CalcEEQBCCharges(TRMatrix &XYZ, TIVector &ATNO, int NAtoms,
+                   int totalcharge, TRVector &q,
+                   bool printerror, bool allowallatoms){
+
+  // initialize the charges to zero
+  q.Init();
 
   TGeomInput molecule;
   molecule.NAtoms=NAtoms;
@@ -217,11 +271,78 @@ int CalcEEQCharges(TRMatrix &XYZ, TIVector &ATNO, int NAtoms, int totalcharge, T
   for (int i=0;i<NAtoms;i++) Index(i)=i;
   TRMatrix dist;
   dist.NewMatrix(NAtoms, NAtoms);
+  multicharge::EEQBCModel chrg_model;
+  // check atomic numbers to guarantee we have all parameters
+  // ghost atoms (ATNO=0) will have no charge
+  for (int i=0;i<NAtoms;i++){
+    if (molecule.ATNO(i)>103){
+      if (!allowallatoms){
+        if (printerror) printMessage("Atomic number %d detected. EEQ-BC charges can not be calculated.\n",ATNO(i));
+        return 1;
+      }
+      // reduce to the previous row
+      else{
+        molecule.ATNO(i) -= 32;
+      }
+    }
+  }
+
+
   dftd4::calc_distances(molecule, Index, dist);
 
   TRMatrix dqdr;
 
-  return dftd4::get_charges(molecule, Index, dist, totalcharge, 25, q, dqdr, false);
+  return chrg_model.get_charges(molecule, Index, dist, totalcharge, 25, q, dqdr, false);
+}
+
+/* --------------------------------------------------------------------------------------
+// Calculates the coordination number and the EEQ charges as done in the D4 paper
+// https://doi.org/10.1063/1.5090222
+// Miquel Garcia-Ratés, 03/2024
+ -------------------------------------------------------------------------------------- */
+int CalcCNEEQ(TOrcaInfo &DAT, TGeomInput &GEO, TRVector &CN, TRVector &Q, TRMatrix &GRADQ,
+              bool lgrad, bool printerror){
+
+  // setup variables
+  int info{0};
+
+  int NAtoms = GEO.NAtoms;
+
+  TGeomInput molecule;
+  molecule.NAtoms = NAtoms;
+  molecule.CC.CopyMat(GEO.CC);
+  molecule.ATNO.CopyVec(GEO.ATNO);
+  TIVector Index(NAtoms);
+  for (int i=0;i<NAtoms;i++) Index(i) = i;
+
+  // Calculate distances
+  TRMatrix dist;
+  dist.NewMatrix(NAtoms, NAtoms);
+  dftd4::calc_distances(molecule, Index, dist);
+  multicharge::EEQModel eeq_model;
+  const double kcn = 7.5, norm_exp = 1.0, cutoff = 30.0;
+  dftd4::NCoordErfD4 ncoord_erf_d4 = dftd4::NCoordErfD4(kcn, norm_exp, cutoff);
+
+  TMatrix<double> dcndr;    // derivative of D4-CN
+
+  int nat = Index.Max() + 1;
+
+  Q.NewVector(NAtoms);
+  if (lgrad) {
+    dcndr.NewMatrix(3*NAtoms, nat);
+    GRADQ.NewMatrix(3*NAtoms, nat);
+  }
+
+  // calculate partial charges from EEQ model
+  info = eeq_model.get_charges(molecule, Index, dist, DAT.Charge, 25.0, Q, GRADQ, lgrad);
+  if (info != EXIT_SUCCESS) return info;
+
+  // get the D4 coordination number
+  info = ncoord_erf_d4.get_ncoord(molecule, Index, dist, CN, dcndr, lgrad);
+  if (info != EXIT_SUCCESS) return info;
+
+  return 0;
+
 }
 
 namespace dftd4 {
diff --git a/scripts/make-qceeq.sh b/scripts/make-qceeq.sh
index 4a6cf63..84f89ef 100755
--- a/scripts/make-qceeq.sh
+++ b/scripts/make-qceeq.sh
@@ -4,8 +4,8 @@ INCLUDE="../include"
 SRC="../src"
 
 dosed(){
-	sed -n '/namespace dftd4 {/,/} \/\/ namespace dftd4/p' $1 | \
-	 sed '/namespace dftd4 {/d; /} \/\/ namespace dftd4/d' | \
+	sed -n '/namespace multicharge {/,/} \/\/ namespace multicharge/p' $1 | \
+	 sed '/namespace multicharge {/d; /} \/\/ namespace multicharge/d' | \
 	 awk '/./,EOF' > $2
 }
 
@@ -25,13 +25,13 @@ cat > qceeq.h << 'EOT'
  * Marvin Friede <friede@thch.uni-bonn.de> (MF121223)
  */
 
-#ifndef QCEEQ_H
-#define QCEEQ_H
+#pragma once
 
 #include "qcinpdat.h" // TGeomInput class
 #include "qcmat1.h"   // TRVector and TRMatrix class
+#include "qcncoord.h"
 
-namespace dftd4 {
+namespace multicharge {
 
 EOT
 
@@ -39,9 +39,7 @@ cat eeq.txt >> qceeq.h
 rm eeq.txt
 
 cat >> qceeq.h << 'EOT'
-} // namespace dftd4
-
-#endif // QCEEQ_H
+} // namespace multicharge
 EOT
 
 ######
@@ -73,16 +71,15 @@ cat > qceeq.cpp << 'EOT'
 // always include self
 #include "qceeq.h"
 
-// wrap everything in the dftd namespace to keep it nicely confined
-namespace dftd4 {
-
+// wrap all charge models in the multicharge namespace
+namespace multicharge {
 EOT
 
 cat eeq.txt >> qceeq.cpp
 rm eeq.txt
 
 cat >> qceeq.cpp << 'EOT'
-} // namespace dftd4
+} // namespace multicharge
 EOT
 
 ##########################################################################
diff --git a/scripts/make-qcncoord.sh b/scripts/make-qcncoord.sh
index c0a18f3..85db1c3 100755
--- a/scripts/make-qcncoord.sh
+++ b/scripts/make-qcncoord.sh
@@ -25,11 +25,11 @@ cat > qcncoord.h << 'EOT'
  * Marvin Friede <friede@thch.uni-bonn.de> (MF121223)
  */
 
-#ifndef QCNCOORD_H
-#define QCNCOORD_H
+#pragma once
 
 #include "qcinpdat.h" // TGeomInput class
 #include "qcmat1.h"   // TRVector and TRMatrix class
+#include "qc_multicharge_param.h"  // Chargemodel Parameters
 
 namespace dftd4 {
 
@@ -40,8 +40,6 @@ rm ncoord.txt
 
 cat >> qcncoord.h << 'EOT'
 } // namespace dftd4
-
-#endif // QCNCOORD_H
 EOT
 
 ######
diff --git a/src/dftd_eeq.cpp b/src/dftd_eeq.cpp
index 62592d1..b4277c7 100644
--- a/src/dftd_eeq.cpp
+++ b/src/dftd_eeq.cpp
@@ -54,7 +54,7 @@ int ChargeModel::get_charges(
   bool lgrad
 ) {
   TIVector realIdx;
-  initializeRealIdx(mol.NAtoms, realIdx);
+  dftd4::initializeRealIdx(mol.NAtoms, realIdx);
 
   return get_charges(mol, realIdx, dist, charge, cutoff, q, dqdr, lgrad);
 };
@@ -99,51 +99,76 @@ int ChargeModel::eeq_chrgeq(
   bool lgrad /*= false*/,
   bool lverbose /*= false*/
 ) {
-  double qtotal = 0.0;
   int info{0};
-  const int n = realIdx.Max() + 1;
-  int m = n + 1;
+  const int n = realIdx.Max() + 1;  // Number of atoms
+  int m = n + 1;  // Number of atoms plus one for constraint
 
-  TMatrix<double> Amat; // Coulomb matrix
-  TVector<double> xvec; // x (chi) vector
+  TMatrix<double> Amat;  // Coulomb matrix
+  TVector<double> xvec;  // x (chi) vector, practically the right-hand side of the set of linear equations
+  TMatrix<double> dxvecdr;  // Derivative of the x vector w.r.t. atom positions
   Amat.NewMat(m, m);
   xvec.NewVec(m);
 
-  TVector<double> dxdcn; // Derivative of chi vector w.r.t. CN
-  if (lgrad) dxdcn.NewVec(m);
+  // EEQ-BC specific variables:
+  TVector<double> qloc;     // Local charge (allocated in get_qloc)
+  TMatrix<double> dqlocdr;  // Derivative of local charge w.r.t. atom positions (allocated in get_ncoord)
+  TMatrix<double> cmat;     // Capacitance matrix (allocated in get_cmat)
+  TMatrix<double> dcmatdr;  // Derivative of the capacitance matrix w.r.t. atom positions (allocated in get_dcmatdr)
 
-  info = get_vrhs(mol, realIdx, charge, dist, cn, xvec, dxdcn, lgrad);
+  info = get_vrhs(mol, realIdx, charge, dist, cn, dcndr, xvec, dxvecdr, qloc, dqlocdr, cmat, dcmatdr, lgrad);
   if (info != EXIT_SUCCESS) return info;
 
-  info = get_amat_0d(mol, realIdx, dist, cn, Amat);
+  info = get_amat_0d(mol, realIdx, dist, cn, qloc, cmat, Amat);
   if (info != EXIT_SUCCESS) return info;
 
-  TVector<double> vrhs;
+  TVector<double> vrhs;  // right-hand side (RHS) of linear equations; finally holds the solution of the linear equations
   vrhs.NewVec(m);
 
-  TMatrix<double> Ainv;
-  Ainv.NewMat(m, m);
-  Ainv.CopyMat(Amat);
+  TMatrix<double> Ainv;  // Inverse of the Coulomb matrix for dqdr
+  // Solve the set of linear equations with the Inverse Ainv
+  //   which we need later for dqdr as requested by lgrad=true
+  if (lgrad) {
+    Ainv.NewMat(m, m);
+    Ainv.CopyMat(Amat);
 
-  // solve: A Q = X (Eq.4) -> Q = Ainv X
-  info = BLAS_InvertMatrix(Ainv);
-  if (info != EXIT_SUCCESS) return info;
+    // solve: A Q = X
+    info = BLAS_InvertMatrix(Ainv);
+    if (info != EXIT_SUCCESS) return info;
 
-  // no return in ORCA
-  BLAS_Add_Mat_x_Vec(vrhs, Ainv, xvec, false, 1.0);
-  // if (info != EXIT_SUCCESS) return info;
+    // no return in ORCA
+    BLAS_Add_Mat_x_Vec(vrhs, Ainv, xvec, false, 1.0);
+    // if (info != EXIT_SUCCESS) return info;
+  }
+  // Solve the set of linear equations directly (faster since we do not need Ainv)
+  else {
+    TVector<double> rhs;
+    rhs.NewVec(m);
+    rhs.CopyVec(xvec);
+
+    // Now solve (A * q = x) for symmetric Coulomb matrix Amat
+    // Note that only the lower triangle of Amat is updated
+    info = BLAS_SolveSymmetric(Amat, rhs);
+    if (info != EXIT_SUCCESS) return info;
+
+    // Extract the solution
+    for (int i = 0, ii = 0; i < mol.NAtoms; i++) {
+        ii = realIdx(i);
+        if (ii < 0) continue;
+        vrhs(ii) = rhs(ii);
+    }
+  }
 
-  // remove charge constraint (make vector smaller by one) and total charge
-  qtotal = 0.0;
+  // Remove charge constraint (make vector smaller by one) and total charge
+  double qtotal = 0.0;  // total charge for checking charge constraint
   for (int iat = 0, ii = 0; iat != mol.NAtoms; iat++) {
     ii = realIdx(iat);
     if (ii < 0) continue;
 
-    qvec(ii) = vrhs(ii);
-    qtotal += qvec(ii);
+    qvec(ii) = vrhs(ii);  // Assign partial charges to the correct variable
+    qtotal += qvec(ii);  // Total charge of the system
   }
 
-  // check total charge and additional printout
+  // Check total charge and additional printout
   if (fabs(qtotal - charge) > 1.0e-8) {
     printf(
       "DFT-D4: EEQ charge constraint error: %14.8f vs. %14d\n", qtotal, charge
@@ -168,35 +193,38 @@ int ChargeModel::eeq_chrgeq(
 
   // Gradient (note that the corresponding gradient flag in Fortran is `cpq`)
   if (lgrad) {
-    TMatrix<double> dAmat;
-    dAmat.NewMat(3 * n, m);
+    TMatrix<double> dAmatdr;
+    dAmatdr.NewMat(3 * n, m);
     TMatrix<double> atrace;
     atrace.NewMat(m, 3);
 
-    info = get_damat_0d(mol, realIdx, dist, vrhs, Amat, dAmat, atrace);
+    // Calculate the derivative of the Coulomb matrix w.r.t. atom positions (dAmatdr)
+    info = get_damat_0d(mol, realIdx, dist, vrhs, Amat, dAmatdr, atrace,
+       cn, dcndr, qloc, dqlocdr, cmat, dcmatdr);
     if (info != EXIT_SUCCESS) return info;
 
     for (int i = 0, ii = 0; i != mol.NAtoms; i++) {
       ii = realIdx(i);
       if (ii < 0) continue;
 
-      dAmat(3 * ii, ii) += atrace(ii, 0);
-      dAmat(3 * ii + 1, ii) += atrace(ii, 1);
-      dAmat(3 * ii + 2, ii) += atrace(ii, 2);
+      dAmatdr(3 * ii, ii) += atrace(ii, 0);
+      dAmatdr(3 * ii + 1, ii) += atrace(ii, 1);
+      dAmatdr(3 * ii + 2, ii) += atrace(ii, 2);
 
       for (int j = 0, jj = 0; j != mol.NAtoms; j++) {
         jj = realIdx(j);
         if (jj < 0) continue;
 
-        dAmat(3 * jj,     ii) -= dcndr(ii, 3 * jj    ) * dxdcn(ii);
-        dAmat(3 * jj + 1, ii) -= dcndr(ii, 3 * jj + 1) * dxdcn(ii);
-        dAmat(3 * jj + 2, ii) -= dcndr(ii, 3 * jj + 2) * dxdcn(ii);
+        // dxvecdr = del C/del R * X + del X / del R * C 
+        // CN contribution
+        dAmatdr(3 * jj,     ii) -= dxvecdr(ii, 3 * jj    );
+        dAmatdr(3 * jj + 1, ii) -= dxvecdr(ii, 3 * jj + 1);
+        dAmatdr(3 * jj + 2, ii) -= dxvecdr(ii, 3 * jj + 2);
       }
     }
 
     // we do not need these gradient-related matrices anymore
     atrace.DelMat();
-    dxdcn.DelVec();
 
     // Ainv with last column removed
     TMatrix<double> A;
@@ -207,17 +235,25 @@ int ChargeModel::eeq_chrgeq(
       }
     }
 
+    // Calculate the charge gradient w.r.t. atom positions (dqdr)
     // no return in ORCA
-    BLAS_Add_Mat_x_Mat(dqdr, dAmat, A, false, false, -1.0);
+    BLAS_Add_Mat_x_Mat(dqdr, dAmatdr, A, false, false, -1.0);
     // if (info != EXIT_SUCCESS) return info;
 
-    dAmat.DelMat();
+    dAmatdr.DelMat();
+    Ainv.DelMat();
+    A.DelMat();
+    dqlocdr.DelMat();
+    dxvecdr.DelMat();
+    dcmatdr.DelMat();
   }
 
   // free all memory
-  Ainv.DelMat();
   Amat.DelMat();
   xvec.DelVec();
+  qloc.DelVec();
+  cmat.DelMat();
+  vrhs.DelVec();
 
   return EXIT_SUCCESS;
 }
@@ -239,15 +275,23 @@ int EEQModel::get_vrhs(
   const int &charge,
   const TMatrix<double> &dist,
   const TVector<double> &cn,
+  const TMatrix<double> &dcndr,
   TVector<double> &xvec,
-  TVector<double> &dxvec,
+  TMatrix<double> &dxvecdr,
+  TVector<double> &qloc,
+  TMatrix<double> &dqlocdr,
+  TMatrix<double> &cmat,
+  TMatrix<double> &dcmatdr,
   bool lgrad
 ) {
+  TVector<double> dxdcn;
   double tmp{0.0};
   int izp;
   const int nat = realIdx.Max() + 1;
 
   if (lgrad) {
+    dxvecdr.NewMat(nat, 3*nat);
+    dxdcn.NewVec(nat);
     for (int i = 0, ii = 0; i != mol.NAtoms; i++) {
       ii = realIdx(i);
       if (ii < 0) continue;
@@ -255,9 +299,16 @@ int EEQModel::get_vrhs(
       izp = mol.ATNO(i);
       tmp = kappa[izp] / std::sqrt(cn(ii) + small);
       xvec(ii) = -xi[izp] + tmp * cn(ii);
-      dxvec(ii) = 0.5 * tmp;
+      dxdcn(ii) = 0.5 * tmp;
+
+      for (int j = 0, jj = 0; j != mol.NAtoms; j++) {
+        jj = realIdx(j);
+        if (jj < 0) continue;
+        dxvecdr(ii, 3*jj  ) = dxdcn(ii) * dcndr(ii, 3 * jj    ); 
+        dxvecdr(ii, 3*jj+1) = dxdcn(ii) * dcndr(ii, 3 * jj + 1); 
+        dxvecdr(ii, 3*jj+2) = dxdcn(ii) * dcndr(ii, 3 * jj + 2); 
+      }
     }
-    dxvec(nat) = 0.0;
   } else {
     for (int i = 0, ii = 0; i != mol.NAtoms; i++) {
       ii = realIdx(i);
@@ -280,6 +331,8 @@ int EEQModel::get_amat_0d(
   const TIVector &realIdx,
   const TMatrix<double> &dist,
   const TVector<double> &cn,
+  const TVector<double> &qloc,
+  const TMatrix<double> &cmat,
   TMatrix<double> &Amat
 ) const {
   double gamij = 0.0;
@@ -323,8 +376,14 @@ int EEQModel::get_damat_0d(
   const TMatrix<double> &dist,
   const TVector<double> &q,
   const TMatrix<double> &Amat,
-  TMatrix<double> &dAmat,
-  TMatrix<double> &atrace
+  TMatrix<double> &dAmatdr,
+  TMatrix<double> &atrace,
+  const TVector<double> &cn,
+  const TMatrix<double> &dcndr,
+  const TVector<double> &qloc,
+  const TMatrix<double> &dqlocdr,
+  const TMatrix<double> &cmat,
+  const TMatrix<double> &dcmatdr
 ) const {
   double alphai, alphaj;
   double rx, ry, rz, r2;
@@ -362,12 +421,12 @@ int EEQModel::get_damat_0d(
       atrace(jj, 1) -= dgy * q(ii);
       atrace(jj, 2) -= dgz * q(ii);
 
-      dAmat(3 * ii,     jj) =  dgx * q(ii);
-      dAmat(3 * ii + 1, jj) =  dgy * q(ii);
-      dAmat(3 * ii + 2, jj) =  dgz * q(ii);
-      dAmat(3 * jj,     ii) = -dgx * q(jj);
-      dAmat(3 * jj + 1, ii) = -dgy * q(jj);
-      dAmat(3 * jj + 2, ii) = -dgz * q(jj);
+      dAmatdr(3 * ii,     jj) =  dgx * q(ii);
+      dAmatdr(3 * ii + 1, jj) =  dgy * q(ii);
+      dAmatdr(3 * ii + 2, jj) =  dgz * q(ii);
+      dAmatdr(3 * jj,     ii) = -dgx * q(jj);
+      dAmatdr(3 * jj + 1, ii) = -dgy * q(jj);
+      dAmatdr(3 * jj + 2, ii) = -dgz * q(jj);
     }
   }
 
@@ -430,26 +489,36 @@ int EEQBCModel::get_vrhs(
   const int &charge,
   const TMatrix<double> &dist,
   const TVector<double> &cn,
+  const TMatrix<double> &dcndr,
   TVector<double> &xvec,
-  TVector<double> &dxvec,
+  TMatrix<double> &dxvecdr,
+  TVector<double> &qloc,
+  TMatrix<double> &dqlocdr,
+  TMatrix<double> &cmat,
+  TMatrix<double> &dcmatdr,
   bool lgrad
 ) {
   int info{0};
   const int n_atoms = realIdx.Max() + 1;
   // calculate the capacitance matrix
-  cmat.NewMatrix(n_atoms + 1, n_atoms + 1);
   info = get_cmat(mol, realIdx, dist, cmat);
   if (info != EXIT_SUCCESS) {
      printf("EEQBCModel::get_vrhs: Failed to calculate the capacitance matrix.");
      return info;
   }
   
-  // calculate the right-hand side
-  info = get_xvec(mol, realIdx, dist, cn, cmat, charge, xvec);
+  if ( ! lgrad) {
+    // calculate the right-hand side (RHS)
+    info = get_xvec(mol, realIdx, dist, cn, cmat, charge, qloc, xvec);
+  } else {
+    // calculate RHS and derivative
+    info = get_xvec_derivs(mol, realIdx, dist, cn, dcndr, cmat, charge, xvec, dxvecdr, qloc, dqlocdr, dcmatdr);
+  }
   if (info != EXIT_SUCCESS) {
      printf("EEQBCModel::get_vrhs: Failed to calculate the right hand side.");
      return info;
   }
+
   return EXIT_SUCCESS;
 };
 
@@ -459,15 +528,17 @@ int EEQBCModel::get_amat_0d(
   const TIVector &realIdx,
   const TMatrix<double> &dist,
   const TVector<double> &cn,
+  const TVector<double> &qloc,
+  const TMatrix<double> &cmat,
   TMatrix<double> &Amat
 ) const {
-  const int nat = mol.NAtoms;
+  const int nat = realIdx.Max() + 1;
   int iat, jat; // atomic numbers
   double norm_cn; // coordination number normalization factor
   double r, radi, radj, gamij2, tmp;
 
 
-  for (int i = 0, ii = 0; i < nat; i++) {
+  for (int i = 0, ii = 0; i < mol.NAtoms; i++) {
     ii = realIdx(i);
     if (ii < 0) continue;
     iat = mol.ATNO(i);
@@ -504,9 +575,107 @@ int EEQBCModel::get_damat_0d(
   const TMatrix<double> &dist,
   const TVector<double> &q,
   const TMatrix<double> &Amat,
-  TMatrix<double> &dAmat,
-  TMatrix<double> &atrace
+  TMatrix<double> &dAmatdr,
+  TMatrix<double> &atrace,
+  const TVector<double> &cn,
+  const TMatrix<double> &dcndr,
+  const TVector<double> &qloc,
+  const TMatrix<double> &dqlocdr,
+  const TMatrix<double> &cmat,
+  const TMatrix<double> &dcmatdr
 ) const {
+
+  int izp, jzp;
+  double norm_cn;
+  double gam;
+  double arg, dtmp1, dtmp2, dtmp3, dtmp4, dtmp5;
+  double radi, dradcn_i, radj, dradcn_j, r2;
+  const int n_atoms = realIdx.Max() + 1;
+  TVector<double> dgamdr;
+  TVector<double> vec, dG;
+  vec.NewVec(3);
+  dG.NewVec(3);
+  dgamdr.NewVec(3*n_atoms);
+
+
+  for (int i = 0, ii = 0; i < mol.NAtoms; i++) {
+    ii = realIdx(i);
+    if (ii < 0) continue;
+    izp = mol.ATNO(i);
+    // Effective charge width of i
+    norm_cn = 1.0 / pow(avg_cn[izp], norm_exp);
+    radi = rad[izp] * (1.0 - kcnrad * cn(ii) * norm_cn);
+    dradcn_i = - rad[izp] * kcnrad * norm_cn;
+    for (int j = 0, jj = 0; j < i; j++) {
+      jj = realIdx(j);
+      if (jj < 0) continue;
+      // Effective charge width of j
+      jzp = mol.ATNO(j);
+      norm_cn = 1.0 / pow(avg_cn[jzp], norm_exp);
+      radj = rad[jzp] * (1.0 - kcnrad * cn(jj) * norm_cn);
+      dradcn_j = - rad[jzp] * kcnrad * norm_cn;
+      gam = 1.0 / (pow(radi,2) + pow(radj,2));
+      for (int k = 0, kk = 0; k < mol.NAtoms; k++) {
+        kk = realIdx(k);
+        if (kk < 0) continue;
+        // Coulomb interaction of Gaussian charges
+        for (int c = 0; c < 3; c++) {
+          vec(c) = mol.CC(jj, c) - mol.CC(ii, c);
+          dgamdr(3*kk+c) =  -(radi * dradcn_i * dcndr(ii, 3*kk + c) + radj * dradcn_j * dcndr(jj, 3*kk + c)) * pow(gam, 3);
+        }
+      }
+      // Explicit derivative
+      r2 = pow(vec(0), 2) + pow(vec(1), 2) + pow(vec(2), 2);
+      arg = gam * gam * r2;
+      dtmp1 = 2.0 * gam * exp(-arg) / (sqrtpi * r2) - erf(sqrt(arg)) / (r2 *sqrt(r2)); 
+      // Effective charge width derivative
+      dtmp2 = 2.0 * exp(-arg) / sqrtpi;
+      // Capacitance derivative off-diagonal
+      dtmp3 = erf(sqrt(r2) * gam) / sqrt(r2);
+      // Capacitance derivative diagonal
+      dtmp4 = (eta[izp] + kqeta[izp] * qloc(ii) + sqrt2pi / radi) * q(ii);
+      dtmp5 = (eta[jzp] + kqeta[jzp] * qloc(jj) + sqrt2pi / radj) * q(jj);
+      for (int c = 0; c < 3; c++) {
+        // Explicit derivative
+        dG(c) = dtmp1 * vec(c);
+        atrace(ii, c) += -dG(c) * q(jj) * cmat(ii, jj);
+        atrace(jj, c) += +dG(c) * q(ii) * cmat(jj, ii);
+        dAmatdr(3*ii+c, jj) += -dG(c) * q(ii) * cmat(jj, ii);
+        dAmatdr(3*jj+c, ii) += +dG(c) * q(jj) * cmat(ii, jj);
+        // Effective charge width derivative
+        atrace(ii, c) += -dtmp2 * q(jj) * dgamdr(3*jj+c) * cmat(ii, jj);
+        atrace(jj, c) += -dtmp2 * q(ii) * dgamdr(3*ii+c) * cmat(jj, ii);
+        dAmatdr(3*ii+c, jj) += +dtmp2 * q(ii) * dgamdr(3*ii+c) * cmat(jj, ii);
+        dAmatdr(3*jj+c, ii) += +dtmp2 * q(jj) * dgamdr(3*jj+c) * cmat(ii, jj);
+        // Capacitance derivative off-diagonal
+        atrace(ii, c) += -dtmp3 * q(jj) * dcmatdr(ii, 3*jj+c);
+        atrace(jj, c) += -dtmp3 * q(ii) * dcmatdr(jj, 3*ii+c);
+        dAmatdr(3*ii+c, jj) += +dtmp3 * q(ii) * dcmatdr(jj, 3*ii+c);
+        dAmatdr(3*jj+c, ii) += +dtmp3 * q(jj) * dcmatdr(ii, 3*jj+c);
+        // Capacitance derivative diagonal
+        dAmatdr(3*jj+c, ii) += -dtmp4 * dcmatdr(ii, 3*jj+c);
+        dAmatdr(3*ii+c, jj) += -dtmp5 * dcmatdr(jj, 3*ii+c);
+      }
+    }  // jj
+      dtmp1 = kqeta[izp] * q(ii) * cmat(ii, ii);  // Hardness
+      dtmp2 = -sqrt2pi * dradcn_i / pow(radi, 2) * q(ii) * cmat(ii, ii);  // Effective charge width
+      for (int k = 0, kk = 0; k < mol.NAtoms; k++) {
+        kk = realIdx(k);
+        if (kk < 0) continue;
+        for (int c = 0; c < 3; c++) {
+          // Hardness derivative
+          dAmatdr(3*kk+c, ii) += +dtmp1 * dqlocdr(ii, 3*kk+c);
+          // Effective charge width derivative
+          dAmatdr(3*kk+c, ii) += +dtmp2 * dcndr(ii, 3*kk+c);
+        }
+      }
+      // Capacitance derivative
+      dtmp3 = (eta[izp] + kqeta[izp] * qloc(ii) + sqrt2pi / radi) * q(ii);
+      for (int c = 0; c < 3; c++) {
+        dAmatdr(3*ii+c, ii) += +dtmp3 * dcmatdr(ii, 3*ii+c);
+      }
+}  // ii
+
   
   return EXIT_SUCCESS;
 };
@@ -517,13 +686,13 @@ int EEQBCModel::get_qloc(
   const TIVector &realIdx,
   const TMatrix<double> &dist,
   const double q_tot,  // total system charge
-  TVector<double> &qloc
+  TVector<double> &qloc,  // Local charge
+  TMatrix<double> &dqlocdr,  // Derivative of local charge w.r.t. atom positions
+  bool lgrad
 ) 
 {
   const double cutoff = 25.0;
-  bool lgrad = false;
   TVector<double> cn;
-  TMatrix<double> dcndr;
   // Electronegativity scaled coordination number with EEQ-BC parameters
   NCoordErfEN ncoord_erf_en(
     ncoorderf_kcn,
@@ -531,10 +700,11 @@ int EEQBCModel::get_qloc(
     ncoorderf_cutoff,
     ncoorderfen_f_directed,
     ncoorderf_cn_max);
-  qloc.NewVector(mol.NAtoms);
-  const double q_tot_norm = q_tot/mol.NAtoms;
+  const int n_atoms = realIdx.Max() + 1;  // Number of atoms
+  qloc.NewVector(n_atoms);
+  const double q_tot_norm = q_tot/n_atoms;
 
-  ncoord_erf_en.get_ncoord(mol, dist, cn, dcndr, lgrad);
+  ncoord_erf_en.get_ncoord(mol, dist, cn, dqlocdr, lgrad);
 
   for (int i = 0, ii = 0; i < mol.NAtoms; i++) {
     ii = realIdx(i);
@@ -588,6 +758,26 @@ int EEQBCModel::get_cpair(
   return EXIT_SUCCESS;
 }
 
+// Get derivative of the capacitance for bond between atoms i and j for EEQ-BC 
+int EEQBCModel::get_dcpair(
+  double dist_ij,  // distance between atom j to atom i
+  double rvdw_ijat,  // pairwise van der Waals radii for ij atom types
+  double cap_ij,  // product of bond capacitances for atom types of i and j
+  TVector<double> &vec, // Vector from j to i
+  TVector<double> &dcdr_ij  // Out: Capacitance for bond ij
+) const {
+  // Calculate the argument of the error function
+  //     arg = -(kbc * (r1      - rvdw       ) / rvdw)**2
+  double arg =   kbc * (dist_ij - rvdw_ijat) / rvdw_ijat;
+  //     dtmp = sqrt(capi * capj) * kbc * exp(arg) / (sqrtpi * rvdw)
+  double dtmp = sqrt(cap_ij) * kbc * exp( - pow(arg, 2)) / (sqrt(pi) * rvdw_ijat);
+  for (int c = 0; c < 3; c++) {
+    dcdr_ij(c) = dtmp * vec(c) / dist_ij;
+  }
+
+  return EXIT_SUCCESS;
+}
+
 // Get the capacitance matrix
 int EEQBCModel::get_cmat(
   const TMolecule &mol,  // molecular geometry
@@ -600,6 +790,7 @@ int EEQBCModel::get_cmat(
   double c_ij;  // Capacitance for bond between atoms i and j
   double dist_ij; // distance between atoms i and j
   const int n_atoms = realIdx.Max() + 1;
+  cmat.NewMatrix(n_atoms + 1, n_atoms + 1);
   for (int i = 0, ii = 0; i < mol.NAtoms; i++)
   {
     ii = realIdx(i);
@@ -624,6 +815,61 @@ int EEQBCModel::get_cmat(
   return EXIT_SUCCESS;
 }
 
+// Get the derivative of the capacitance matrix
+int EEQBCModel::get_dcmatdr(
+  const TMolecule &mol,  // molecular geometry
+  const TIVector &realIdx,  // The real atom indices (for excluding dummy atoms)
+  const TMatrix<double> &dist, // atom distances
+  TMatrix<double> &dcmatdr  // Out: Capacitance matrix
+) {
+  int iat;  // atom type of i
+  int jat;  // atom type of j
+  int ij_min, ij_max; // min and max from i and j
+  int ij_at;  // pairwise index for atom types of i and j
+  double c_ij;  // Capacitance for bond between atoms i and j
+  double dist_ij; // distance between atoms i and j
+  double rvdw_ijat;  // pairwise van der Waals radii for ij atom types
+  double cap_ij;  // product of bond capacitances for atom types of i and j
+  TVector<double> vec; // Vector from i to j
+  vec.NewVec(3);
+  TVector<double> dcdr_ij; // Part ij of capacitance derivative
+  dcdr_ij.NewVec(3);
+  const int n_atoms = realIdx.Max() + 1;
+  dcmatdr.NewMat(n_atoms, 3 * n_atoms);
+  for (int i = 0, ii = 0; i < mol.NAtoms; i++)
+  {
+    ii = realIdx(i);
+    if (ii < 0) continue;
+    iat = mol.ATNO(i);
+    for (int j = 0, jj = 0; j < i; j++)
+    {
+      jj = realIdx(j);
+      if (jj < 0) continue;
+      jat = mol.ATNO(j);
+      for (int c = 0; c < 3; c++) {
+        vec(c) = mol.CC(jj, c) - mol.CC(ii, c);
+      }
+      dist_ij = dist(ii,jj);
+      ij_min = std::min(iat, jat);
+      ij_max = std::max(iat, jat);
+      ij_at = ij_min + ij_max * (ij_max - 1)/2 - 1;
+      rvdw_ijat = rvdw[ij_at];
+      cap_ij = cap[iat] * cap[jat];
+      get_dcpair( dist_ij, rvdw_ijat, cap_ij, vec, dcdr_ij);
+      for (int c = 0; c < 3; c++) {
+      // Calculate Off-diagonal elements; bond capacitances
+      dcmatdr(jj, 3*ii+c) = - dcdr_ij(c);
+      dcmatdr(ii, 3*jj+c) = + dcdr_ij(c);
+      // Calculate diagonal elements; self-capacitance as the negative sum of bond capacitances
+      dcmatdr(ii, 3*ii+c) = dcmatdr(ii, 3*ii+c) + dcdr_ij(c);
+      dcmatdr(jj, 3*jj+c) = dcmatdr(jj, 3*jj+c) - dcdr_ij(c);
+      }
+    }
+  }
+
+  return EXIT_SUCCESS;
+}
+
 // Get the right-hand side vector b
 // A * q = b with Coulomb matrix A and partial charges q
 // b = X * C with electronegativity vector X and capacitance matrix C
@@ -634,16 +880,19 @@ int EEQBCModel::get_xvec(
   const TVector<double> &cn,
   TMatrix<double> &cmat, // capacitance matrix
   int charge,  // total charge of the system
+  TVector<double> &qloc,
   TVector<double> &xvec  // Out: electronegativity vector
 ) {
   int info{0};
   const int n_atoms = realIdx.Max() + 1;
   TVector<double> x_tmp;  // dummy for xvec, has dimension N+1 including the constraint
+  TMatrix<double> dqlocdr;
   x_tmp.NewVector(n_atoms+1);
   int i_atno;  // atomic number of atom i
 
   // get local charge
-  info = get_qloc(mol, realIdx, dist, charge, qloc);
+  info = get_qloc(mol, realIdx, dist, charge, qloc, dqlocdr, false);
+  if (info != EXIT_SUCCESS) return info;
 
   for (int i = 0, ii = 0; i < mol.NAtoms; i++)
   {
@@ -655,65 +904,112 @@ int EEQBCModel::get_xvec(
   
   xvec.NewVector(n_atoms + 1);
   xvec(n_atoms) = charge;
-  for (int i = 0; i < n_atoms; i++)
+  for (int i = 0, ii = 0; i < mol.NAtoms; i++)
   {
-    for (int j = 0; j < n_atoms; j++)
+    ii = realIdx(i);
+    if (ii < 0) continue;
+    for (int j = 0, jj = 0; j < mol.NAtoms; j++)
     {
-      xvec(i) = xvec(i) + cmat(i, j) * x_tmp(j);
+      jj = realIdx(j);
+      if (jj < 0) continue;
+      xvec(ii) = xvec(ii) + cmat(ii, jj) * x_tmp(jj);
     }
   }
 
   return EXIT_SUCCESS; 
 }
 
-// numerical gradient of partial charges w.r.t. atom positions
-int EEQBCModel::num_grad_dqdr(
-  TMolecule &mol,  // molecular geometry
+// Get the right-hand side vector b and its derivative w.r.t. CN
+// A * q = b with Coulomb matrix A and partial charges q
+// b = X * C with electronegativity vector X and capacitance matrix C
+int EEQBCModel::get_xvec_derivs(
+  const TMolecule &mol,  // molecular geometry
   const TIVector &realIdx,  // The real atom indices (for excluding dummy atoms)
-  int charge, // total charge of the system
-  TMatrix<double> &num_dqdr // numerical gradient
+  const TMatrix<double> &dist,  // atom distances
+  const TVector<double> &cn,  // coordination number (CN)
+  const TMatrix<double> &dcndr, // coordination number derivative
+  TMatrix<double> &cmat, // capacitance matrix
+  int charge,  // total charge of the system
+  TVector<double> &xvec,  // Out: electronegativity vector
+  TMatrix<double> &dxvecdr,  // derivative of the entire right hand side w.r.t. atom positions
+  TVector<double> &qloc,
+  TMatrix<double> &dqlocdr,  // derivative of qloc w.r.t. atom positions
+  TMatrix<double> &dcmatdr
 ) {
-  TVector<double> q_r, q_l;  // forward and backward point
-  double step{1.0e-6};  // step size of finite differences
-  const int nat = realIdx.Max() +1;  // number of atoms
-  TVector<double> cn;  // coordination number
-  TMatrix<double> dcndr;  // derivative of the coordination number
-  TMatrix<double> dqdr;  // dummy variable for analytical derivative
-  dqdr.NewMat(3 * nat, nat);
-  num_dqdr.NewMat(3 * nat, nat);
-  TMatrix<double> dist;  // matrix with pairwise atomic distances
-  multicharge::EEQBCModel eeqbc_model;  // instance of the EEQ-BC model
+  int info{0};
+  const int n_atoms = realIdx.Max() + 1;
+  int i_atno;  // atomic number of atom i
+  TVector<double> x_tmp;  // dummy for xvec, has dimension N+1 including the constraint
+  TMatrix<double> dxvecdr_tmp;
+  TMatrix<double> cmat_tmp;
+  x_tmp.NewVector(n_atoms+1);
+  dxvecdr_tmp.NewMat(n_atoms, 3*n_atoms);
+  cmat_tmp.NewMat(n_atoms, n_atoms);
+  dxvecdr.NewMat(n_atoms, 3*n_atoms);
 
-  // calculate numerical gradient via finite difference method
-  for (int i = 0, ii = 0; i < mol.NAtoms; i++) {
+  // calculate derivative of the capacitance
+  info = get_dcmatdr(mol, realIdx, dist, dcmatdr);
+  if (info != EXIT_SUCCESS) return info;
+
+  // get local charge
+  bool lgrad = true; // calculate dqlocdr
+  info = get_qloc(mol, realIdx, dist, charge, qloc, dqlocdr, lgrad);
+  if (info != EXIT_SUCCESS) return info;
+  for (int i=0, ii = 0; i < mol.NAtoms; i++)
+  {
+    ii = realIdx(i);
+    if (ii < 0) continue;
+    i_atno = mol.ATNO(i);
+    for (int j = 0, jj = 0; j < mol.NAtoms; j++)
+    {
+      jj = realIdx(j);
+      if (jj < 0) continue;
+      // setup C-matrix with correct size
+      cmat_tmp(ii,jj) = cmat(ii,jj);
+      for (int c = 0; c < 3; c++)
+      {
+        dxvecdr_tmp(ii, 3*jj+c) = kcnchi[i_atno] * dcndr(ii, 3*jj+c) + dxvecdr_tmp(ii, 3*jj+c);
+        dxvecdr_tmp(ii, 3*jj+c) = kqchi[i_atno] * dqlocdr(ii, 3*jj+c) + dxvecdr_tmp(ii, 3*jj+c);
+      }
+    }
+  }
+  
+  BLAS_Add_Mat_x_Mat(dxvecdr, cmat_tmp, dxvecdr_tmp, false, false, 1.0);
+
+  for (int i = 0, ii = 0; i < mol.NAtoms; i++)
+  {
     ii = realIdx(i);
     if (ii < 0) continue;
-    for (int c = 0; c < 3; c++) {
-      // calculate forward point
-      mol.CC(i, c) += step;
-      dist.NewMatrix(nat, nat);
-      calc_distances(mol, realIdx, dist);
-      eeqbc_model.get_cn(mol, realIdx, dist, cn, dcndr, false);
-      q_r.NewVec(nat);
-      eeqbc_model.eeq_chrgeq(mol, realIdx, dist, cn, dcndr, charge, q_r, dqdr, false, false);
-
-      // calculate backward point
-      mol.CC(i, c) = mol.CC(i, c) - 2 * step;
-      dist.NewMatrix(nat, nat);
-      calc_distances(mol, realIdx, dist);
-      eeqbc_model.get_cn(mol, realIdx, dist, cn, dcndr, false);
-      q_l.NewVec(nat);
-      eeqbc_model.eeq_chrgeq(mol, realIdx, dist, cn, dcndr, charge, q_l, dqdr, false, false);
-
-      // calculate numerical gradient as finite difference
-      mol.CC(i, c) = mol.CC(i, c) + step;
-      for (int j = 0, jj = 0; j < nat; j++) {
-        num_dqdr(3 * ii + c, j) = 0.5 * (q_r(j) - q_l(j)) / step;
+    i_atno = mol.ATNO(i);
+    x_tmp(ii) = - chi[i_atno] + kcnchi[i_atno]*cn(ii) + kqchi[i_atno]*qloc(ii);
+  }
+  
+  xvec.NewVector(n_atoms + 1);
+  xvec(n_atoms) = charge;
+  BLAS_Add_Mat_x_Vec(xvec, cmat, x_tmp, false, 1.0);
 
+  for (int i = 0, ii = 0; i < mol.NAtoms; i++)
+  {
+    ii = realIdx(i);
+    if (ii < 0) continue;
+    for (int j = 0, jj = 0; j < i; j++)
+    {
+      jj = realIdx(j);
+      if (jj < 0) continue;
+      for (int c = 0; c < 3; c++) {
+        // setup dCij/dR * Xj
+        dxvecdr(ii, 3*ii+c) += x_tmp(jj) * dcmatdr(jj, 3*ii+c);
+        dxvecdr(jj, 3*jj+c) += x_tmp(ii) * dcmatdr(ii, 3*jj+c);
+        dxvecdr(jj, 3*ii+c) += (x_tmp(ii)-x_tmp(jj)) * dcmatdr(jj, 3*ii+c);
+        dxvecdr(ii, 3*jj+c) += (x_tmp(jj)-x_tmp(ii)) * dcmatdr(ii, 3*jj+c);
       }
     }
+    dxvecdr(ii, 3*ii  ) += x_tmp(ii) * dcmatdr(ii, 3*ii  );
+    dxvecdr(ii, 3*ii+1) += x_tmp(ii) * dcmatdr(ii, 3*ii+1);
+    dxvecdr(ii, 3*ii+2) += x_tmp(ii) * dcmatdr(ii, 3*ii+2);
   }
-  return EXIT_SUCCESS;
+
+  return EXIT_SUCCESS; 
 }
 
 } // namespace multicharge
diff --git a/src/dftd_ncoord.cpp b/src/dftd_ncoord.cpp
index 6b0eb8a..ec229d5 100644
--- a/src/dftd_ncoord.cpp
+++ b/src/dftd_ncoord.cpp
@@ -267,18 +267,18 @@ int NCoordBase::dr_ncoord_base(
       cn(ii) += countf;
       cn(jj) += countf * f_directed;
 
-      dcountf = f_en * dr_count_fct(r, rcovij) / rcovij;
+      dcountf = f_en * dr_count_fct(r, rcovij) / pow(rcovij, norm_exp);
       dcndr(jj, 3 * jj    ) -= dcountf * rx * f_directed;
       dcndr(jj, 3 * jj + 1) -= dcountf * ry * f_directed;
       dcndr(jj, 3 * jj + 2) -= dcountf * rz * f_directed;
 
-      dcndr(jj, 3 * ii    ) += dcountf * rx;
-      dcndr(jj, 3 * ii + 1) += dcountf * ry;
-      dcndr(jj, 3 * ii + 2) += dcountf * rz;
+      dcndr(jj, 3 * ii    ) += dcountf * rx * f_directed;
+      dcndr(jj, 3 * ii + 1) += dcountf * ry * f_directed;
+      dcndr(jj, 3 * ii + 2) += dcountf * rz * f_directed;
 
-      dcndr(ii, 3 * jj    ) -= dcountf * rx * f_directed;
-      dcndr(ii, 3 * jj + 1) -= dcountf * ry * f_directed;
-      dcndr(ii, 3 * jj + 2) -= dcountf * rz * f_directed;
+      dcndr(ii, 3 * jj    ) -= dcountf * rx;
+      dcndr(ii, 3 * jj + 1) -= dcountf * ry;
+      dcndr(ii, 3 * jj + 2) -= dcountf * rz;
 
       dcndr(ii, 3 * ii    ) += dcountf * rx;
       dcndr(ii, 3 * ii + 1) += dcountf * ry;
diff --git a/test/molecules.h b/test/molecules.h
index ba819b1..4b27cf6 100644
--- a/test/molecules.h
+++ b/test/molecules.h
@@ -363,4 +363,125 @@ static const double actinides_coord[actinides_n * 3]{
   +2.01082807362334, +1.34838807211157, -4.70482633508447,
 };
 
+// Amalgam
+static const int amalgam_n{56};
+static const int amalgam_charge{0};
+static const char amalgam_atoms[amalgam_n][3]{
+ "he",
+ "ne",
+ "ar",
+ "kr",
+ "xe",
+ "rn",
+ "li",
+ "b",
+ "be",
+ "c",
+ "h",
+ "o",
+ "h",
+ "n",
+ "h",
+ "o",
+ "c",
+ "c",
+ "f",
+ "cl",
+ "c",
+ "c",
+ "c",
+ "h",
+ "c",
+ "br",
+ "at",
+ "i",
+ "na",
+ "al",
+ "mg",
+ "h",
+ "si",
+ "p",
+ "h",
+ "s",
+ "h",
+ "c",
+ "c",
+ "h",
+ "h",
+ "o",
+ "h",
+ "h",
+ "o",
+ "h",
+ "h",
+ "se",
+ "h",
+ "h",
+ "te",
+ "h",
+ "h",
+ "po",
+ "h",
+ "h"
+};
+
+static const double amalgam_coord[amalgam_n * 3]{
+   -8.24412055,   -6.58780979,   -0.91570650,
+   -6.42286818,    3.90712547,    4.86500198,
+   -4.90737238,   -1.95282939,    7.86825624,
+   -0.37726965,    4.92890868,    0.95588204,
+   13.38555226,    2.11718711,    1.81147578,
+   13.15321224,    2.33063559,   -5.87115854,
+   -7.65859324,   -3.86785454,   -3.03173103,
+   -8.90082657,   -1.51503593,    1.95474829,
+   -7.15048995,   -0.24936176,   -0.62038223,
+   -4.29597130,   -1.27062930,    0.62721393,
+  -10.72806234,   -1.93304209,    3.29657114,
+   -2.08637070,   -0.49174028,    0.81639801,
+   -7.07902488,   -3.71571695,    4.64117407,
+   -4.59992498,   -3.76766281,    1.88857834,
+   -3.52261371,   -3.71637066,    3.48719798,
+   -7.19030864,   -3.91148428,    2.82109087,
+    0.40738064,   -2.42265922,   -3.49890723,
+    1.60194914,   -0.09413544,   -3.65055172,
+    2.22710681,    3.76932919,   -5.71850741,
+   -1.36786746,    3.07862751,   -9.79436860,
+    0.99258997,    1.54015982,   -5.57686673,
+   -0.78697416,    0.96614072,   -7.40879667,
+   -2.02138722,   -1.34900108,   -7.21789800,
+    2.97566733,    0.47460381,   -2.23605219,
+   -1.43901063,   -3.02645297,   -5.27319702,
+   -4.52048514,   -2.13126179,   -9.69487380,
+    1.50705279,   -5.00133236,   -0.44613628,
+   -3.34847853,   -6.54426878,   -4.98924265,
+    7.39326107,   -4.78317478,   -2.73955769,
+    6.73805630,    1.57365902,   -5.61518652,
+    6.57487539,    4.43501257,   -1.01996941,
+    7.67481939,   -1.07984981,   -2.56613001,
+    7.07714209,   -0.94539154,   -9.82369772,
+    8.06046833,   -4.55945339,   -8.32690067,
+    4.24002482,   -1.09444619,   -9.48623334,
+    4.83030003,   -6.50824697,   -6.76694609,
+    3.16857808,   -6.16789369,   -8.60564090,
+    5.65038102,  -12.11407383,   -2.16358858,
+    7.49828789,  -10.84757354,   -1.04222812,
+    3.98947118,  -11.17048194,   -2.88947999,
+    5.72477568,  -14.14203461,   -2.34561989,
+    7.47130878,   -8.31368664,   -0.71866358,
+    9.36082193,   -7.70479742,   -0.60920006,
+    9.17502500,  -11.80483919,   -0.30762578,
+   -6.62601245,    3.21948915,   -0.72736379,
+   -6.32322321,    3.63571554,   -2.46098961,
+   -5.00474256,    3.33318462,    0.10681509,
+  -11.53711695,   -5.76950992,   -5.74618428,
+  -11.11536735,   -4.21494504,   -7.93442044,
+  -10.46118337,   -7.94640025,   -6.97071464,
+   12.82621494,   -6.62177000,   -1.60795607,
+   11.38819319,   -4.20140874,   -2.84890169,
+   13.91845142,   -4.62346343,    0.64687720,
+   10.24369296,   -9.23027806,   -6.33500557,
+   12.98362119,   -8.51243817,   -7.95590708,
+   11.88048032,  -12.04935553,   -5.38433128
+};
+
 #endif // MOLECULES_H
diff --git a/test/test_grad.cpp b/test/test_grad.cpp
index d2ecdd5..11d9226 100644
--- a/test/test_grad.cpp
+++ b/test/test_grad.cpp
@@ -169,8 +169,7 @@ int test_tpss0d4mbd_rost61m1() {
   return test_numgrad(mol, charge, par);
 }
 
-
-int test_numgrad_dqdr(
+int test_numgrad_dqdr_eeq(
   int n,
   const char atoms[][3],
   const double coord[]
@@ -243,7 +242,7 @@ int test_numgrad_dqdr(
     for (int c = 0; c < 3; c++) {
       for (int j = 0; j < mol.NAtoms; j++) {
         if (check(analytic_dqdr(3 * i + c, j), num_dqdr(3 * i + c, j), thr) != EXIT_SUCCESS) {
-          print_fail("Gradient mismatch for dqdr\n", analytic_dqdr(3 * i + c, j), num_dqdr(3 * i + c, j));
+          print_fail("Gradient mismatch for dqdr in EEQ.\n", analytic_dqdr(3 * i + c, j), num_dqdr(3 * i + c, j));
           return EXIT_FAILURE;
         }
       }
@@ -253,11 +252,95 @@ int test_numgrad_dqdr(
   return EXIT_SUCCESS;
 }
 
+int test_numgrad_dqdr_eeqbc(
+  int n,
+  const char atoms[][3],
+  const double coord[]
+) {
+  // assemble molecule
+  int info;
+  TMolecule mol;
+  info = get_molecule(n, atoms, coord, mol);
+  if (info != EXIT_SUCCESS) return info;
+  TVector<double> q_r, q_l;
+  double step{1.0e-6}; // accurate up to 1.0E-8
+  double thr{4.0e-8};
+
+  TMatrix<double> dist;
+  dist.NewMat(mol.NAtoms, mol.NAtoms);
+  TMatrix<double> num_dqdr;  // numerical gradient of the partial charges
+  num_dqdr.NewMat(3 * mol.NAtoms, mol.NAtoms);
+  TMatrix<double> analytic_dqdr;  // analytical gradient of the partial charges
+  analytic_dqdr.NewMat(3 * mol.NAtoms, mol.NAtoms);
+  multicharge::EEQBCModel eeqbc_model;
+  TVector<double> q;
+  TMatrix<double> dqdr;
+  q.NewVec(mol.NAtoms);
+  dqdr.NewMat(3 * mol.NAtoms, mol.NAtoms);
+
+  TCutoff cutoff;
+
+  // masking (nothing excluded)
+  TVector<int> realIdx;
+  realIdx.NewVec(mol.NAtoms);
+  int nat = 0;
+  for (int i = 0; i != mol.NAtoms; i++) {
+    realIdx(i) = nat;
+    nat++;
+  }
+
+  // analytical gradient
+  calc_distances(mol, realIdx, dist);
+  info =
+    eeqbc_model.get_charges(mol, realIdx, dist, 0, cutoff.cn_eeq, q, analytic_dqdr, true);
+  if (info != EXIT_SUCCESS) return info;
+
+  // calculate numerical gradient via finite difference method
+  for (int i = 0; i < mol.NAtoms; i++) {
+    for (int c = 0; c < 3; c++) {
+      // calculate forward point
+      q_r.NewVec(n);
+      mol.CC(i, c) += step;
+      calc_distances(mol, realIdx, dist);
+      eeqbc_model.get_charges(mol, realIdx, dist, 0, cutoff.cn_eeq, q_r, dqdr, false);
+
+      // calculate backward point
+      q_l.NewVec(n);
+      mol.CC(i, c) = mol.CC(i, c) - 2 * step;
+      calc_distances(mol, realIdx, dist);
+      eeqbc_model.get_charges(mol, realIdx, dist, 0, cutoff.cn_eeq, q_l, dqdr, false);
+
+      // calculate numerical gradient as finite difference
+      mol.CC(i, c) = mol.CC(i, c) + step;
+      for (int j = 0; j < mol.NAtoms; j++) {
+        num_dqdr(3 * i + c, j) = 0.5 * (q_r(j) - q_l(j)) / step;
+      }
+    }
+  }
+
+  // compare against numerical gradient
+  for (int i = 0; i < mol.NAtoms; i++) {
+    for (int c = 0; c < 3; c++) {
+      for (int j = 0; j < mol.NAtoms; j++) {
+        if (check(analytic_dqdr(3 * i + c, j), num_dqdr(3 * i + c, j), thr) != EXIT_SUCCESS) {
+          print_fail("Gradient mismatch for dqdr in EEQ-BC.\n", analytic_dqdr(3 * i + c, j), num_dqdr(3 * i + c, j));
+          return EXIT_FAILURE;
+        }
+      }
+    }
+  }
+
+  return EXIT_SUCCESS;
+}
 
 int test_grad() {
   int info{0};
 
-  info = test_numgrad_dqdr(
+  info = test_numgrad_dqdr_eeq(
+    mb16_43_01_n, mb16_43_01_atoms, mb16_43_01_coord);
+  if (info != EXIT_SUCCESS) return info;
+
+  info = test_numgrad_dqdr_eeqbc(
     mb16_43_01_n, mb16_43_01_atoms, mb16_43_01_coord);
   if (info != EXIT_SUCCESS) return info;
 
diff --git a/test/test_multi.cpp b/test/test_multi.cpp
index 502d281..5f17064 100644
--- a/test/test_multi.cpp
+++ b/test/test_multi.cpp
@@ -54,6 +54,7 @@ int test_multi_functions(){
   info = get_molecule(mb16_43_01_n, mb16_43_01_atoms, mb16_43_01_coord, mol);
   if (info != EXIT_SUCCESS) {
      printf("Multicharge: Functions, Failed to set up molecule.");
+     fflush(stdout);
      return info;
   }
   
@@ -71,12 +72,13 @@ int test_multi_functions(){
   TVector<double> q;
   TMatrix<double> dqdr;
   q.NewVector(mb16_43_01_n);
-  dqdr.NewMatrix(mb16_43_01_n, mb16_43_01_n);
+  dqdr.NewMatrix(3*mb16_43_01_n, mb16_43_01_n);
 
-  eeqbc_model.get_cn(mol, realIdx, dist, cn, dcndr, false);
-  info = eeqbc_model.eeq_chrgeq(mol, realIdx, dist, cn, dcndr, mb16_43_01_charge, q, dqdr, false, false);
+  eeqbc_model.get_cn(mol, realIdx, dist, cn, dcndr, true);
+  info = eeqbc_model.eeq_chrgeq(mol, realIdx, dist, cn, dcndr, mb16_43_01_charge, q, dqdr, true, false);
   if (info != EXIT_SUCCESS) {
      printf("Multicharge: Functions, Failed to calculate charges.");
+     fflush(stdout);
      return info;
   }
 
@@ -91,6 +93,54 @@ int test_multi_functions(){
   return EXIT_SUCCESS;
 }
 
+// Test member functions of ChargeModel and derived classes
+int test_eeqbc_amalgam(){
+  int info;
+  // assemble molecule
+  TMolecule mol;
+  info = get_molecule(amalgam_n, amalgam_atoms, amalgam_coord, mol);
+  if (info != EXIT_SUCCESS) {
+     printf("Multicharge: Functions, Failed to set up molecule.");
+     fflush(stdout);
+     return info;
+  }
+  
+  // Test EEQ-BC member functions
+  multicharge::EEQBCModel eeqbc_model;
+  TIVector realIdx;
+  TMatrix<double> dist;
+  dftd4::initializeRealIdx(amalgam_n, realIdx);
+  dist.NewMatrix(amalgam_n, amalgam_n);
+  dftd4::calc_distances(mol, realIdx, dist);
+  TVector<double> cn;
+  TMatrix<double> dcndr;
+
+  // Calculate partial charges
+  TVector<double> q;
+  TMatrix<double> dqdr;
+  q.NewVector(amalgam_n);
+  dqdr.NewMatrix(3*amalgam_n, amalgam_n);
+
+  eeqbc_model.get_cn(mol, realIdx, dist, cn, dcndr, true);
+
+  info = eeqbc_model.eeq_chrgeq(mol, realIdx, dist, cn, dcndr, amalgam_charge, q, dqdr, false, false);
+  if (info != EXIT_SUCCESS) {
+     printf("Multicharge: Failed to calculate charges.");
+     fflush(stdout);
+     return info;
+  }
+
+  // Check against multicharge reference calculation
+  for (int i=0; i < mol.NAtoms; i++) {
+    if (check(q(i), qvec_amalgam_reference[i], 1.0E-8) == EXIT_FAILURE) {
+        print_fail("Multicharge: Functions, Partial charge differs from reference for amalgam.", q(i), qvec_amalgam_reference[i]);
+        return EXIT_FAILURE;
+    }
+  }
+
+  return EXIT_SUCCESS;
+}
+
 // Test for multicharge models
 int test_multi(){
   int info;
@@ -103,5 +153,9 @@ int test_multi(){
   info = test_multi_functions();
   if (info != EXIT_SUCCESS) return info;
 
+  // Test amalgam structure
+  info = test_eeqbc_amalgam();
+  if (info != EXIT_SUCCESS) return info;
+
   return EXIT_SUCCESS;
 }
diff --git a/test/test_multi.h b/test/test_multi.h
index 2935c13..caaf1aa 100644
--- a/test/test_multi.h
+++ b/test/test_multi.h
@@ -20,6 +20,340 @@ static const double qvec_reference[]{
     -0.0144595425,  0.2577820828, -0.1117775795,  0.4834861246
   };
 
+// reference gradient dqdr for mb16_43_01
+static const double dqdr_reference[]{
+     0.000000000,     0.000000000,     0.000000000, // (0-2,0,0)
+     0.119414865,     0.119414865,     0.119414865,
+    -0.009766387,    -0.009766387,    -0.009766387,
+    -0.129430372,    -0.129430372,    -0.129430372,
+    -0.026797289,    -0.026797289,    -0.026797289,
+     0.034144659,     0.034144659,     0.034144659,
+    -0.007388839,    -0.007388839,    -0.007388839,
+     0.032304539,     0.032304539,     0.032304539,
+     0.010374760,     0.010374760,     0.010374760,
+    -0.037425485,    -0.037425485,    -0.037425485,
+    -0.003313457,    -0.003313457,    -0.003313457,
+     0.017425270,     0.017425270,     0.017425270,
+    -0.006572642,    -0.006572642,    -0.006572642,
+     0.001862799,     0.001862799,     0.001862799,
+    -0.004669860,    -0.004669860,    -0.004669860,
+     0.006015602,     0.006015602,     0.006015602,
+
+     0.003821836,     0.003821836,     0.003821836, // (0-2,0,1)
+     0.004036443,     0.004036443,     0.004036443, // (0-2,1,1)
+     0.122929791,     0.122929791,     0.122929791,
+     0.001205792,     0.001205792,     0.001205792,
+    -0.026345244,    -0.026345244,    -0.026345244,
+     0.001145827,     0.001145827,     0.001145827,
+     0.000982307,     0.000982307,     0.000982307,
+     0.035008816,     0.035008816,     0.035008816,
+     0.008471227,     0.008471227,     0.008471227,
+     0.000693925,     0.000693925,     0.000693925,
+    -0.002006854,    -0.002006854,    -0.002006854,
+    -0.007359077,    -0.007359077,    -0.007359077,
+    -0.010110234,    -0.010110234,    -0.010110234,
+    -0.149754592,    -0.149754592,    -0.149754592,
+    -0.001580143,    -0.001580143,    -0.001580143,
+     0.005455892,     0.005455892,     0.005455892,
+
+     0.017226125,     0.017226125,     0.017226125, // (0-2,0,2)
+    -0.032435186,    -0.032435186,    -0.032435186,
+    -0.002954316,    -0.002954316,    -0.002954316,
+     0.091498968,     0.091498968,     0.091498968,
+    -0.003224121,    -0.003224121,    -0.003224121,
+     0.091854054,     0.091854054,     0.091854054,
+    -0.003845758,    -0.003845758,    -0.003845758,
+     0.003182059,     0.003182059,     0.003182059,
+     0.006322413,     0.006322413,     0.006322413,
+    -0.022201436,    -0.022201436,    -0.022201436,
+    -0.002890625,    -0.002890625,    -0.002890625,
+     0.009411540,     0.009411540,     0.009411540,
+    -0.000828985,    -0.000828985,    -0.000828985,
+     0.003915955,     0.003915955,     0.003915955,
+    -0.140659022,    -0.140659022,    -0.140659022,
+     0.006111833,     0.006111833,     0.006111833,
+
+    -0.003257371,    -0.003257371,    -0.003257371, // (0-2,0,3)
+    -0.012890896,    -0.012890896,    -0.012890896,
+    -0.021780184,    -0.021780184,    -0.021780184,
+     0.016465486,     0.016465486,     0.016465486,
+     0.139385123,     0.139385123,     0.139385123,
+     0.002945519,     0.002945519,     0.002945519,
+     0.000291555,     0.000291555,     0.000291555,
+     0.025585789,     0.025585789,     0.025585789,
+     0.006272205,     0.006272205,     0.006272205,
+    -0.002388125,    -0.002388125,    -0.002388125,
+    -0.002249996,    -0.002249996,    -0.002249996,
+    -0.004658705,    -0.004658705,    -0.004658705,
+    -0.009508024,    -0.009508024,    -0.009508024,
+    -0.156869719,    -0.156869719,    -0.156869719,
+     0.000302461,     0.000302461,     0.000302461,
+     0.005672318,     0.005672318,     0.005672318,
+
+     0.013425194,     0.013425194,     0.013425194, // (0-2,0,4)
+    -0.005822358,    -0.005822358,    -0.005822358,
+    -0.002234564,    -0.002234564,    -0.002234564,
+    -0.006363872,    -0.006363872,    -0.006363872,
+    -0.001465590,    -0.001465590,    -0.001465590,
+    -0.273912363,    -0.273912363,    -0.273912363,
+    -0.012920615,    -0.012920615,    -0.012920615,
+    -0.002148003,    -0.002148003,    -0.002148003,
+     0.005297600,     0.005297600,     0.005297600,
+    -0.015193613,    -0.015193613,    -0.015193613,
+    -0.002463893,    -0.002463893,    -0.002463893,
+     0.052157319,     0.052157319,     0.052157319,
+    -0.000202900,    -0.000202900,    -0.000202900,
+     0.005540465,     0.005540465,     0.005540465,
+     0.295144040,     0.295144040,     0.295144040,
+    -0.028845350,    -0.028845350,    -0.028845350,
+
+    -0.006566304,    -0.006566304,    -0.006566304, // (0-2,0,5)
+    -0.001105205,    -0.001105205,    -0.001105205,
+    -0.001685836,    -0.001685836,    -0.001685836,
+    -0.003194386,    -0.003194386,    -0.003194386,
+    -0.001199854,    -0.001199854,    -0.001199854,
+     0.015421924,     0.015421924,     0.015421924,
+     0.115192196,     0.115192196,     0.115192196,
+    -0.000231987,    -0.000231987,    -0.000231987,
+     0.008316320,     0.008316320,     0.008316320,
+    -0.016758746,    -0.016758746,    -0.016758746,
+    -0.002445727,    -0.002445727,    -0.002445727,
+     0.080670024,     0.080670024,     0.080670024,
+    -0.003115501,    -0.003115501,    -0.003115501,
+     0.002587791,     0.002587791,     0.002587791,
+    -0.010221708,    -0.010221708,    -0.010221708,
+    -0.177149646,    -0.177149646,    -0.177149646,
+
+    -0.005079660,    -0.005079660,    -0.005079660, // (0-2,0,6)
+    -0.011925192,    -0.011925192,    -0.011925192,
+    -0.014270148,    -0.014270148,    -0.014270148,
+     0.007805712,     0.007805712,     0.007805712,
+    -0.013772958,    -0.013772958,    -0.013772958,
+     0.004219976,     0.004219976,     0.004219976,
+    -0.000954442,    -0.000954442,    -0.000954442,
+    -0.255694347,    -0.255694347,    -0.255694347,
+     0.015893669,     0.015893669,     0.015893669,
+    -0.005276456,    -0.005276456,    -0.005276456,
+    -0.006550628,    -0.006550628,    -0.006550628,
+     0.000240351,     0.000240351,     0.000240351,
+    -0.008600808,    -0.008600808,    -0.008600808,
+     0.286812002,     0.286812002,     0.286812002,
+     0.000429281,     0.000429281,     0.000429281,
+     0.006868732,     0.006868732,     0.006868732,
+
+    -0.005224742,    -0.005224742,    -0.005224742, // (0-2,0,7)
+     0.006672921,     0.006672921,     0.006672921,
+    -0.003714513,    -0.003714513,    -0.003714513,
+    -0.007324966,    -0.007324966,    -0.007324966,
+    -0.002267798,    -0.002267798,    -0.002267798,
+     0.006697736,     0.006697736,     0.006697736,
+    -0.006381599,    -0.006381599,    -0.006381599,
+     0.027388953,     0.027388953,     0.027388953,
+    -0.053364033,    -0.053364033,    -0.053364033,
+    -0.006528629,    -0.006528629,    -0.006528629,
+     0.005615009,     0.005615009,     0.005615009,
+     0.024071156,     0.024071156,     0.024071156,
+    -0.011923458,    -0.011923458,    -0.011923458,
+    -0.014965471,    -0.014965471,    -0.014965471,
+    -0.003707340,    -0.003707340,    -0.003707340,
+    -0.001355106,    -0.001355106,    -0.001355106,
+
+     0.041087138,     0.041087138,     0.041087138, // (0-2,0,8)
+    -0.023177438,    -0.023177438,    -0.023177438,
+    -0.002370227,    -0.002370227,    -0.002370227,
+     0.000900653,     0.000900653,     0.000900653,
+    -0.003519113,    -0.003519113,    -0.003519113,
+     0.047562942,     0.047562942,     0.047562942,
+    -0.025927233,    -0.025927233,    -0.025927233,
+     0.015878930,     0.015878930,     0.015878930,
+     0.016866294,     0.016866294,     0.016866294,
+     0.115275958,     0.115275958,     0.115275958,
+    -0.005375032,    -0.005375032,    -0.005375032,
+     0.070287906,     0.070287906,     0.070287906,
+    -0.007716735,    -0.007716735,    -0.007716735,
+    -0.007394883,    -0.007394883,    -0.007394883,
+    -0.085831905,    -0.085831905,    -0.085831905,
+    -0.107035018,    -0.107035018,    -0.107035018,
+
+     0.001574902,     0.001574902,     0.001574902, // (0-2,0,9)
+     0.006476725,     0.006476725,     0.006476725,
+    -0.004368797,    -0.004368797,    -0.004368797,
+    -0.005996139,    -0.005996139,    -0.005996139,
+    -0.001961175,    -0.001961175,    -0.001961175,
+     0.001843756,     0.001843756,     0.001843756,
+    -0.000795551,    -0.000795551,    -0.000795551,
+     0.017728537,     0.017728537,     0.017728537,
+    -0.046142506,    -0.046142506,    -0.046142506,
+    -0.000581687,    -0.000581687,    -0.000581687,
+     0.062923124,     0.062923124,     0.062923124,
+     0.001159794,     0.001159794,     0.001159794,
+    -0.006605954,    -0.006605954,    -0.006605954,
+    -0.010304509,    -0.010304509,    -0.010304509,
+    -0.001999150,    -0.001999150,    -0.001999150,
+     0.004188529,     0.004188529,     0.004188529,
+
+    -0.015564996,    -0.015564996,    -0.015564996, // (0-2,0,10)
+    -0.000042665,    -0.000042665,    -0.000042665,
+    -0.001554026,    -0.001554026,    -0.001554026,
+    -0.004663795,    -0.004663795,    -0.004663795,
+    -0.001018132,    -0.001018132,    -0.001018132,
+     0.030833945,     0.030833945,     0.030833945,
+    -0.033435117,    -0.033435117,    -0.033435117,
+    -0.000226070,    -0.000226070,    -0.000226070,
+     0.010605614,     0.010605614,     0.010605614,
+    -0.011962138,    -0.011962138,    -0.011962138,
+    -0.003437893,    -0.003437893,    -0.003437893,
+    -0.303378514,    -0.303378514,    -0.303378514,
+    -0.002764909,    -0.002764909,    -0.002764909,
+     0.002696539,     0.002696539,     0.002696539,
+    -0.025076301,    -0.025076301,    -0.025076301,
+     0.353690658,     0.353690658,     0.353690658,
+
+    -0.010267196,    -0.010267196,    -0.010267196, // (0-2,0,11)
+     0.000709221,     0.000709221,     0.000709221,
+    -0.015230567,    -0.015230567,    -0.015230567,
+    -0.000187393,    -0.000187393,    -0.000187393,
+    -0.014756803,    -0.014756803,    -0.014756803,
+     0.004259217,     0.004259217,     0.004259217,
+    -0.005529434,    -0.005529434,    -0.005529434,
+     0.020828686,     0.020828686,     0.020828686,
+     0.015871064,     0.015871064,     0.015871064,
+    -0.009486153,    -0.009486153,    -0.009486153,
+    -0.005849931,    -0.005849931,    -0.005849931,
+     0.008209615,     0.008209615,     0.008209615,
+     0.066306147,     0.066306147,     0.066306147,
+     0.033531397,     0.033531397,     0.033531397,
+    -0.000569532,    -0.000569532,    -0.000569532,
+     0.009270397,     0.009270397,     0.009270397,
+
+    -0.107375933,    -0.107375933,    -0.107375933, // (0-2,0,12)
+    -0.013475484,    -0.013475484,    -0.013475484,
+    -0.031861757,    -0.031861757,    -0.031861757,
+     0.012619124,     0.012619124,     0.012619124,
+    -0.031773440,    -0.031773440,    -0.031773440,
+     0.006559095,     0.006559095,     0.006559095,
+    -0.001108022,    -0.001108022,    -0.001108022,
+     0.043803544,     0.043803544,     0.043803544,
+     0.017342703,     0.017342703,     0.017342703,
+    -0.005718482,    -0.005718482,    -0.005718482,
+    -0.005896076,    -0.005896076,    -0.005896076,
+    -0.001850527,    -0.001850527,    -0.001850527,
+    -0.032412625,    -0.032412625,    -0.032412625,
+     0.016598570,     0.016598570,     0.016598570,
+     0.000011074,     0.000011074,     0.000011074,
+     0.009770531,     0.009770531,     0.009770531,
+
+     0.017391771,     0.017391771,     0.017391771, // (0-2,0,13)
+    -0.024588804,    -0.024588804,    -0.024588804,
+    -0.001995858,    -0.001995858,    -0.001995858,
+     0.031076870,     0.031076870,     0.031076870,
+    -0.002307055,    -0.002307055,    -0.002307055,
+    -0.027399437,    -0.027399437,    -0.027399437,
+    -0.010985254,    -0.010985254,    -0.010985254,
+     0.003758121,     0.003758121,     0.003758121,
+     0.007718910,     0.007718910,     0.007718910,
+    -0.002458795,    -0.002458795,    -0.002458795,
+    -0.002719305,    -0.002719305,    -0.002719305,
+     0.034241928,     0.034241928,     0.034241928,
+    -0.002471401,    -0.002471401,    -0.002471401,
+     0.000681513,     0.000681513,     0.000681513,
+     0.012984889,     0.012984889,     0.012984889,
+    -0.010239657,    -0.010239657,    -0.010239657,
+
+    -0.005296665,    -0.005296665,    -0.005296665, // (0-2,0,14)
+    -0.004941121,    -0.004941121,    -0.004941121,
+    -0.001760473,    -0.001760473,    -0.001760473,
+    -0.003270176,    -0.003270176,    -0.003270176,
+    -0.001497442,    -0.001497442,    -0.001497442,
+     0.033680126,     0.033680126,     0.033680126,
+     0.009405484,     0.009405484,     0.009405484,
+     0.005154229,     0.005154229,     0.005154229,
+     0.017761506,     0.017761506,     0.017761506,
+     0.050065382,     0.050065382,     0.050065382,
+    -0.005866127,    -0.005866127,    -0.005866127,
+    -0.026835990,    -0.026835990,    -0.026835990,
+    -0.006425759,    -0.006425759,    -0.006425759,
+    -0.001463026,    -0.001463026,    -0.001463026,
+    -0.024222777,    -0.024222777,    -0.024222777,
+    -0.043074035,    -0.043074035,    -0.043074035,
+
+     0.003290198,     0.003290198,     0.003290198, // (0-2,0,15)
+    -0.006905825,    -0.006905825,    -0.006905825,
+    -0.007382138,    -0.007382138,    -0.007382138,
+    -0.001141505,    -0.001141505,    -0.001141505,
+    -0.007479111,    -0.007479111,    -0.007479111,
+     0.020143024,     0.020143024,     0.020143024,
+    -0.016599679,    -0.016599679,    -0.016599679,
+     0.027678204,     0.027678204,     0.027678204,
+    -0.047607746,    -0.047607746,    -0.047607746,
+    -0.030055519,    -0.030055519,    -0.030055519,
+    -0.017472588,    -0.017472588,    -0.017472588,
+     0.046207911,     0.046207911,     0.046207911,
+     0.042953788,     0.042953788,     0.042953788,
+    -0.013474832,    -0.013474832,    -0.013474832,
+    -0.010334008,    -0.010334008,    -0.010334008,
+    -0.039345679,    -0.039345679,    -0.039345679,
+};
+
+static const double qvec_amalgam_reference[]{
+  0.032868508,
+  0.072854108,
+  0.075831106,
+  0.110432030,
+  0.096887982,
+  0.094833734,
+  0.201333715,
+ -0.080608908,
+  0.068731821,
+  0.047872189,
+ -0.123961630,
+ -0.292706072,
+  0.127453659,
+ -0.069640537,
+  0.102473385,
+ -0.123128698,
+ -0.072745013,
+ -0.056184508,
+ -0.126128189,
+  0.006018975,
+  0.082870682,
+  0.035737788,
+ -0.010427644,
+  0.032408281,
+ -0.037079002,
+ -0.005864535,
+  0.032073163,
+  0.024264188,
+  0.425432489,
+  0.107423386,
+  0.248688911,
+ -0.297062332,
+  0.020683455,
+ -0.096287867,
+ -0.126368074,
+ -0.086839237,
+  0.075028379,
+ -0.121742062,
+  0.032952950,
+  0.031287018,
+  0.039686807,
+ -0.228757265,
+  0.113827034,
+  0.017100373,
+ -0.237251612,
+  0.120001511,
+  0.115043736,
+ -0.064053710,
+  0.010428225,
+  0.008714721,
+  0.019227400,
+ -0.043394797,
+ -0.083676541,
+ -0.090495470,
+ -0.058955600,
+ -0.097112410
+};
+
 extern int test_multi();
 
 #endif  // TEST_MULTICHARGE_H
\ No newline at end of file