diff --git a/include/dftd_cblas.h b/include/dftd_cblas.h
index 840d64b..2188575 100644
--- a/include/dftd_cblas.h
+++ b/include/dftd_cblas.h
@@ -24,7 +24,7 @@ extern "C" {
 #include "lapacke.h"
 }
 
-namespace dftd {
+namespace dftd4 {
 
 inline int BLAS_Add_Mat_x_Vec(TVector<double>& C, TMatrix<double>& A,
                               TVector<double>& V, bool Transpose,
@@ -167,4 +167,4 @@ inline int BLAS_InvertMatrix(TMatrix<double>& a) {
   return EXIT_SUCCESS;
 };
 
-}  // namespace dftd
+}  // namespace dftd4
diff --git a/include/dftd_cutoff.h b/include/dftd_cutoff.h
index d5f78eb..7321e46 100644
--- a/include/dftd_cutoff.h
+++ b/include/dftd_cutoff.h
@@ -17,7 +17,7 @@
  */
 #pragma once
 
-namespace dftd {
+namespace dftd4 {
 
 // Collection of real space cutoffs.
 class TCutoff {
@@ -36,4 +36,4 @@ class TCutoff {
 
 };
 
-}; // namespace dftd
\ No newline at end of file
+}; // namespace dftd4
\ No newline at end of file
diff --git a/include/dftd_damping.h b/include/dftd_damping.h
index 30f76be..f473ce8 100644
--- a/include/dftd_damping.h
+++ b/include/dftd_damping.h
@@ -26,7 +26,7 @@
 
 #include "dftd_dispersion.h"
 
-namespace dftd {
+namespace dftd4 {
 
 /**
  * @brief Collect the D4 parameters.
@@ -38,7 +38,7 @@ namespace dftd {
  */
 extern int d4par(
   const std::string func,
-  dftd::dparam &par,
+  dftd4::dparam &par,
   const bool latm = true
 );
 
diff --git a/include/dftd_dispersion.h b/include/dftd_dispersion.h
index 26b178e..f33736d 100644
--- a/include/dftd_dispersion.h
+++ b/include/dftd_dispersion.h
@@ -25,11 +25,7 @@
 #include "dftd_geometry.h"
 #include "dftd_matrix.h"
 
-namespace dftd {
-
-extern inline double fdmpr_bj(const int n, const double r, const double c);
-extern inline double fdmprdr_bj(const int n, const double r, const double c);
-
+namespace dftd4 {
 
 class dparam {
  public:
@@ -42,6 +38,26 @@ class dparam {
   int alp;
 };
 
+/**
+ * @brief Wrapper to handle the evaluation of dispersion energy and derivatives.
+ * 
+ * @param Dat Molecular information (only for consistency with ORCA).
+ * @param mol Molecular geometry.
+ * @param par DFT-D4 parameters.
+ * @param cutoff Real-space cutoffs for CN and dispersion.
+ * @param energy Dispersion energy.
+ * @param GRAD Dispersion gradient.
+ * @return Exit status.
+ */
+extern int DFTVDW_D4(
+  const TMolInfo &Dat,
+  const TMolecule &mol,
+  const dparam &par,
+  TCutoff cutoff,
+  double &energy,
+  double *GRAD
+);
+
 // Generic wrappers for two- and three-body dispersion
 
 extern int get_dispersion2(
@@ -150,19 +166,7 @@ extern int get_atm_dispersion_derivs(
 
 extern double triple_scale(int ii, int jj, int kk);
 
-/**
- * @brief Wrapper to handle the evaluation of dispersion energy and derivatives.
- * 
- * @param mol Molecular geometry.
- * @param par DFT-D4 parameters.
- * @param charge Charge of the molecule.
- * @param cutoff Real-space cutoffs for CN and dispersion.
- * @param energy Dispersion energy.
- * @param GRAD Dispersion gradient.
- * @return Exit status.
- */
-extern
-int get_dispersion(const TMolecule &mol, const dparam &par, const int &charge,
-                   TCutoff cutoff, double &energy, double *GRAD);
+extern inline double fdmpr_bj(const int n, const double r, const double c);
+extern inline double fdmprdr_bj(const int n, const double r, const double c);
 
-}  // namespace dftd
+}  // namespace dftd4
diff --git a/include/dftd_eeq.h b/include/dftd_eeq.h
index 4b0e0e1..77e2fff 100644
--- a/include/dftd_eeq.h
+++ b/include/dftd_eeq.h
@@ -25,7 +25,7 @@
 #include "dftd_geometry.h"
 #include "dftd_matrix.h"
 
-namespace dftd {
+namespace dftd4 {
 
 extern int get_charges(
   const TMolecule& mol,
diff --git a/include/dftd_geometry.h b/include/dftd_geometry.h
index 347f519..be69a1e 100644
--- a/include/dftd_geometry.h
+++ b/include/dftd_geometry.h
@@ -20,7 +20,7 @@
 
 #include "dftd_matrix.h"
 
-namespace dftd {
+namespace dftd4 {
 // Input of the molecular geometry
 class TMolecule {
  public:
@@ -44,4 +44,13 @@ class TMolecule {
   }
 };
 
-}  // namespace dftd
+class TMolInfo {
+  public:
+    int Charge;
+
+    TMolInfo(int c = 0) {
+      Charge = c;
+    }
+};
+
+}  // namespace dftd4
diff --git a/include/dftd_matrix.h b/include/dftd_matrix.h
index 62904a3..1ae7f9e 100644
--- a/include/dftd_matrix.h
+++ b/include/dftd_matrix.h
@@ -22,7 +22,7 @@
 #include <cstring>
 #include <iostream>
 
-namespace dftd {
+namespace dftd4 {
 
 // Define a vector
 template <class T>
@@ -40,7 +40,7 @@ class TVector {
   ~TVector() {
     if (p != 0) Delete();
   }
-  void New(int VectorLength) {
+  void NewVector(int VectorLength) {
     if (VectorLength < 0) {
       std::exit(EXIT_FAILURE);
     }
@@ -60,6 +60,11 @@ class TVector {
       Init();
     }
   }
+
+  // alias for NewVector
+  void New(int VectorLength) { return NewVector(VectorLength); }
+  void NewVec(int VectorLength) { return NewVector(VectorLength); }
+
   void Delete(void) {
     if (p != 0 && N != 0) {
       delete[] p;
@@ -116,7 +121,8 @@ class TMatrix {
     if (p != 0) Delete();
   }
 
-  void New(int r, int c) {
+
+  void NewMatrix(int r, int c) {
     if (r < 0 || c < 0) std::exit(EXIT_FAILURE);
     if (p != 0 && r == rows && c == cols) {
       Init();
@@ -135,7 +141,11 @@ class TMatrix {
     return;
   }
 
-  void New(const TMatrix& v) { New(v.rows, v.cols); }
+  void NewMatrix(const TMatrix& v) { NewMatrix(v.rows, v.cols); }
+
+  // alias for NewMatrix
+  void New(int r, int c) { return NewMatrix(r, c); }
+  void NewMat(int r, int c) { return NewMatrix(r, c); }
 
   void Delete(void) {
     if (p != 0 && rows * cols != 0) {
@@ -172,7 +182,7 @@ class TMatrix {
         // for non-square matrix, we need an additional copy
         TMatrix<T> temp;
         temp.CopyMat(*this);
-        New(cols, rows);
+        NewMatrix(cols, rows);
         for (i = 0; i < rows; i++) {
           for (j = 0; j < cols; j++) {
             p[i * cols + j] = temp.p[j * cols + i];
@@ -214,4 +224,4 @@ class TMatrix {
   inline T* operator[](int i) { return p + i * cols; }
 };
 
-}  // namespace dftd
+}  // namespace dftd4
diff --git a/include/dftd_model.h b/include/dftd_model.h
index 4e5a2d2..7782928 100644
--- a/include/dftd_model.h
+++ b/include/dftd_model.h
@@ -21,7 +21,7 @@
 #include "dftd_geometry.h"
 #include "dftd_matrix.h"
 
-namespace dftd {
+namespace dftd4 {
 
 extern inline double trapzd(const double a[23], const double b[23]);
 
@@ -32,6 +32,8 @@ extern inline double zeta(const double a, const double c, const double qref,
 extern inline double dzeta(const double a, const double c, const double qref,
                     const double qmod);
 
+extern int get_max_ref(const TMolecule &mol, int &mref);
+
 extern int weight_references(
   const TMolecule& mol,
   const TVector<double>& cn,
@@ -57,4 +59,4 @@ extern int set_refalpha_eeq(const TMolecule& mol, TMatrix<double>& alpha);
 
 extern bool is_exceptional(double val);
 
-} // namespace dftd
\ No newline at end of file
+} // namespace dftd4
\ No newline at end of file
diff --git a/include/dftd_ncoord.h b/include/dftd_ncoord.h
index 1c79d41..02dc26a 100644
--- a/include/dftd_ncoord.h
+++ b/include/dftd_ncoord.h
@@ -26,10 +26,10 @@
 #include "dftd_geometry.h"
 #include "dftd_matrix.h"
 
-namespace dftd {
+namespace dftd4 {
 
 /**
- * Calculate all distance pairs and store in matrix
+ * Calculate all distance pairs and store in matrix.
  * 
  * @param mol Molecule object.
  * @param dist Distance matrix (inout).
@@ -186,4 +186,4 @@ extern inline double log_cn_cut(const double cn_max, const double cn);
 
 extern inline double dlog_cn_cut(const double cn_max, const double cn);
 
-};  // namespace dftd
+};  // namespace dftd4
diff --git a/include/dftd_parameters.h b/include/dftd_parameters.h
index baf1ba3..3da1aa7 100644
--- a/include/dftd_parameters.h
+++ b/include/dftd_parameters.h
@@ -26,7 +26,7 @@
  */
 #pragma once
 
-namespace dftd
+namespace dftd4
 {
 
 static const double zeff[119] { 0,
diff --git a/include/dftd_readxyz.h b/include/dftd_readxyz.h
index 4d5b8c2..03f8309 100644
--- a/include/dftd_readxyz.h
+++ b/include/dftd_readxyz.h
@@ -21,6 +21,6 @@
 
 #include "dftd_geometry.h"
 
-extern void read_xyzfile(const std::string&, dftd::TMolecule&);
+extern void read_xyzfile(const std::string&, dftd4::TMolecule&);
 
 extern int element(const std::string&);
diff --git a/src/dftd_cutoff.cpp b/src/dftd_cutoff.cpp
index a7cf33c..c388233 100644
--- a/src/dftd_cutoff.cpp
+++ b/src/dftd_cutoff.cpp
@@ -17,7 +17,7 @@
  */
 #include "dftd_cutoff.h"
 
-namespace dftd {
+namespace dftd4 {
 
 // Real space cutoff for CN within D4
 static const double cn_default = 30.0;
@@ -47,4 +47,4 @@ void TCutoff::disable() {
   cn_eeq = new_cutoff;
 };
 
-}; // namespace dftd
\ No newline at end of file
+}; // namespace dftd4
\ No newline at end of file
diff --git a/src/dftd_damping.cpp b/src/dftd_damping.cpp
index 60e325b..5256560 100644
--- a/src/dftd_damping.cpp
+++ b/src/dftd_damping.cpp
@@ -22,7 +22,7 @@
 #include "dftd_damping.h"
 #include "dftd_dispersion.h"
 
-namespace dftd {
+namespace dftd4 {
 
 enum dfunc {
   none = 0,
@@ -1536,7 +1536,7 @@ dparam get_d4eeqbjmbd_2019_parameter(dfunc num) {
 
 int d4par(
   const std::string func,
-  dftd::dparam &par,
+  dftd4::dparam &par,
   const bool latm/* = true*/
 ) {
   auto num = get_dfunc(func);
@@ -1548,4 +1548,4 @@ int d4par(
   return EXIT_SUCCESS;
 }
 
-}  // namespace dftd
+}  // namespace dftd4
diff --git a/src/dftd_dispersion.cpp b/src/dftd_dispersion.cpp
index 621c686..3d17063 100644
--- a/src/dftd_dispersion.cpp
+++ b/src/dftd_dispersion.cpp
@@ -31,7 +31,7 @@
 #include "dftd_ncoord.h"
 #include "dftd_parameters.h"
 
-namespace dftd {
+namespace dftd4 {
 
 inline double fdmpr_bj(const int n, const double r, const double c) {
   return 1.0 / (pow(r, n) + pow(c, n));
@@ -486,8 +486,14 @@ int get_atm_dispersion_derivs(
 }
 
 
-int get_dispersion(const TMolecule &mol, const dparam &par, const int &charge,
-                   TCutoff cutoff, double &energy, double *GRAD) {
+int DFTVDW_D4(
+  const TMolInfo &Dat,
+  const TMolecule &mol,
+  const dparam &par,
+  TCutoff cutoff,
+  double &energy,
+  double *GRAD
+) {
   // setup variables
   bool lmbd = (par.s9 != 0.0);
   bool lgrad = !!GRAD;
@@ -495,7 +501,7 @@ int get_dispersion(const TMolecule &mol, const dparam &par, const int &charge,
 
   // distances
   TMatrix<double> dist;
-  dist.New(mol.NAtoms, mol.NAtoms);
+  dist.NewMat(mol.NAtoms, mol.NAtoms);
   calc_distances(mol, dist);
 
   TVector<double> cn;        // D4 coordination number
@@ -504,16 +510,16 @@ int get_dispersion(const TMolecule &mol, const dparam &par, const int &charge,
   TMatrix<double> dqdr;      // derivative of partial charges
   TVector<double> gradient;  // derivative of dispersion energy
 
-  cn.New(mol.NAtoms);
-  q.New(mol.NAtoms);
+  cn.NewVec(mol.NAtoms);
+  q.NewVec(mol.NAtoms);
   if (lgrad) {
-    dcndr.New(3 * mol.NAtoms, mol.NAtoms);
-    dqdr.New(3 * mol.NAtoms, mol.NAtoms);
-    gradient.New(3 * mol.NAtoms);
+    dcndr.NewMat(3 * mol.NAtoms, mol.NAtoms);
+    dqdr.NewMat(3 * mol.NAtoms, mol.NAtoms);
+    gradient.NewVec(3 * mol.NAtoms);
   }
 
-  // get charges from EEQ model
-  info = get_charges(mol, dist, charge, cutoff.cn_eeq, q, dqdr, lgrad);
+  // calculate partial charges from EEQ model
+  info = get_charges(mol, dist, Dat.Charge, cutoff.cn_eeq, q, dqdr, lgrad);
   if (!info == EXIT_SUCCESS) return info;
 
   // get the D4 coordination number
@@ -522,19 +528,17 @@ int get_dispersion(const TMolecule &mol, const dparam &par, const int &charge,
 
 
   // maximum number of reference systems
-  int mref = refn[mol.at(0)];
-  for (int i = 1; i != mol.NAtoms; i++) {
-    int val = refn[mol.at(i)];
-    if (val > mref) mref = val;
-  }
+  int mref{0};
+  info = get_max_ref(mol, mref); 
+  if (!info == EXIT_SUCCESS) return info;
 
   TMatrix<double> gwvec;
   TMatrix<double> dgwdcn;
   TMatrix<double> dgwdq;
-  gwvec.New(mref, mol.NAtoms);
+  gwvec.NewMat(mref, mol.NAtoms);
   if (lgrad) {
-    dgwdcn.New(mref, mol.NAtoms);
-    dgwdq.New(mref, mol.NAtoms);    
+    dgwdcn.NewMat(mref, mol.NAtoms);
+    dgwdq.NewMat(mref, mol.NAtoms);    
   }
   info = weight_references(mol, cn, q, gwvec, dgwdcn, dgwdq, lgrad);
   if (!info == EXIT_SUCCESS) return info;
@@ -542,10 +546,10 @@ int get_dispersion(const TMolecule &mol, const dparam &par, const int &charge,
   TMatrix<double> c6;
   TMatrix<double> dc6dcn;
   TMatrix<double> dc6dq;
-  c6.New(mol.NAtoms, mol.NAtoms);
+  c6.NewMat(mol.NAtoms, mol.NAtoms);
   if (lgrad) {
-    dc6dcn.New(mol.NAtoms, mol.NAtoms);
-    dc6dq.New(mol.NAtoms, mol.NAtoms);
+    dc6dcn.NewMat(mol.NAtoms, mol.NAtoms);
+    dc6dq.NewMat(mol.NAtoms, mol.NAtoms);
   }
 
   info = get_atomic_c6(mol, gwvec, dgwdcn, dgwdq, c6, dc6dcn, dc6dq, lgrad);
@@ -558,10 +562,10 @@ int get_dispersion(const TMolecule &mol, const dparam &par, const int &charge,
   TVector<double> dEdcn;
   TVector<double> dEdq;
   TVector<double> energies;
-  energies.New(mol.NAtoms);
+  energies.NewVec(mol.NAtoms);
   if (lgrad) {
-    dEdcn.New(mol.NAtoms);
-    dEdq.New(mol.NAtoms);    
+    dEdcn.NewVec(mol.NAtoms);
+    dEdq.NewVec(mol.NAtoms);    
   }
 
   info = get_dispersion2(
@@ -588,10 +592,10 @@ int get_dispersion(const TMolecule &mol, const dparam &par, const int &charge,
     }
 
     // calculate weight references
-    gwvec.New(mref, mol.NAtoms);
+    gwvec.NewMat(mref, mol.NAtoms);
     if (lgrad) {
-      dgwdcn.New(mref, mol.NAtoms);
-      dgwdq.New(mref, mol.NAtoms);    
+      dgwdcn.NewMat(mref, mol.NAtoms);
+      dgwdq.NewMat(mref, mol.NAtoms);    
     }
     info = weight_references(mol, cn, q, gwvec, dgwdcn, dgwdq, lgrad);
     if (!info == EXIT_SUCCESS) return info;
@@ -600,10 +604,10 @@ int get_dispersion(const TMolecule &mol, const dparam &par, const int &charge,
     q.Delete();
 
     // calculate reference C6 coefficients
-    c6.New(mol.NAtoms, mol.NAtoms);
+    c6.NewMat(mol.NAtoms, mol.NAtoms);
     if (lgrad) {
-      dc6dcn.New(mol.NAtoms, mol.NAtoms);
-      dc6dq.New(mol.NAtoms, mol.NAtoms);
+      dc6dcn.NewMat(mol.NAtoms, mol.NAtoms);
+      dc6dq.NewMat(mol.NAtoms, mol.NAtoms);
     }
     info = get_atomic_c6(mol, gwvec, dgwdcn, dgwdq, c6, dc6dcn, dc6dq, lgrad);
     if (!info == EXIT_SUCCESS) return info;
@@ -657,4 +661,4 @@ int get_dispersion(const TMolecule &mol, const dparam &par, const int &charge,
   return EXIT_SUCCESS;
 }
 
-}  // namespace dftd
+}  // namespace dftd4
diff --git a/src/dftd_eeq.cpp b/src/dftd_eeq.cpp
index 8baaacb..8bc8a1a 100644
--- a/src/dftd_eeq.cpp
+++ b/src/dftd_eeq.cpp
@@ -29,84 +29,87 @@
 #include "dftd_ncoord.h"
 
 // wrap everything in the dftd namespace to keep it nicely confined
-namespace dftd {
+namespace dftd4 {
 
 static const double xi[87]{
-    0.0,        1.23695041, 1.26590957, 0.54341808, 0.99666991, 1.26691604,
-    1.40028282, 1.55819364, 1.56866440, 1.57540015, 1.15056627, 0.55936220,
-    0.72373742, 1.12910844, 1.12306840, 1.52672442, 1.40768172, 1.48154584,
-    1.31062963, 0.40374140, 0.75442607, 0.76482096, 0.98457281, 0.96702598,
-    1.05266584, 0.93274875, 1.04025281, 0.92738624, 1.07419210, 1.07900668,
-    1.04712861, 1.15018618, 1.15388455, 1.36313743, 1.36485106, 1.39801837,
-    1.18695346, 0.36273870, 0.58797255, 0.71961946, 0.96158233, 0.89585296,
-    0.81360499, 1.00794665, 0.92613682, 1.09152285, 1.14907070, 1.13508911,
-    1.08853785, 1.11005982, 1.12452195, 1.21642129, 1.36507125, 1.40340000,
-    1.16653482, 0.34125098, 0.58884173, 0.68441115, 0.56999999, 0.56999999,
-    0.56999999, 0.56999999, 0.56999999, 0.56999999, 0.56999999, 0.56999999,
-    0.56999999, 0.56999999, 0.56999999, 0.56999999, 0.56999999, 0.56999999,
-    0.87936784, 1.02761808, 0.93297476, 1.10172128, 0.97350071, 1.16695666,
-    1.23997927, 1.18464453, 1.14191734, 1.12334192, 1.01485321, 1.12950808,
-    1.30804834, 1.33689961, 1.27465977};
+  0.0,        1.23695041, 1.26590957, 0.54341808, 0.99666991, 1.26691604,
+  1.40028282, 1.55819364, 1.56866440, 1.57540015, 1.15056627, 0.55936220,
+  0.72373742, 1.12910844, 1.12306840, 1.52672442, 1.40768172, 1.48154584,
+  1.31062963, 0.40374140, 0.75442607, 0.76482096, 0.98457281, 0.96702598,
+  1.05266584, 0.93274875, 1.04025281, 0.92738624, 1.07419210, 1.07900668,
+  1.04712861, 1.15018618, 1.15388455, 1.36313743, 1.36485106, 1.39801837,
+  1.18695346, 0.36273870, 0.58797255, 0.71961946, 0.96158233, 0.89585296,
+  0.81360499, 1.00794665, 0.92613682, 1.09152285, 1.14907070, 1.13508911,
+  1.08853785, 1.11005982, 1.12452195, 1.21642129, 1.36507125, 1.40340000,
+  1.16653482, 0.34125098, 0.58884173, 0.68441115, 0.56999999, 0.56999999,
+  0.56999999, 0.56999999, 0.56999999, 0.56999999, 0.56999999, 0.56999999,
+  0.56999999, 0.56999999, 0.56999999, 0.56999999, 0.56999999, 0.56999999,
+  0.87936784, 1.02761808, 0.93297476, 1.10172128, 0.97350071, 1.16695666,
+  1.23997927, 1.18464453, 1.14191734, 1.12334192, 1.01485321, 1.12950808,
+  1.30804834, 1.33689961, 1.27465977
+};
 static const double gam[87]{
-    0.0,         -0.35015861, 1.04121227,  0.09281243,  0.09412380,
-    0.26629137,  0.19408787,  0.05317918,  0.03151644,  0.32275132,
-    1.30996037,  0.24206510,  0.04147733,  0.11634126,  0.13155266,
-    0.15350650,  0.15250997,  0.17523529,  0.28774450,  0.42937314,
-    0.01896455,  0.07179178,  -0.01121381, -0.03093370, 0.02716319,
-    -0.01843812, -0.15270393, -0.09192645, -0.13418723, -0.09861139,
-    0.18338109,  0.08299615,  0.11370033,  0.19005278,  0.10980677,
-    0.12327841,  0.25345554,  0.58615231,  0.16093861,  0.04548530,
-    -0.02478645, 0.01909943,  0.01402541,  -0.03595279, 0.01137752,
-    -0.03697213, 0.08009416,  0.02274892,  0.12801822,  -0.02078702,
-    0.05284319,  0.07581190,  0.09663758,  0.09547417,  0.07803344,
-    0.64913257,  0.15348654,  0.05054344,  0.11000000,  0.11000000,
-    0.11000000,  0.11000000,  0.11000000,  0.11000000,  0.11000000,
-    0.11000000,  0.11000000,  0.11000000,  0.11000000,  0.11000000,
-    0.11000000,  0.11000000,  -0.02786741, 0.01057858,  -0.03892226,
-    -0.04574364, -0.03874080, -0.03782372, -0.07046855, 0.09546597,
-    0.21953269,  0.02522348,  0.15263050,  0.08042611,  0.01878626,
-    0.08715453,  0.10500484};
+  0.0,         -0.35015861, 1.04121227,  0.09281243,  0.09412380,
+  0.26629137,  0.19408787,  0.05317918,  0.03151644,  0.32275132,
+  1.30996037,  0.24206510,  0.04147733,  0.11634126,  0.13155266,
+  0.15350650,  0.15250997,  0.17523529,  0.28774450,  0.42937314,
+  0.01896455,  0.07179178,  -0.01121381, -0.03093370, 0.02716319,
+  -0.01843812, -0.15270393, -0.09192645, -0.13418723, -0.09861139,
+  0.18338109,  0.08299615,  0.11370033,  0.19005278,  0.10980677,
+  0.12327841,  0.25345554,  0.58615231,  0.16093861,  0.04548530,
+  -0.02478645, 0.01909943,  0.01402541,  -0.03595279, 0.01137752,
+  -0.03697213, 0.08009416,  0.02274892,  0.12801822,  -0.02078702,
+  0.05284319,  0.07581190,  0.09663758,  0.09547417,  0.07803344,
+  0.64913257,  0.15348654,  0.05054344,  0.11000000,  0.11000000,
+  0.11000000,  0.11000000,  0.11000000,  0.11000000,  0.11000000,
+  0.11000000,  0.11000000,  0.11000000,  0.11000000,  0.11000000,
+  0.11000000,  0.11000000,  -0.02786741, 0.01057858,  -0.03892226,
+  -0.04574364, -0.03874080, -0.03782372, -0.07046855, 0.09546597,
+  0.21953269,  0.02522348,  0.15263050,  0.08042611,  0.01878626,
+  0.08715453,  0.10500484
+};
 static const double kappa[87]{
-    0.0,         0.04916110,  0.10937243,  -0.12349591, -0.02665108,
-    -0.02631658, 0.06005196,  0.09279548,  0.11689703,  0.15704746,
-    0.07987901,  -0.10002962, -0.07712863, -0.02170561, -0.04964052,
-    0.14250599,  0.07126660,  0.13682750,  0.14877121,  -0.10219289,
-    -0.08979338, -0.08273597, -0.01754829, -0.02765460, -0.02558926,
-    -0.08010286, -0.04163215, -0.09369631, -0.03774117, -0.05759708,
-    0.02431998,  -0.01056270, -0.02692862, 0.07657769,  0.06561608,
-    0.08006749,  0.14139200,  -0.05351029, -0.06701705, -0.07377246,
-    -0.02927768, -0.03867291, -0.06929825, -0.04485293, -0.04800824,
-    -0.01484022, 0.07917502,  0.06619243,  0.02434095,  -0.01505548,
-    -0.03030768, 0.01418235,  0.08953411,  0.08967527,  0.07277771,
-    -0.02129476, -0.06188828, -0.06568203, -0.11000000, -0.11000000,
-    -0.11000000, -0.11000000, -0.11000000, -0.11000000, -0.11000000,
-    -0.11000000, -0.11000000, -0.11000000, -0.11000000, -0.11000000,
-    -0.11000000, -0.11000000, -0.03585873, -0.03132400, -0.05902379,
-    -0.02827592, -0.07606260, -0.02123839, 0.03814822,  0.02146834,
-    0.01580538,  -0.00894298, -0.05864876, -0.01817842, 0.07721851,
-    0.07936083,  0.05849285};
+  0.0,         0.04916110,  0.10937243,  -0.12349591, -0.02665108,
+  -0.02631658, 0.06005196,  0.09279548,  0.11689703,  0.15704746,
+  0.07987901,  -0.10002962, -0.07712863, -0.02170561, -0.04964052,
+  0.14250599,  0.07126660,  0.13682750,  0.14877121,  -0.10219289,
+  -0.08979338, -0.08273597, -0.01754829, -0.02765460, -0.02558926,
+  -0.08010286, -0.04163215, -0.09369631, -0.03774117, -0.05759708,
+  0.02431998,  -0.01056270, -0.02692862, 0.07657769,  0.06561608,
+  0.08006749,  0.14139200,  -0.05351029, -0.06701705, -0.07377246,
+  -0.02927768, -0.03867291, -0.06929825, -0.04485293, -0.04800824,
+  -0.01484022, 0.07917502,  0.06619243,  0.02434095,  -0.01505548,
+  -0.03030768, 0.01418235,  0.08953411,  0.08967527,  0.07277771,
+  -0.02129476, -0.06188828, -0.06568203, -0.11000000, -0.11000000,
+  -0.11000000, -0.11000000, -0.11000000, -0.11000000, -0.11000000,
+  -0.11000000, -0.11000000, -0.11000000, -0.11000000, -0.11000000,
+  -0.11000000, -0.11000000, -0.03585873, -0.03132400, -0.05902379,
+  -0.02827592, -0.07606260, -0.02123839, 0.03814822,  0.02146834,
+  0.01580538,  -0.00894298, -0.05864876, -0.01817842, 0.07721851,
+  0.07936083,  0.05849285
+};
 static const double alp[87]{
-    0.0,        0.55159092, 0.66205886, 0.90529132, 1.51710827, 2.86070364,
-    1.88862966, 1.32250290, 1.23166285, 1.77503721, 1.11955204, 1.28263182,
-    1.22344336, 1.70936266, 1.54075036, 1.38200579, 2.18849322, 1.36779065,
-    1.27039703, 1.64466502, 1.58859404, 1.65357953, 1.50021521, 1.30104175,
-    1.46301827, 1.32928147, 1.02766713, 1.02291377, 0.94343886, 1.14881311,
-    1.47080755, 1.76901636, 1.98724061, 2.41244711, 2.26739524, 2.95378999,
-    1.20807752, 1.65941046, 1.62733880, 1.61344972, 1.63220728, 1.60899928,
-    1.43501286, 1.54559205, 1.32663678, 1.37644152, 1.36051851, 1.23395526,
-    1.65734544, 1.53895240, 1.97542736, 1.97636542, 2.05432381, 3.80138135,
-    1.43893803, 1.75505957, 1.59815118, 1.76401732, 1.63999999, 1.63999999,
-    1.63999999, 1.63999999, 1.63999999, 1.63999999, 1.63999999, 1.63999999,
-    1.63999999, 1.63999999, 1.63999999, 1.63999999, 1.63999999, 1.63999999,
-    1.47055223, 1.81127084, 1.40189963, 1.54015481, 1.33721475, 1.57165422,
-    1.04815857, 1.78342098, 2.79106396, 1.78160840, 2.47588882, 2.37670734,
-    1.76613217, 2.66172302, 2.82773085};
+  0.0,        0.55159092, 0.66205886, 0.90529132, 1.51710827, 2.86070364,
+  1.88862966, 1.32250290, 1.23166285, 1.77503721, 1.11955204, 1.28263182,
+  1.22344336, 1.70936266, 1.54075036, 1.38200579, 2.18849322, 1.36779065,
+  1.27039703, 1.64466502, 1.58859404, 1.65357953, 1.50021521, 1.30104175,
+  1.46301827, 1.32928147, 1.02766713, 1.02291377, 0.94343886, 1.14881311,
+  1.47080755, 1.76901636, 1.98724061, 2.41244711, 2.26739524, 2.95378999,
+  1.20807752, 1.65941046, 1.62733880, 1.61344972, 1.63220728, 1.60899928,
+  1.43501286, 1.54559205, 1.32663678, 1.37644152, 1.36051851, 1.23395526,
+  1.65734544, 1.53895240, 1.97542736, 1.97636542, 2.05432381, 3.80138135,
+  1.43893803, 1.75505957, 1.59815118, 1.76401732, 1.63999999, 1.63999999,
+  1.63999999, 1.63999999, 1.63999999, 1.63999999, 1.63999999, 1.63999999,
+  1.63999999, 1.63999999, 1.63999999, 1.63999999, 1.63999999, 1.63999999,
+  1.47055223, 1.81127084, 1.40189963, 1.54015481, 1.33721475, 1.57165422,
+  1.04815857, 1.78342098, 2.79106396, 1.78160840, 2.47588882, 2.37670734,
+  1.76613217, 2.66172302, 2.82773085
+};
 
 static const double small = 1e-14;
 static const double pi = 3.1415926535897932384626433832795029;
 static const double sqrtpi = std::sqrt(pi);
 static const double sqrt2pi = std::sqrt(2.0 / pi);
-static const char _p_up = 'u', _p_lo = 'l';
 
 
 int get_charges(
@@ -124,8 +127,8 @@ int get_charges(
   TVector<double> cn;        // EEQ cordination number
   TMatrix<double> dcndr;     // Derivative of EEQ-CN
 
-  cn.New(mol.NAtoms);
-  if (lgrad) dcndr.New(mol.NAtoms, 3 * mol.NAtoms);
+  cn.NewVec(mol.NAtoms);
+  if (lgrad) dcndr.NewMat(mol.NAtoms, 3 * mol.NAtoms);
 
   // get the EEQ coordination number
   info = get_ncoord_erf(mol, dist, cutoff, cn, dcndr, lgrad);
@@ -269,13 +272,11 @@ int eeq_chrgeq(const TMolecule& mol, const TMatrix<double>& dist,
 
   TMatrix<double> Amat;   // Coulomb matrix
   TVector<double> xvec;   // x (chi) vector
-  Amat.New(m, m);
-  xvec.New(m);
+  Amat.NewMat(m, m);
+  xvec.NewVec(m);
 
   TVector<double> dxdcn;  // Derivative of chi vector w.r.t. CN
-  if (lgrad) {
-    dxdcn.New(m);
-  }
+  if (lgrad) dxdcn.NewVec(m);
 
   info = get_vrhs(mol, charge, cn, xvec, dxdcn, lgrad);
   if (!info == EXIT_SUCCESS) return info;
@@ -284,10 +285,10 @@ int eeq_chrgeq(const TMolecule& mol, const TMatrix<double>& dist,
   if (!info == EXIT_SUCCESS) return info;
 
   TVector<double> vrhs;
-  vrhs.New(m);
+  vrhs.NewVec(m);
 
   TMatrix<double> Ainv;
-  Ainv.New(m, m);
+  Ainv.NewMat(m, m);
   Ainv.CopyMat(Amat);
 
   // solve: A Q = X (Eq.4) -> Q = Ainv X
@@ -337,9 +338,9 @@ int eeq_chrgeq(const TMolecule& mol, const TMatrix<double>& dist,
     int ThreeN = 3 * mol.NAtoms;
 
     TMatrix<double> dAmat;
-    dAmat.New(ThreeN, m);
+    dAmat.NewMat(ThreeN, m);
     TMatrix<double> atrace;
-    atrace.New(m, 3);
+    atrace.NewMat(m, 3);
 
     info = get_damat_0d(mol, dist, vrhs, Amat, dAmat, atrace);
     if (!info == EXIT_SUCCESS) return info;
@@ -362,7 +363,7 @@ int eeq_chrgeq(const TMolecule& mol, const TMatrix<double>& dist,
 
     // Ainv with last column removed
     TMatrix<double> A;
-    A.New(Ainv.rows, Ainv.cols - 1);
+    A.NewMat(Ainv.rows, Ainv.cols - 1);
     for (int i = 0; i < Ainv.rows; i++) {
       for (int j = 0; j < Ainv.cols - 1; j++) {
         A(i, j) = Ainv(i, j);
@@ -383,4 +384,4 @@ int eeq_chrgeq(const TMolecule& mol, const TMatrix<double>& dist,
   return EXIT_SUCCESS;
 }
 
-}  // namespace dftd
+}  // namespace dftd4
diff --git a/src/dftd_model.cpp b/src/dftd_model.cpp
index c15dc5e..88eb9ef 100644
--- a/src/dftd_model.cpp
+++ b/src/dftd_model.cpp
@@ -24,7 +24,7 @@
 #include "dftd_model.h"
 #include "dftd_parameters.h"
 
-namespace dftd {
+namespace dftd4 {
 
 // Parameters
 const double wf = 6.0;
@@ -94,6 +94,17 @@ inline double dzeta(const double a, const double c, const double qref,
   }
 }
 
+int get_max_ref(const TMolecule &mol, int &mref) {
+  mref = refn[mol.at(0)];
+  for (int i = 1; i != mol.NAtoms; i++) {
+    int val = refn[mol.at(i)];
+    if (val > mref) mref = val;
+  }
+
+  if (mref == 0.0) return EXIT_FAILURE;
+  return EXIT_SUCCESS;
+}
+
 int weight_references(
   const TMolecule& mol,
   const TVector<double>& cn,
@@ -103,9 +114,8 @@ int weight_references(
   TMatrix<double>& dgwdq,
   bool lgrad /*= false*/
 ) {
-  int info{0};
-  int iat = 0, jat = 0, izp = 0, jzp = 0;
-  double gw = 0.0, twf = 0.0, c6 = 0.0, maxcn = 0.0, tmp = 0.0;
+  int izp{0};
+  double gw{0.0}, twf{0.0}, maxcn{0.0};
   double norm{0.0}, dnorm{0.0};
   double zi{0.0}, gi{0.0};
   double expw{0.0}, dexpw{0.0}, gwk{0.0}, dgwk{0.0};
@@ -207,18 +217,16 @@ int get_atomic_c6(
   TMatrix<double>& dc6dq,
   bool lgrad /*= false*/
 ) {
-  int iat{0}, izp{0}, jzp{0}, info{0};
-  double iz{0.0}, refc6{0.0}, dc6{0.0};
+  int izp{0}, jzp{0}, info{0};
+  double refc6{0.0}, dc6{0.0};
 
-  // find maximum number of references
-  int mref = refn[mol.at(0)];
-  for (int i = 1; i != mol.NAtoms; i++) {
-    int val = refn[mol.at(i)];
-    if (val > mref) mref = val;
-  }
+  // maximum number of reference systems
+  int mref{0};
+  info = get_max_ref(mol, mref); 
+  if (!info == EXIT_SUCCESS) return info;
 
   TMatrix<double> alpha;
-  alpha.New(mol.NAtoms, 23*mref);
+  alpha.NewMat(mol.NAtoms, 23*mref);
   info = set_refalpha_eeq(mol, alpha);
   if (!info == EXIT_SUCCESS) return info;
 
@@ -307,4 +315,4 @@ bool is_exceptional(double val) {
   return std::isnan(val) || (fabs(val) > std::numeric_limits<double>::max());
 }
 
-} // namespace dftd
\ No newline at end of file
+} // namespace dftd4
\ No newline at end of file
diff --git a/src/dftd_ncoord.cpp b/src/dftd_ncoord.cpp
index 281ff7f..09dcc46 100644
--- a/src/dftd_ncoord.cpp
+++ b/src/dftd_ncoord.cpp
@@ -31,7 +31,7 @@
 #include "dftd_geometry.h"
 #include "dftd_matrix.h"
 
-namespace dftd {
+namespace dftd4 {
 
 // convert bohr (a.u.) to Ångström and back
 static const double autoaa = 0.52917721090449243;
@@ -45,25 +45,25 @@ static const double aatoau = 1.0 / autoaa;
  * Only the scaled values below are used (`rad`).
 */
 static const double covalent_rad_d3[119]{ 0.0,
-    0.32, 0.46,                                                  //  H,He
-    1.20, 0.94, 0.77, 0.75, 0.71, 0.63, 0.64, 0.67,              //  Li-Ne
-    1.40, 1.25, 1.13, 1.04, 1.10, 1.02, 0.99, 0.96,              //  Na-Ar
-    1.76, 1.54,                                                  //  K,Ca
-    1.33, 1.22, 1.21, 1.10, 1.07, 1.04, 1.00, 0.99, 1.01, 1.09,  //  Sc-Zn
-    1.12, 1.09, 1.15, 1.10, 1.14, 1.17,                          //  Ga-Kr
-    1.89, 1.67,                                                  //  Rb,Sr
-    1.47, 1.39, 1.32, 1.24, 1.15, 1.13, 1.13, 1.08, 1.15, 1.23,  //  Y-Cd
-    1.28, 1.26, 1.26, 1.23, 1.32, 1.31,                          //  In-Xe
-    2.09, 1.76,                                                  //  Cs,Ba
-    1.62, 1.47, 1.58, 1.57, 1.56, 1.55, 1.51,                    //  La-Eu
-    1.52, 1.51, 1.50, 1.49, 1.49, 1.48, 1.53,                    //  Gd-Yb
-    1.46, 1.37, 1.31, 1.23, 1.18, 1.16, 1.11, 1.12, 1.13, 1.32,  //  Lu-Hg
-    1.30, 1.30, 1.36, 1.31, 1.38, 1.42,                          //  Tl-Rn
-    2.01, 1.81,                                                  //  Fr,Ra
-    1.67, 1.58, 1.52, 1.53, 1.54, 1.55, 1.49,                    //  Ac-Am
-    1.49, 1.51, 1.51, 1.48, 1.50, 1.56, 1.58,                    //  Cm-No
-    1.45, 1.41, 1.34, 1.29, 1.27, 1.21, 1.16, 1.15, 1.09, 1.22,  //  Lr-Cn
-    1.22, 1.29, 1.46, 1.58, 1.48, 1.41                           //  Nh-Og
+  0.32, 0.46,                                                  //  H,He
+  1.20, 0.94, 0.77, 0.75, 0.71, 0.63, 0.64, 0.67,              //  Li-Ne
+  1.40, 1.25, 1.13, 1.04, 1.10, 1.02, 0.99, 0.96,              //  Na-Ar
+  1.76, 1.54,                                                  //  K,Ca
+  1.33, 1.22, 1.21, 1.10, 1.07, 1.04, 1.00, 0.99, 1.01, 1.09,  //  Sc-Zn
+  1.12, 1.09, 1.15, 1.10, 1.14, 1.17,                          //  Ga-Kr
+  1.89, 1.67,                                                  //  Rb,Sr
+  1.47, 1.39, 1.32, 1.24, 1.15, 1.13, 1.13, 1.08, 1.15, 1.23,  //  Y-Cd
+  1.28, 1.26, 1.26, 1.23, 1.32, 1.31,                          //  In-Xe
+  2.09, 1.76,                                                  //  Cs,Ba
+  1.62, 1.47, 1.58, 1.57, 1.56, 1.55, 1.51,                    //  La-Eu
+  1.52, 1.51, 1.50, 1.49, 1.49, 1.48, 1.53,                    //  Gd-Yb
+  1.46, 1.37, 1.31, 1.23, 1.18, 1.16, 1.11, 1.12, 1.13, 1.32,  //  Lu-Hg
+  1.30, 1.30, 1.36, 1.31, 1.38, 1.42,                          //  Tl-Rn
+  2.01, 1.81,                                                  //  Fr,Ra
+  1.67, 1.58, 1.52, 1.53, 1.54, 1.55, 1.49,                    //  Ac-Am
+  1.49, 1.51, 1.51, 1.48, 1.50, 1.56, 1.58,                    //  Cm-No
+  1.45, 1.41, 1.34, 1.29, 1.27, 1.21, 1.16, 1.15, 1.09, 1.22,  //  Lr-Cn
+  1.22, 1.29, 1.46, 1.58, 1.48, 1.41                           //  Nh-Og
 };
 
 /**
@@ -71,69 +71,69 @@ static const double covalent_rad_d3[119]{ 0.0,
  * rad = covalent_rad_d3 * 4.0/3.0 * aatoau
  */
 static const double rad[119]{ 0.0,
-    0.8062831465047213,     1.1590320231005369,     3.0235617993927044,
-    2.3684567428576182,     1.9401188212769855,     1.8897261246204402,
-    1.7889407313073502,     1.5873699446811698,     1.6125662930094427,
-    1.6881553379942602,     3.5274887659581555,     3.1495435410340673,
-    2.8471873610947966,     2.6204202261403440,     2.7715983161099791,
-    2.5700275294837986,     2.4944384844989811,     2.4188494395141635,
-    4.4345573057759662,     3.8802376425539711,     3.3511143276602473,
-    3.0739544960492493,     3.0487581477209771,     2.7715983161099791,
-    2.6960092711251620,     2.6204202261403440,     2.5196348328272538,
-    2.4944384844989811,     2.5448311811555264,     2.7464019677817069,
-    2.8219910127665244,     2.7464019677817069,     2.8975800577513415,
-    2.7715983161099791,     2.8723837094230693,     2.9479727544078864,
-    4.7621098340435095,     4.2077901708215135,     3.7038632042560633,
-    3.5022924176298824,     3.3259179793319751,     3.1243471927057946,
-    2.8975800577513415,     2.8471873610947966,     2.8471873610947966,
-    2.7212056194534342,     2.8975800577513415,     3.0991508443775224,
-    3.2251325860188853,     3.1747398893623395,     3.1747398893623395,
-    3.0991508443775224,     3.3259179793319751,     3.3007216310037024,
-    5.2660368006089602,     4.4345573057759662,     4.0818084291801515,
-    3.7038632042560633,     3.9810230358670613,     3.9558266875387886,
-    3.9306303392105164,     3.9054339908822433,     3.8046485975691535,
-    3.8298449458974257,     3.8046485975691535,     3.7794522492408804,
-    3.7542559009126082,     3.7542559009126082,     3.7290595525843355,
-    3.8550412942256984,     3.6786668559277906,     3.4518997209733380,
-    3.3007216310037024,     3.0991508443775224,     2.9731691027361595,
-    2.9227764060796142,     2.7967946644382522,     2.8219910127665244,
-    2.8471873610947966,     3.3259179793319751,     3.2755252826754302,
-    3.2755252826754302,     3.4267033726450653,     3.3007216310037024,
-    3.4770960693016097,     3.5778814626147004,     5.0644660139827797,
-    4.5605390474173291,     4.2077901708215135,     3.9810230358670613,
-    3.8298449458974257,     3.8550412942256984,     3.8802376425539711,
-    3.9054339908822433,     3.7542559009126082,     3.7542559009126082,
-    3.8046485975691535,     3.8046485975691535,     3.7290595525843355,
-    3.7794522492408804,     3.9306303392105164,     3.9810230358670613,
-    3.6534705075995175,     3.5526851142864277,     3.3763106759885204,
-    3.2503289343471575,     3.1999362376906122,     3.0487581477209771,
-    2.9227764060796142,     2.8975800577513415,     2.7464019677817069,
-    3.0739544960492493,     3.4267033726450653,     3.6030778109429726,
-    3.6786668559277906,     3.9810230358670613,     3.7290595525843355,
-    3.9558266875387886 
+  0.8062831465047213,     1.1590320231005369,     3.0235617993927044,
+  2.3684567428576182,     1.9401188212769855,     1.8897261246204402,
+  1.7889407313073502,     1.5873699446811698,     1.6125662930094427,
+  1.6881553379942602,     3.5274887659581555,     3.1495435410340673,
+  2.8471873610947966,     2.6204202261403440,     2.7715983161099791,
+  2.5700275294837986,     2.4944384844989811,     2.4188494395141635,
+  4.4345573057759662,     3.8802376425539711,     3.3511143276602473,
+  3.0739544960492493,     3.0487581477209771,     2.7715983161099791,
+  2.6960092711251620,     2.6204202261403440,     2.5196348328272538,
+  2.4944384844989811,     2.5448311811555264,     2.7464019677817069,
+  2.8219910127665244,     2.7464019677817069,     2.8975800577513415,
+  2.7715983161099791,     2.8723837094230693,     2.9479727544078864,
+  4.7621098340435095,     4.2077901708215135,     3.7038632042560633,
+  3.5022924176298824,     3.3259179793319751,     3.1243471927057946,
+  2.8975800577513415,     2.8471873610947966,     2.8471873610947966,
+  2.7212056194534342,     2.8975800577513415,     3.0991508443775224,
+  3.2251325860188853,     3.1747398893623395,     3.1747398893623395,
+  3.0991508443775224,     3.3259179793319751,     3.3007216310037024,
+  5.2660368006089602,     4.4345573057759662,     4.0818084291801515,
+  3.7038632042560633,     3.9810230358670613,     3.9558266875387886,
+  3.9306303392105164,     3.9054339908822433,     3.8046485975691535,
+  3.8298449458974257,     3.8046485975691535,     3.7794522492408804,
+  3.7542559009126082,     3.7542559009126082,     3.7290595525843355,
+  3.8550412942256984,     3.6786668559277906,     3.4518997209733380,
+  3.3007216310037024,     3.0991508443775224,     2.9731691027361595,
+  2.9227764060796142,     2.7967946644382522,     2.8219910127665244,
+  2.8471873610947966,     3.3259179793319751,     3.2755252826754302,
+  3.2755252826754302,     3.4267033726450653,     3.3007216310037024,
+  3.4770960693016097,     3.5778814626147004,     5.0644660139827797,
+  4.5605390474173291,     4.2077901708215135,     3.9810230358670613,
+  3.8298449458974257,     3.8550412942256984,     3.8802376425539711,
+  3.9054339908822433,     3.7542559009126082,     3.7542559009126082,
+  3.8046485975691535,     3.8046485975691535,     3.7290595525843355,
+  3.7794522492408804,     3.9306303392105164,     3.9810230358670613,
+  3.6534705075995175,     3.5526851142864277,     3.3763106759885204,
+  3.2503289343471575,     3.1999362376906122,     3.0487581477209771,
+  2.9227764060796142,     2.8975800577513415,     2.7464019677817069,
+  3.0739544960492493,     3.4267033726450653,     3.6030778109429726,
+  3.6786668559277906,     3.9810230358670613,     3.7290595525843355,
+  3.9558266875387886 
 };
 
 // pauling EN's
 static const double en[119]{
-    0.0, 2.20, 3.00,                                             //  H,He
-    0.98, 1.57, 2.04, 2.55, 3.04, 3.44, 3.98, 4.50,              //  Li-Ne
-    0.93, 1.31, 1.61, 1.90, 2.19, 2.58, 3.16, 3.50,              //  Na-Ar
-    0.82, 1.00,                                                  //  K,Ca
-    1.36, 1.54, 1.63, 1.66, 1.55, 1.83, 1.88, 1.91, 1.90, 1.65,  //  Sc-Zn
-    1.81, 2.01, 2.18, 2.55, 2.96, 3.00,                          //  Ga-Kr
-    0.82, 0.95,                                                  //  Rb,Sr
-    1.22, 1.33, 1.60, 2.16, 1.90, 2.20, 2.28, 2.20, 1.93, 1.69,  //  Y-Cd
-    1.78, 1.96, 2.05, 2.10, 2.66, 2.60,                          //  In-Xe
-    0.79, 0.89,                                                  //  Cs,Ba
-    1.10, 1.12, 1.13, 1.14, 1.15, 1.17, 1.18,                    //  La-Eu
-    1.20, 1.21, 1.22, 1.23, 1.24, 1.25, 1.26,                    //  Gd-Yb
-    1.27, 1.30, 1.50, 2.36, 1.90, 2.20, 2.20, 2.28, 2.54, 2.00,  //  Lu-Hg
-    1.62, 2.33, 2.02, 2.00, 2.20, 2.20,                          //  Tl-Rn
-    // only dummies below
-    1.50, 1.50,                                                  //  Fr,Ra
-    1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50,                    //  Ac-Am
-    1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50,                    //  Cm-No
-    1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50,  //  Rf-Cn
+  0.0, 2.20, 3.00,                                             //  H,He
+  0.98, 1.57, 2.04, 2.55, 3.04, 3.44, 3.98, 4.50,              //  Li-Ne
+  0.93, 1.31, 1.61, 1.90, 2.19, 2.58, 3.16, 3.50,              //  Na-Ar
+  0.82, 1.00,                                                  //  K,Ca
+  1.36, 1.54, 1.63, 1.66, 1.55, 1.83, 1.88, 1.91, 1.90, 1.65,  //  Sc-Zn
+  1.81, 2.01, 2.18, 2.55, 2.96, 3.00,                          //  Ga-Kr
+  0.82, 0.95,                                                  //  Rb,Sr
+  1.22, 1.33, 1.60, 2.16, 1.90, 2.20, 2.28, 2.20, 1.93, 1.69,  //  Y-Cd
+  1.78, 1.96, 2.05, 2.10, 2.66, 2.60,                          //  In-Xe
+  0.79, 0.89,                                                  //  Cs,Ba
+  1.10, 1.12, 1.13, 1.14, 1.15, 1.17, 1.18,                    //  La-Eu
+  1.20, 1.21, 1.22, 1.23, 1.24, 1.25, 1.26,                    //  Gd-Yb
+  1.27, 1.30, 1.50, 2.36, 1.90, 2.20, 2.20, 2.28, 2.54, 2.00,  //  Lu-Hg
+  1.62, 2.33, 2.02, 2.00, 2.20, 2.20,                          //  Tl-Rn
+  // only dummies below
+  1.50, 1.50,                                                  //  Fr,Ra
+  1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50,                    //  Ac-Am
+  1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50,                    //  Cm-No
+  1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50, 1.50,  //  Rf-Cn
     1.50, 1.50, 1.50, 1.50, 1.50, 1.50                           //  Nh-Og
 };
 
@@ -375,15 +375,15 @@ int cut_coordination_number(
   }
 
   return EXIT_SUCCESS;
-};
+}
 
 
 inline double log_cn_cut(const double cn_max, const double cn) {
   return log(1.0 + exp(cn_max)) - log(1.0 + exp(cn_max - cn));
-};
+}
 
 inline double dlog_cn_cut(const double cn_max, const double cn) {
   return exp(cn_max) / (exp(cn_max) + exp(cn));
-};
+}
 
-}  // namespace dftd
+}  // namespace dftd4
diff --git a/src/dftd_readxyz.cpp b/src/dftd_readxyz.cpp
index 37fe4e2..746bb7a 100644
--- a/src/dftd_readxyz.cpp
+++ b/src/dftd_readxyz.cpp
@@ -26,7 +26,7 @@
 #include "dftd_econv.h"
 #include "dftd_geometry.h"
 
-void read_xyzfile(const std::string& name, dftd::TMolecule& mol) {
+void read_xyzfile(const std::string& name, dftd4::TMolecule& mol) {
   std::ifstream geo;
   std::stringstream buffer;
   std::string line;
diff --git a/src/program_dftd.cpp b/src/program_dftd.cpp
index 0655b53..805943a 100644
--- a/src/program_dftd.cpp
+++ b/src/program_dftd.cpp
@@ -113,11 +113,11 @@ int main(int argc, char **argv) {
   std::string func;
   bool lverbose{false};
   bool lmbd{true}, lgrad{false};
-  dftd::dparam par;  // damping parameter for DFT-D4 calculation
-  dftd::TMolecule mol;
-  int info;
+  dftd4::dparam par;  // damping parameter for DFT-D4 calculation
+  dftd4::TMolecule mol;
+  int info{0};
+  int charge{0};
   double energy{0.0};
-  int charge;
 
   // check for complete command line
   // we need at least the program name and an input file
@@ -140,7 +140,7 @@ int main(int argc, char **argv) {
   }
   if (args.getflag("--func")) {
     func = args.getopt("--func");
-    dftd::d4par(func, par, lmbd);
+    dftd4::d4par(func, par, lmbd);
   }
   // last argument is assumed to filename since
   // dftd4 [options] <file>
@@ -150,9 +150,12 @@ int main(int argc, char **argv) {
   read_xyzfile(fname, mol);
 
   // initialize default cutoffs
-  dftd::TCutoff cutoff;
+  dftd4::TCutoff cutoff;
+
+  // molecular information
+  dftd4::TMolInfo dat = dftd4::TMolInfo(charge);
   
-  info = dftd::get_dispersion(mol, par, charge, cutoff, energy, nullptr);
+  info = dftd4::DFTVDW_D4(dat, mol, par, cutoff, energy, nullptr);
   if (info != 0) return EXIT_FAILURE;
 
   std::cout << "Dispersion energy: " << energy << " Eh\n";
diff --git a/test/test_disp.cpp b/test/test_disp.cpp
index e120e6f..4bd6acf 100644
--- a/test/test_disp.cpp
+++ b/test/test_disp.cpp
@@ -23,7 +23,7 @@
 #include "test_disp.h"
 #include "util.h"
 
-using namespace dftd;
+using namespace dftd4;
 
 
 int test_energy( 
@@ -33,15 +33,16 @@ int test_energy(
   const int charge,
   const double ref
 ) {
-  int info = 0;
-  bool lmbd = true;
-  double energy = 0.0;
+  int info{0};
+  bool latm{true};
+  double energy{0.0};
 
   // BP86 parameters
   dparam par;
-  d4par("bp86", par, lmbd);
+  d4par("bp86", par, latm);
 
   // assemble molecule
+  TMolInfo dat = TMolInfo(charge);
   TMolecule mol;
   info = get_molecule(n, atoms, coord, mol);
   if (!info == EXIT_SUCCESS) return info;
@@ -50,7 +51,7 @@ int test_energy(
   cutoff.disable(); // do not use cutoffs
 
   // dispersion main function
-  info = get_dispersion(mol, par, charge, cutoff, energy, nullptr);
+  info = DFTVDW_D4(dat, mol, par, cutoff, energy, nullptr);
   if (!info == EXIT_SUCCESS) return info;
 
   if (check(energy, ref) == EXIT_FAILURE) {
@@ -67,10 +68,10 @@ int test_disp() {
   info = test_energy(water_n, water_atoms, water_coord, water_charge, water_ref_energy);
   if (!info == EXIT_SUCCESS) return info;
   
-  info = test_energy(16, mb16_43_01_atoms, mb16_43_01_coord, mb16_43_01_charge, mb16_43_01_ref_energy);
+  info = test_energy(mb16_43_01_n, mb16_43_01_atoms, mb16_43_01_coord, mb16_43_01_charge, mb16_43_01_ref_energy);
   if (!info == EXIT_SUCCESS) return info;
 
-  info = test_energy(22, rost61_m1_atoms, rost61_m1_coord, rost61_m1_charge, rost61_m1_ref_energy);
+  info = test_energy(rost61_m1_n, rost61_m1_atoms, rost61_m1_coord, rost61_m1_charge, rost61_m1_ref_energy);
   if (!info == EXIT_SUCCESS) return info;
 
   return EXIT_SUCCESS;
diff --git a/test/test_grad.cpp b/test/test_grad.cpp
index 55b3ef3..8090f8f 100644
--- a/test/test_grad.cpp
+++ b/test/test_grad.cpp
@@ -24,9 +24,9 @@
 #include "test_grad.h"
 #include "util.h"
 
-using namespace dftd;
+using namespace dftd4;
 
-int test_numgrad(TMolecule &mol, const int charge, const dparam &par) {
+int test_numgrad(const TMolInfo &dat, TMolecule &mol, const dparam &par) {
   int info{0};
   double energy{0.0}, er{0.0}, el{0.0};
   double step{1.0e-6};
@@ -44,10 +44,10 @@ int test_numgrad(TMolecule &mol, const int charge, const dparam &par) {
       el = 0.0;
 
       mol.xyz(i, c) += step;
-      get_dispersion(mol, par, charge, cutoff, er, nullptr);
+      DFTVDW_D4(dat, mol, par, cutoff, er, nullptr);
 
       mol.xyz(i, c) = mol.xyz(i, c) - 2*step;
-      get_dispersion(mol, par, charge, cutoff, el, nullptr);
+      DFTVDW_D4(dat, mol, par, cutoff, el, nullptr);
 
       mol.xyz(i, c) = mol.xyz(i, c) + step;
       numgrad(i, c) = 0.5 * (er - el) / step;
@@ -57,7 +57,7 @@ int test_numgrad(TMolecule &mol, const int charge, const dparam &par) {
   // analytical gradient
   double* d4grad = new double[3*mol.NAtoms];
   for (int i = 0; i < 3*mol.NAtoms; i++) d4grad[i] = 0.0;
-  info = get_dispersion(mol, par, charge, cutoff, energy, d4grad);
+  info = DFTVDW_D4(dat, mol, par, cutoff, energy, d4grad);
   if (!info == EXIT_SUCCESS) return info;
 
   // check translational invariance of analytical gradient
@@ -104,9 +104,6 @@ int is_trans_invar(const TMolecule& mol, double gradient[]) {
 
 
 int test_pbed4_mb01() {
-  int charge{0};
-  int info{0};
-
   // PBE-D4(EEQ) parameters
   dparam par;
   par.s6 = 1.0;
@@ -117,46 +114,40 @@ int test_pbed4_mb01() {
   par.a2 = 4.80688534;
 
   // assemble molecule
+  TMolInfo dat = TMolInfo(mb16_43_01_charge);
   TMolecule mol;
-  info = get_molecule(mb16_43_01_n, mb16_43_01_atoms, mb16_43_01_coord, mol);
+  int info = get_molecule(mb16_43_01_n, mb16_43_01_atoms, mb16_43_01_coord, mol);
   if (!info == EXIT_SUCCESS) return info;
-  charge = mb16_43_01_charge; 
 
-  return test_numgrad(mol, charge, par);
+  return test_numgrad(dat, mol, par);
 };
 
 int test_bp86d4atm_water() {
-  int charge{0};
-  int info{0};
-
   // PBE-D4(EEQ)-ATM parameters
   dparam par;
   d4par("bp86", par, true);
 
   // assemble molecule
+  TMolInfo dat = TMolInfo(water_charge);
   TMolecule mol;
-  info = get_molecule(water_n, water_atoms, water_coord, mol);
+  int info = get_molecule(water_n, water_atoms, water_coord, mol);
   if (!info == EXIT_SUCCESS) return info;
-  charge = water_charge; 
 
-  return test_numgrad(mol, charge, par);
+  return test_numgrad(dat, mol, par);
 }
 
 int test_tpss0d4mbd_rost61m1() {
-  int charge{0};
-  int info{0};
-
   // TPSS0-D4(EEQ)-MBD parameters
   dparam par;
   d4par("tpss0", par, false);
 
   // assemble molecule
+  TMolInfo dat = TMolInfo(rost61_m1_charge);
   TMolecule mol;
-  info = get_molecule(rost61_m1_n, rost61_m1_atoms, rost61_m1_coord, mol);
+  int info = get_molecule(rost61_m1_n, rost61_m1_atoms, rost61_m1_coord, mol);
   if (!info == EXIT_SUCCESS) return info;
-  charge = rost61_m1_charge; 
 
-  return test_numgrad(mol, charge, par);
+  return test_numgrad(dat, mol, par);
 }
 
 int test_grad() {
diff --git a/test/test_grad.h b/test/test_grad.h
index bbdc303..9d0dc13 100644
--- a/test/test_grad.h
+++ b/test/test_grad.h
@@ -21,9 +21,9 @@
 #include "dftd_dispersion.h"
 #include "molecules.h"
 
-using namespace dftd;
+using namespace dftd4;
 
-extern int test_numgrad(TMolecule &mol, const int charge, const dparam &par);
+extern int test_numgrad(const TMolInfo &dat, TMolecule &mol, const dparam &par);
 extern int is_trans_invar(const TMolecule &mol, double gradient[]);
 
 extern int test_bp86d4atm_water(void);
diff --git a/test/test_ncoord.cpp b/test/test_ncoord.cpp
index 58e9ceb..f6b0e78 100644
--- a/test/test_ncoord.cpp
+++ b/test/test_ncoord.cpp
@@ -25,7 +25,7 @@
 #include "test_ncoord.h"
 #include "util.h"
 
-using namespace dftd;
+using namespace dftd4;
 
 
 int test_cn(
diff --git a/test/test_param.cpp b/test/test_param.cpp
index 5089eb9..8c9533d 100644
--- a/test/test_param.cpp
+++ b/test/test_param.cpp
@@ -27,7 +27,7 @@
 #include "test_param.h"
 #include "util.h"
 
-using namespace dftd;
+using namespace dftd4;
 
 int test_rational_damping(const double ref[], TCutoff cutoff) {
   int info;
@@ -36,18 +36,17 @@ int test_rational_damping(const double ref[], TCutoff cutoff) {
   dparam par;
 
   // assemble molecule
+  TMolInfo dat = TMolInfo(upu23_0a_charge);
   TMolecule mol;
   info = get_molecule(upu23_0a_n, upu23_0a_atoms, upu23_0a_coord, mol);
   if (!info == EXIT_SUCCESS) return info;
-  charge = upu23_0a_charge;
-
 
   for (int i = 0; i < nfuncs; i++) {
     std::string func = funcs[i];
     d4par(func, par, true);
 
     energy = 0.0;
-    info = get_dispersion(mol, par, charge, cutoff, energy, nullptr);
+    info = DFTVDW_D4(dat, mol, par, cutoff, energy, nullptr);
     if (!info == EXIT_SUCCESS) return info;
 
     if (check(energy, ref[i]) == EXIT_FAILURE) {
diff --git a/test/test_param.h b/test/test_param.h
index 6de5cc0..37a333d 100644
--- a/test/test_param.h
+++ b/test/test_param.h
@@ -22,7 +22,7 @@
 #include <dftd_geometry.h>
 #include <dftd_dispersion.h>
 
-using namespace dftd;
+using namespace dftd4;
 
 extern int test_param(void);
 extern int test_rational_damping(const double ref[], TCutoff cutoff);
diff --git a/test/util.cpp b/test/util.cpp
index e152d26..88df121 100644
--- a/test/util.cpp
+++ b/test/util.cpp
@@ -25,7 +25,7 @@
 
 #include "util.h"
 
-using namespace dftd;
+using namespace dftd4;
 
 
 int get_molecule(int n, const char atoms[][4], const double coord[], TMolecule& mol) {
diff --git a/test/util.h b/test/util.h
index f36f439..8d476c7 100644
--- a/test/util.h
+++ b/test/util.h
@@ -23,7 +23,7 @@
 #include <dftd_geometry.h>
 
 
-extern int get_molecule(int n, const char atoms[][4], const double coord[], dftd::TMolecule& mol);
+extern int get_molecule(int n, const char atoms[][4], const double coord[], dftd4::TMolecule& mol);
 
 extern bool check(double actual, double expected, double epsilon = 1e-12, bool rel = false);
 extern bool check(float actual, float expected, float epsilon = 1e-6, bool rel = false);