diff --git a/include/dftd_cblas.h b/include/dftd_cblas.h index 95f8b26..a9ed3a3 100644 --- a/include/dftd_cblas.h +++ b/include/dftd_cblas.h @@ -100,7 +100,8 @@ inline int BLAS_Add_Mat_x_Mat( const double alpha ) { // check for size 0 matrices - if (A.cols == 0 || A.rows == 0 || B.cols == 0 || B.rows == 0 || C.cols == 0 || C.rows == 0) + if (A.cols == 0 || A.rows == 0 || B.cols == 0 || B.rows == 0 || C.cols == 0 || + C.rows == 0) exit(EXIT_FAILURE); // check for transpositions @@ -150,7 +151,7 @@ inline int BLAS_Add_Mat_x_Mat( C.cols ); }; // B transposed - } // A not transposed + } // A not transposed else { if (!TransposeB) { // check dimensions for C=AT*B diff --git a/include/dftd_damping.h b/include/dftd_damping.h index 70be056..f425fca 100644 --- a/include/dftd_damping.h +++ b/include/dftd_damping.h @@ -36,7 +36,6 @@ namespace dftd4 { * @param latm Switch for D4-ATM (true) or D4-MBD (false) parameters. * @returns Exit status. */ -extern int - d4par(const std::string func, dftd4::dparam &par, bool latm = true); +extern int d4par(const std::string func, dftd4::dparam &par, bool latm = true); } // namespace dftd4 diff --git a/include/dftd_matrix.h b/include/dftd_matrix.h index a47db62..210825b 100644 --- a/include/dftd_matrix.h +++ b/include/dftd_matrix.h @@ -26,189 +26,189 @@ namespace dftd4 { // Define a vector template class TVector { - public: - int N; // Dimension of the vector - int ElementSize; // Size of each element in the vector - T *p; // The pointer to the vector - - TVector() { - N = 0; - p = nullptr; - ElementSize = sizeof(T); - } - ~TVector() { - if (p != nullptr) Delete(); - } - void NewVector(int VectorLength) { - if (VectorLength < 0) { std::exit(EXIT_FAILURE); } - if (p != nullptr && N == VectorLength) { - Init(); - } else { - Delete(); - if (VectorLength == 0) { return; } - // get new memory - p = new T[VectorLength]; - if (!p) { std::exit(EXIT_FAILURE); } - N = VectorLength; - Init(); + public: + int N; // Dimension of the vector + int ElementSize; // Size of each element in the vector + T *p; // The pointer to the vector + + TVector() { + N = 0; + p = nullptr; + ElementSize = sizeof(T); } - } - - // alias for NewVector - void New(int VectorLength) { return NewVector(VectorLength); } - void NewVec(int VectorLength) { return NewVector(VectorLength); } - - void Delete() { - if (p != nullptr && N != 0) { delete[] p; } - p = nullptr; - N = 0; - } - void CopyVec(const TVector &v) { - long int mem; - if (N != v.N) { - Delete(); - New(v.N); + ~TVector() { + if (p != nullptr) Delete(); } - if (v.N == 0) return; - mem = (long int)N * ElementSize; - std::memcpy(p, v.p, mem); - } - void Init() { - if (p != nullptr) { - long int mem = (long int)N * ElementSize; - std::memset(p, 0, mem); + void NewVector(int VectorLength) { + if (VectorLength < 0) { std::exit(EXIT_FAILURE); } + if (p != nullptr && N == VectorLength) { + Init(); + } else { + Delete(); + if (VectorLength == 0) { return; } + // get new memory + p = new T[VectorLength]; + if (!p) { std::exit(EXIT_FAILURE); } + N = VectorLength; + Init(); + } } - } - void Print(char name[]) { - printf("Vector printed: %s (%d)\n", name, N); - for (int i = 0; i < N; i++) { - printf("%+23.15e\n", p[i]); + // alias for NewVector + void New(int VectorLength) { return NewVector(VectorLength); } + void NewVec(int VectorLength) { return NewVector(VectorLength); } + + void Delete() { + if (p != nullptr && N != 0) { delete[] p; } + p = nullptr; + N = 0; + } + void CopyVec(const TVector &v) { + long int mem; + if (N != v.N) { + Delete(); + New(v.N); + } + if (v.N == 0) return; + mem = (long int)N * ElementSize; + std::memcpy(p, v.p, mem); + } + void Init() { + if (p != nullptr) { + long int mem = (long int)N * ElementSize; + std::memset(p, 0, mem); + } } - printf("\n"); - } - inline T &operator()(int i) { return p[i]; } - inline const T &operator()(int i) const { return p[i]; } - inline T &operator[](int i) { return p[i]; } - inline const T &operator[](int i) const { return p[i]; } + void Print(char name[]) { + printf("Vector printed: %s (%d)\n", name, N); + for (int i = 0; i < N; i++) { + printf("%+23.15e\n", p[i]); + } + printf("\n"); + } + + inline T &operator()(int i) { return p[i]; } + inline const T &operator()(int i) const { return p[i]; } + inline T &operator[](int i) { return p[i]; } + inline const T &operator[](int i) const { return p[i]; } }; // Define a normal matrix template class TMatrix { - public: - int rows, cols; // dimensions - int ElementSize; // Size of elements in matrix - T *p; // pointer to dynamic memory - - TMatrix() { - cols = 0; - rows = 0; - p = nullptr; - ElementSize = sizeof(T); - } - ~TMatrix() { - if (p != nullptr) Delete(); - } - - void NewMatrix(int r, int c) { - if (r < 0 || c < 0) std::exit(EXIT_FAILURE); - if (p != nullptr && r == rows && c == cols) { - Init(); - } else { - long int mem = (long int)r * (long int)c; - if (p != nullptr) Delete(); // Eventually delete old matrix - - if (mem == 0) return; // don't touch pointer if no memory is allocated - - p = new T[mem]; - if (!p) std::exit(EXIT_FAILURE); - rows = r; - cols = c; - Init(); + public: + int rows, cols; // dimensions + int ElementSize; // Size of elements in matrix + T *p; // pointer to dynamic memory + + TMatrix() { + cols = 0; + rows = 0; + p = nullptr; + ElementSize = sizeof(T); + } + ~TMatrix() { + if (p != nullptr) Delete(); + } + + void NewMatrix(int r, int c) { + if (r < 0 || c < 0) std::exit(EXIT_FAILURE); + if (p != nullptr && r == rows && c == cols) { + Init(); + } else { + long int mem = (long int)r * (long int)c; + if (p != nullptr) Delete(); // Eventually delete old matrix + + if (mem == 0) return; // don't touch pointer if no memory is allocated + + p = new T[mem]; + if (!p) std::exit(EXIT_FAILURE); + rows = r; + cols = c; + Init(); + } + return; } - return; - } - - 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() { - if (p != nullptr && rows * cols != 0) { delete[] p; } - rows = 0; - cols = 0; - p = nullptr; - } - - void Init() { - long int mem; - if (p != nullptr) { - mem = (long int)cols * (long int)rows * ElementSize; - std::memset(p, 0, mem); + + 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() { + if (p != nullptr && rows * cols != 0) { delete[] p; } + rows = 0; + cols = 0; + p = nullptr; } - } - - void Transpose() { - T x; - int i, j; - - if (p != nullptr) { - if (rows == cols) { - for (i = 0; i < rows; i++) { - for (j = 0; j < i; j++) { - x = p[i * cols + j]; - p[i * cols + j] = p[j * cols + i]; - p[j * cols + i] = x; - } // j - } // i - } // if NxN - else { - // for non-square matrix, we need an additional copy - TMatrix temp; - temp.CopyMat(*this); - NewMatrix(cols, rows); - for (i = 0; i < rows; i++) { - for (j = 0; j < cols; j++) { - p[i * cols + j] = temp.p[j * cols + i]; - } // j - } // i + + void Init() { + long int mem; + if (p != nullptr) { + mem = (long int)cols * (long int)rows * ElementSize; + std::memset(p, 0, mem); } - } // if data is loaded - } // for NxN matrices transpose elements + } + + void Transpose() { + T x; + int i, j; + + if (p != nullptr) { + if (rows == cols) { + for (i = 0; i < rows; i++) { + for (j = 0; j < i; j++) { + x = p[i * cols + j]; + p[i * cols + j] = p[j * cols + i]; + p[j * cols + i] = x; + } // j + } // i + } // if NxN + else { + // for non-square matrix, we need an additional copy + TMatrix temp; + temp.CopyMat(*this); + NewMatrix(cols, rows); + for (i = 0; i < rows; i++) { + for (j = 0; j < cols; j++) { + p[i * cols + j] = temp.p[j * cols + i]; + } // j + } // i + } + } // if data is loaded + } // for NxN matrices transpose elements - void CopyMat(const TMatrix &m) { - long int mem; + void CopyMat(const TMatrix &m) { + long int mem; - if ((m.rows != rows) || (m.cols != cols)) { - Delete(); - New(m.rows, m.cols); + if ((m.rows != rows) || (m.cols != cols)) { + Delete(); + New(m.rows, m.cols); + } + mem = (long int)rows * (long int)cols * ElementSize; + if (mem == 0) return; + std::memcpy(p, m.p, mem); } - mem = (long int)rows * (long int)cols * ElementSize; - if (mem == 0) return; - std::memcpy(p, m.p, mem); - } - - void Print(const char name[] = "unknown") { - printf("Matrix printed: %s (%d, %d)\n", name, rows, cols); - for (int i = 0; i < rows; i++) { - for (int j = 0; j < cols; j++) { - printf("%+23.15e", p[i * cols + j]); - if (j == cols - 1) { - printf("\n"); - } else { - printf(" "); + + void Print(const char name[] = "unknown") { + printf("Matrix printed: %s (%d, %d)\n", name, rows, cols); + for (int i = 0; i < rows; i++) { + for (int j = 0; j < cols; j++) { + printf("%+23.15e", p[i * cols + j]); + if (j == cols - 1) { + printf("\n"); + } else { + printf(" "); + } } } + printf("\n"); } - printf("\n"); - } - inline T &operator()(int i, int j) { return p[i * cols + j]; } - inline const T &operator()(int i, int j) const { return p[i * cols + j]; } - inline T *operator[](int i) { return p + i * cols; } + inline T &operator()(int i, int j) { return p[i * cols + j]; } + inline const T &operator()(int i, int j) const { return p[i * cols + j]; } + inline T *operator[](int i) { return p + i * cols; } }; } // namespace dftd4 diff --git a/include/dftd_model.h b/include/dftd_model.h index 4e109c4..12c9b24 100644 --- a/include/dftd_model.h +++ b/include/dftd_model.h @@ -76,13 +76,10 @@ class TD4Model { extern inline double trapzd(const double a[23], const double b[23]); -extern inline double - weight_cn(double wf, double cn, double cnref); +extern inline double weight_cn(double wf, double cn, double cnref); -extern inline double - zeta(double a, double c, double qref, double qmod); -extern inline double - dzeta(double a, double c, double qref, double qmod); +extern inline double zeta(double a, double c, double qref, double qmod); +extern inline double dzeta(double a, double c, double qref, double qmod); extern int get_max_ref(const TMolecule &mol, int &mref); diff --git a/include/dftd_parameters.h b/include/dftd_parameters.h index a89e5cd..112dc2b 100644 --- a/include/dftd_parameters.h +++ b/include/dftd_parameters.h @@ -48,7 +48,6 @@ static const double zeff[119]{ 9,10,11,30,31,32,33,34,35,36,37,38,39,40,41,42,43, // Fr-Lr 12,13,14,15,16,17,18,19,20,21,22,23,24,25,26 // Rf-Og }; -// clang-format on /** * PBE0/def2-QZVP atomic values calculated by S. Grimme in Gaussian (2010) @@ -63,10 +62,10 @@ static const double zeff[119]{ * These values are actually never used in the code. Only r4r2 is used. */ static const double r2r4[119]{ - 0.0, // dummy - 8.0589, 3.4698, // H,He + 0.0, // dummy + 8.0589, 3.4698, // H,He 29.0974, 14.8517, 11.8799, 7.8715, 5.5588, 4.7566, 3.8025, - 3.1036, // Li-Ne + 3.1036, // Li-Ne 26.1552, 17.2304, 17.7210, 12.7442, 9.5361, 8.1652, 6.7463, 5.6004, // Na-Ar 29.2012, 22.3934, // K,Ca @@ -90,6 +89,7 @@ static const double r2r4[119]{ 7.5667, 6.9456, 6.3946, 5.9159, 5.4929, // -Cn 6.7286, 6.5144, 10.9169, 10.3600, 9.4723, 8.6641 // Nh-Og }; +// clang-format on // r4/r2 expectation value static const double r4r2[119] = { @@ -6149,18 +6149,24 @@ static const double *refcovcn[MAXELEMENT]{ refcovcn102, refcovcn103, }; static const double *refq_eeq[MAXELEMENT]{ - nullptr, refq_eeq1, refq_eeq2, refq_eeq3, refq_eeq4, refq_eeq5, refq_eeq6, refq_eeq7, refq_eeq8, - refq_eeq9, refq_eeq10, refq_eeq11, refq_eeq12, refq_eeq13, refq_eeq14, refq_eeq15, refq_eeq16, refq_eeq17, - refq_eeq18, refq_eeq19, refq_eeq20, refq_eeq21, refq_eeq22, refq_eeq23, refq_eeq24, refq_eeq25, refq_eeq26, - refq_eeq27, refq_eeq28, refq_eeq29, refq_eeq30, refq_eeq31, refq_eeq32, refq_eeq33, refq_eeq34, refq_eeq35, - refq_eeq36, refq_eeq37, refq_eeq38, refq_eeq39, refq_eeq40, refq_eeq41, refq_eeq42, refq_eeq43, refq_eeq44, - refq_eeq45, refq_eeq46, refq_eeq47, refq_eeq48, refq_eeq49, refq_eeq50, refq_eeq51, refq_eeq52, refq_eeq53, - refq_eeq54, refq_eeq55, refq_eeq56, refq_eeq57, refq_eeq58, refq_eeq59, refq_eeq60, refq_eeq61, refq_eeq62, - refq_eeq63, refq_eeq64, refq_eeq65, refq_eeq66, refq_eeq67, refq_eeq68, refq_eeq69, refq_eeq70, refq_eeq71, - refq_eeq72, refq_eeq73, refq_eeq74, refq_eeq75, refq_eeq76, refq_eeq77, refq_eeq78, refq_eeq79, refq_eeq80, - refq_eeq81, refq_eeq82, refq_eeq83, 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, + nullptr, refq_eeq1, refq_eeq2, refq_eeq3, refq_eeq4, refq_eeq5, + refq_eeq6, refq_eeq7, refq_eeq8, refq_eeq9, refq_eeq10, refq_eeq11, + refq_eeq12, refq_eeq13, refq_eeq14, refq_eeq15, refq_eeq16, refq_eeq17, + refq_eeq18, refq_eeq19, refq_eeq20, refq_eeq21, refq_eeq22, refq_eeq23, + refq_eeq24, refq_eeq25, refq_eeq26, refq_eeq27, refq_eeq28, refq_eeq29, + refq_eeq30, refq_eeq31, refq_eeq32, refq_eeq33, refq_eeq34, refq_eeq35, + refq_eeq36, refq_eeq37, refq_eeq38, refq_eeq39, refq_eeq40, refq_eeq41, + refq_eeq42, refq_eeq43, refq_eeq44, refq_eeq45, refq_eeq46, refq_eeq47, + refq_eeq48, refq_eeq49, refq_eeq50, refq_eeq51, refq_eeq52, refq_eeq53, + refq_eeq54, refq_eeq55, refq_eeq56, refq_eeq57, refq_eeq58, refq_eeq59, + refq_eeq60, refq_eeq61, refq_eeq62, refq_eeq63, refq_eeq64, refq_eeq65, + refq_eeq66, refq_eeq67, refq_eeq68, refq_eeq69, refq_eeq70, refq_eeq71, + refq_eeq72, refq_eeq73, refq_eeq74, refq_eeq75, refq_eeq76, refq_eeq77, + refq_eeq78, refq_eeq79, refq_eeq80, refq_eeq81, refq_eeq82, refq_eeq83, + 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, }; static const double *refsq[MAXELEMENT]{ nullptr, refsq1, refsq2, refsq3, refsq4, refsq5, refsq6, refsq7, diff --git a/src/damping/dftd_atm.cpp b/src/damping/dftd_atm.cpp index cc1e2e6..d0faa1e 100644 --- a/src/damping/dftd_atm.cpp +++ b/src/damping/dftd_atm.cpp @@ -20,63 +20,73 @@ * @brief Three-body (ATM) dispersion */ - #include +#include "damping/dftd_atm.h" #include "dftd_cblas.h" -#include "dftd_eeq.h" #include "dftd_dispersion.h" +#include "dftd_eeq.h" #include "dftd_geometry.h" #include "dftd_matrix.h" #include "dftd_ncoord.h" #include "dftd_parameters.h" -#include "damping/dftd_atm.h" namespace dftd4 { int get_atm_dispersion( - const TMolecule& mol, - const TMatrix& dist, + const TMolecule &mol, + const TMatrix &dist, const double cutoff, const double s9, const double a1, const double a2, const double alp, - const TMatrix& c6, - const TMatrix& dc6dcn, - const TMatrix& dc6dq, - TVector& energy, - TVector& dEdcn, - TVector& dEdq, - TVector& gradient, + const TMatrix &c6, + const TMatrix &dc6dcn, + const TMatrix &dc6dq, + TVector &energy, + TVector &dEdcn, + TVector &dEdq, + TVector &gradient, bool lgrad /*= false*/ ) { int info{0}; if (lgrad) { info = get_atm_dispersion_derivs( - mol, dist, cutoff, s9, a1, a2, alp, c6, dc6dcn, dc6dq, energy, - dEdcn, dEdq, gradient + mol, + dist, + cutoff, + s9, + a1, + a2, + alp, + c6, + dc6dcn, + dc6dq, + energy, + dEdcn, + dEdq, + gradient ); } else { - info = get_atm_dispersion_energy( - mol, dist, cutoff, s9, a1, a2, alp, c6, energy - ); + info = + get_atm_dispersion_energy(mol, dist, cutoff, s9, a1, a2, alp, c6, energy); } return info; } int get_atm_dispersion_energy( - const TMolecule& mol, - const TMatrix& dist, + const TMolecule &mol, + const TMatrix &dist, const double cutoff, const double s9, const double a1, const double a2, const double alp, - const TMatrix& c6, - TVector& energy + const TMatrix &c6, + TVector &energy ) { int izp{0}, jzp{0}, kzp{0}; double r0ij{0.0}, r0ik{0.0}, r0jk{0.0}, r0ijk{0.0}; @@ -91,7 +101,7 @@ int get_atm_dispersion_energy( izp = mol.at(iat); for (int jat = 0; jat != iat; jat++) { rij = dist(iat, jat); - if (rij > cutoff) continue; + if (rij > cutoff) continue; r2ij = pow(rij, 2); jzp = mol.at(jat); @@ -100,7 +110,7 @@ int get_atm_dispersion_energy( for (int kat = 0; kat != jat; kat++) { rik = dist(iat, kat); - if (rik > cutoff) continue; + if (rik > cutoff) continue; rjk = dist(jat, kat); if (rjk > cutoff) continue; @@ -124,7 +134,9 @@ int get_atm_dispersion_energy( fdmp = 1.0 / (1.0 + 6.0 * pow(r0ijk / rijk, alp / 3.0)); ang = ((0.375 * (r2ij + r2jk - r2ik) * (r2ij + r2ik - r2jk) * - (r2ik + r2jk - r2ij) / r2ijk) + 1.0) / r3ijk; + (r2ik + r2jk - r2ij) / r2ijk) + + 1.0) / + r3ijk; e = ang * fdmp * c9ijk / 3.0 * triple; energy(iat) += e; @@ -137,22 +149,21 @@ int get_atm_dispersion_energy( return EXIT_SUCCESS; } - int get_atm_dispersion_derivs( - const TMolecule& mol, - const TMatrix& dist, + const TMolecule &mol, + const TMatrix &dist, const double cutoff, const double s9, const double a1, const double a2, const double alp, - const TMatrix& c6, - const TMatrix& dc6dcn, - const TMatrix& dc6dq, - TVector& energy, - TVector& dEdcn, - TVector& dEdq, - TVector& gradient + const TMatrix &c6, + const TMatrix &dc6dcn, + const TMatrix &dc6dq, + TVector &energy, + TVector &dEdcn, + TVector &dEdq, + TVector &gradient ) { int izp{0}, jzp{0}, kzp{0}; double r0ij{0.0}, r0ik{0.0}, r0jk{0.0}, r0ijk{0.0}; @@ -174,7 +185,7 @@ int get_atm_dispersion_derivs( izp = mol.at(iat); for (int jat = 0; jat != iat; jat++) { rij = dist(iat, jat); - if (rij > cutoff) continue; + if (rij > cutoff) continue; r2ij = pow(rij, 2); jzp = mol.at(jat); @@ -183,7 +194,7 @@ int get_atm_dispersion_derivs( for (int kat = 0; kat != jat; kat++) { rik = dist(iat, kat); - if (rik > cutoff) continue; + if (rik > cutoff) continue; rjk = dist(jat, kat); if (rjk > cutoff) continue; @@ -209,7 +220,9 @@ int get_atm_dispersion_derivs( tmp = pow(r0ijk / rijk, alp / 3.0); fdmp = 1.0 / (1.0 + 6.0 * tmp); ang = ((0.375 * (r2ij + r2jk - r2ik) * (r2ij + r2ik - r2jk) * - (r2ik + r2jk - r2ij) / r2ijk) + 1.0) / r3ijk; + (r2ik + r2jk - r2ij) / r2ijk) + + 1.0) / + r3ijk; e = ang * fdmp * c9ijk * triple; energy(iat) -= e / 3.0; @@ -232,55 +245,61 @@ int get_atm_dispersion_derivs( dfdmp = -2.0 * alp * tmp * pow(fdmp, 2); // d/drij - dang = -0.375 * (pow(r2ij, 3) + pow(r2ij, 2) * (r2jk + r2ik) - + r2ij * (3.0 * pow(r2jk, 2) + 2.0 * r2jk*r2ik - + 3.0 * pow(r2ik, 2)) - - 5.0 * pow((r2jk - r2ik), 2) * (r2jk + r2ik)) / r5ijk; - dgxij = c9ijk * (-dang*fdmp + ang*dfdmp) / r2ij * xij; - dgyij = c9ijk * (-dang*fdmp + ang*dfdmp) / r2ij * yij; - dgzij = c9ijk * (-dang*fdmp + ang*dfdmp) / r2ij * zij; + dang = -0.375 * + (pow(r2ij, 3) + pow(r2ij, 2) * (r2jk + r2ik) + + r2ij * (3.0 * pow(r2jk, 2) + 2.0 * r2jk * r2ik + + 3.0 * pow(r2ik, 2)) - + 5.0 * pow((r2jk - r2ik), 2) * (r2jk + r2ik)) / + r5ijk; + dgxij = c9ijk * (-dang * fdmp + ang * dfdmp) / r2ij * xij; + dgyij = c9ijk * (-dang * fdmp + ang * dfdmp) / r2ij * yij; + dgzij = c9ijk * (-dang * fdmp + ang * dfdmp) / r2ij * zij; // d/drik - dang = -0.375 * (pow(r2ik, 3) + pow(r2ik, 2) * (r2jk + r2ij) - + r2ik * (3.0 * pow(r2jk, 2) + 2.0 * r2jk * r2ij - + 3.0 * pow(r2ij, 2)) - - 5.0 * pow((r2jk - r2ij), 2) * (r2jk + r2ij)) / r5ijk; + dang = -0.375 * + (pow(r2ik, 3) + pow(r2ik, 2) * (r2jk + r2ij) + + r2ik * (3.0 * pow(r2jk, 2) + 2.0 * r2jk * r2ij + + 3.0 * pow(r2ij, 2)) - + 5.0 * pow((r2jk - r2ij), 2) * (r2jk + r2ij)) / + r5ijk; dgxik = c9ijk * (-dang * fdmp + ang * dfdmp) / r2ik * xik; dgyik = c9ijk * (-dang * fdmp + ang * dfdmp) / r2ik * yik; dgzik = c9ijk * (-dang * fdmp + ang * dfdmp) / r2ik * zik; // d/drjk - dang = -0.375 * (pow(r2jk, 3) + pow(r2jk, 2)*(r2ik + r2ij) - + r2jk * (3.0 * pow(r2ik, 2) + 2.0 * r2ik * r2ij - + 3.0 * pow(r2ij, 2)) - - 5.0 * pow((r2ik - r2ij), 2) * (r2ik + r2ij)) / r5ijk; + dang = -0.375 * + (pow(r2jk, 3) + pow(r2jk, 2) * (r2ik + r2ij) + + r2jk * (3.0 * pow(r2ik, 2) + 2.0 * r2ik * r2ij + + 3.0 * pow(r2ij, 2)) - + 5.0 * pow((r2ik - r2ij), 2) * (r2ik + r2ij)) / + r5ijk; dgxjk = c9ijk * (-dang * fdmp + ang * dfdmp) / r2jk * xjk; dgyjk = c9ijk * (-dang * fdmp + ang * dfdmp) / r2jk * yjk; dgzjk = c9ijk * (-dang * fdmp + ang * dfdmp) / r2jk * zjk; - - gradient(3*iat + 0) += - dgxij - dgxik; - gradient(3*iat + 1) += - dgyij - dgyik; - gradient(3*iat + 2) += - dgzij - dgzik; - gradient(3*jat + 0) += dgxij - dgxjk; - gradient(3*jat + 1) += dgyij - dgyjk; - gradient(3*jat + 2) += dgzij - dgzjk; - gradient(3*kat + 0) += dgxik + dgxjk; - gradient(3*kat + 1) += dgyik + dgyjk; - gradient(3*kat + 2) += dgzik + dgzjk; - - dEdcn(iat) -= e * 0.5 - * (dc6dcn(iat, jat) / c6ij + dc6dcn(iat, kat) / c6ik); - dEdcn(jat) -= e * 0.5 - * (dc6dcn(jat, iat) / c6ij + dc6dcn(jat, kat) / c6jk); - dEdcn(kat) -= e * 0.5 - * (dc6dcn(kat, iat) / c6ik + dc6dcn(kat, jat) / c6jk); - - dEdq(iat) -= e * 0.5 - * (dc6dq(iat, jat) / c6ij + dc6dq(iat, kat) / c6ik); - dEdq(jat) -= e * 0.5 - * (dc6dq(jat, iat) / c6ij + dc6dq(jat, kat) / c6jk); - dEdq(kat) -= e * 0.5 - * (dc6dq(kat, iat) / c6ik + dc6dq(kat, jat) / c6jk); + + gradient(3 * iat + 0) += -dgxij - dgxik; + gradient(3 * iat + 1) += -dgyij - dgyik; + gradient(3 * iat + 2) += -dgzij - dgzik; + gradient(3 * jat + 0) += dgxij - dgxjk; + gradient(3 * jat + 1) += dgyij - dgyjk; + gradient(3 * jat + 2) += dgzij - dgzjk; + gradient(3 * kat + 0) += dgxik + dgxjk; + gradient(3 * kat + 1) += dgyik + dgyjk; + gradient(3 * kat + 2) += dgzik + dgzjk; + + dEdcn(iat) -= + e * 0.5 * (dc6dcn(iat, jat) / c6ij + dc6dcn(iat, kat) / c6ik); + dEdcn(jat) -= + e * 0.5 * (dc6dcn(jat, iat) / c6ij + dc6dcn(jat, kat) / c6jk); + dEdcn(kat) -= + e * 0.5 * (dc6dcn(kat, iat) / c6ik + dc6dcn(kat, jat) / c6jk); + + dEdq(iat) -= + e * 0.5 * (dc6dq(iat, jat) / c6ij + dc6dq(iat, kat) / c6ik); + dEdq(jat) -= + e * 0.5 * (dc6dq(jat, iat) / c6ij + dc6dq(jat, kat) / c6jk); + dEdq(kat) -= + e * 0.5 * (dc6dq(kat, iat) / c6ik + dc6dq(kat, jat) / c6jk); } } } @@ -312,4 +331,4 @@ double triple_scale(int ii, int jj, int kk) { return triple; } -} // namespace dftd4 +} // namespace dftd4 diff --git a/src/damping/dftd_rational.cpp b/src/damping/dftd_rational.cpp index d1c4de6..ac04cac 100644 --- a/src/damping/dftd_rational.cpp +++ b/src/damping/dftd_rational.cpp @@ -18,16 +18,15 @@ #include +#include "damping/dftd_atm.h" +#include "damping/dftd_rational.h" #include "dftd_cblas.h" -#include "dftd_eeq.h" #include "dftd_dispersion.h" +#include "dftd_eeq.h" #include "dftd_geometry.h" #include "dftd_matrix.h" #include "dftd_ncoord.h" #include "dftd_parameters.h" -#include "damping/dftd_atm.h" -#include "damping/dftd_rational.h" - namespace dftd4 { @@ -39,17 +38,17 @@ inline double fdmprdr_bj(const int n, const double r, const double c) { } int get_dispersion2( - const TMolecule& mol, - const TMatrix& dist, + const TMolecule &mol, + const TMatrix &dist, const double cutoff, - const dparam& par, - const TMatrix& c6, - const TMatrix& dc6dcn, - const TMatrix& dc6dq, - TVector& energy, - TVector& dEdcn, - TVector& dEdq, - TVector& gradient, + const dparam &par, + const TMatrix &c6, + const TMatrix &dc6dcn, + const TMatrix &dc6dq, + TVector &energy, + TVector &dEdcn, + TVector &dEdq, + TVector &gradient, bool lgrad /*= false*/ ) { int info{0}; @@ -66,12 +65,12 @@ int get_dispersion2( } int get_dispersion2_energy( - const TMolecule& mol, - const TMatrix& dist, + const TMolecule &mol, + const TMatrix &dist, const double cutoff, - const dparam& par, - const TMatrix& c6, - TVector& energy + const dparam &par, + const TMatrix &c6, + TVector &energy ) { int izp{0}, jzp{0}; double r{0.0}, r0ij{0.0}, r4r2ij{0.0}, c6ij{0.0}; @@ -84,7 +83,7 @@ int get_dispersion2_energy( jzp = mol.at(jat); r = dist(iat, jat); if (r > cutoff) continue; - + r4r2ij = 3.0 * r4r2[izp] * r4r2[jzp]; r0ij = par.a1 * sqrt(r4r2ij) + par.a2; c6ij = c6(iat, jat); @@ -99,7 +98,7 @@ int get_dispersion2_energy( edisp += par.s10 * 49.0 / 40.0 * pow(r4r2ij, 2) * t10; } - e = -c6ij*edisp * 0.5; + e = -c6ij * edisp * 0.5; energy(iat) += e; energy(jat) += e; } @@ -109,17 +108,17 @@ int get_dispersion2_energy( } int get_dispersion2_derivs( - const TMolecule& mol, - const TMatrix& dist, + const TMolecule &mol, + const TMatrix &dist, const double cutoff, - const dparam& par, - const TMatrix& c6, - const TMatrix& dc6dcn, - const TMatrix& dc6dq, - TVector& energy, - TVector& dEdcn, - TVector& dEdq, - TVector& gradient + const dparam &par, + const TMatrix &c6, + const TMatrix &dc6dcn, + const TMatrix &dc6dq, + TVector &energy, + TVector &dEdcn, + TVector &dEdq, + TVector &gradient ) { int izp{0}, jzp{0}; double r{0.0}, r0ij{0.0}, r4r2ij{0.0}, c6ij{0.0}; @@ -135,7 +134,7 @@ int get_dispersion2_derivs( jzp = mol.at(jat); r = dist(iat, jat); if (r > cutoff) continue; - + r4r2ij = 3.0 * r4r2[izp] * r4r2[jzp]; r0ij = par.a1 * sqrt(r4r2ij) + par.a2; c6ij = c6(iat, jat); @@ -156,7 +155,7 @@ int get_dispersion2_derivs( gdisp += par.s10 * 49.0 / 40.0 * pow(r4r2ij, 2) * dt10; } - e = -c6ij*edisp * 0.5; + e = -c6ij * edisp * 0.5; energy(iat) += e; energy(jat) += e; @@ -173,36 +172,48 @@ int get_dispersion2_derivs( dEdq(iat) -= dc6dq(iat, jat) * edisp; dEdq(jat) -= dc6dq(jat, iat) * edisp; - gradient(3*iat) += dgx; - gradient(3*iat + 1) += dgy; - gradient(3*iat + 2) += dgz; - gradient(3*jat) -= dgx; - gradient(3*jat + 1) -= dgy; - gradient(3*jat + 2) -= dgz; + gradient(3 * iat) += dgx; + gradient(3 * iat + 1) += dgy; + gradient(3 * iat + 2) += dgz; + gradient(3 * jat) -= dgx; + gradient(3 * jat + 1) -= dgy; + gradient(3 * jat + 2) -= dgz; } } return EXIT_SUCCESS; } - int get_dispersion3( - const TMolecule& mol, - const TMatrix& dist, + const TMolecule &mol, + const TMatrix &dist, const double cutoff, - const dparam& par, - const TMatrix& c6, - const TMatrix& dc6dcn, - const TMatrix& dc6dq, - TVector& energy, - TVector& dEdcn, - TVector& dEdq, - TVector& gradient, + const dparam &par, + const TMatrix &c6, + const TMatrix &dc6dcn, + const TMatrix &dc6dq, + TVector &energy, + TVector &dEdcn, + TVector &dEdq, + TVector &gradient, bool lgrad /*= false*/ ) { return get_atm_dispersion( - mol, dist, cutoff, par.s9, par.a1, par.a2, par.alp, c6, dc6dcn, dc6dq, - energy, dEdcn, dEdq, gradient, lgrad + mol, + dist, + cutoff, + par.s9, + par.a1, + par.a2, + par.alp, + c6, + dc6dcn, + dc6dq, + energy, + dEdcn, + dEdq, + gradient, + lgrad ); } -} // namespace dftd4 +} // namespace dftd4 diff --git a/test/test_ncoord.cpp b/test/test_ncoord.cpp index 4247bdd..31dd7bf 100644 --- a/test/test_ncoord.cpp +++ b/test/test_ncoord.cpp @@ -27,7 +27,6 @@ using namespace dftd4; - int test_cn( const int n, const char atoms[][4], @@ -49,7 +48,7 @@ int test_cn( } // distances - TMatrix dist; + TMatrix dist; dist.New(n, n); calc_distances(mol, dist); @@ -62,9 +61,7 @@ int test_cn( // compare to ref for (int i = 0; i != n; i++) { - if (check(cn(i), ref(i)) == EXIT_FAILURE) { - return EXIT_FAILURE; - } + if (check(cn(i), ref(i)) == EXIT_FAILURE) { return EXIT_FAILURE; } } return EXIT_SUCCESS; @@ -72,15 +69,14 @@ int test_cn( int test_ncoord() { int info; - + info = test_cn( mb16_43_01_n, mb16_43_01_atoms, mb16_43_01_coord, mb16_43_01_ref_cn ); if (info != EXIT_SUCCESS) return info; - info = test_cn( - rost61_m1_n, rost61_m1_atoms, rost61_m1_coord, rost61_m1_ref_cn - ); + info = + test_cn(rost61_m1_n, rost61_m1_atoms, rost61_m1_coord, rost61_m1_ref_cn); if (info != EXIT_SUCCESS) return info; return EXIT_SUCCESS; diff --git a/test/test_ncoord.h b/test/test_ncoord.h index 6a0fe70..ed235be 100644 --- a/test/test_ncoord.h +++ b/test/test_ncoord.h @@ -18,23 +18,33 @@ #ifndef TEST_NCOORD_H #define TEST_NCOORD_H -static const double mb16_43_01_ref_cn[16] { - +3.0734967711040162E+00, +9.3146160511610288E-01, +1.4370943937583889E+00, - +1.3330943158196045E+00, +7.2074352703033706E-01, +8.5965900477098189E-01, - +1.3578215817792136E+00, +1.5394000699602508E+00, +3.1940036819525854E+00, - +8.1216211163134222E-01, +8.5953344378485375E-01, +1.5334710815558710E+00, - +4.2331498952572053E+00, +3.0304850456739563E+00, +3.4522931948830609E+00, +static const double mb16_43_01_ref_cn[16]{ + +3.0734967711040162E+00, + +9.3146160511610288E-01, + +1.4370943937583889E+00, + +1.3330943158196045E+00, + +7.2074352703033706E-01, + +8.5965900477098189E-01, + +1.3578215817792136E+00, + +1.5394000699602508E+00, + +3.1940036819525854E+00, + +8.1216211163134222E-01, + +8.5953344378485375E-01, + +1.5334710815558710E+00, + +4.2331498952572053E+00, + +3.0304850456739563E+00, + +3.4522931948830609E+00, +4.2847828965226427E+00 }; -static const double rost61_m1_ref_cn[22] { - +3.7743890839860823E+00, +3.7533525980157583E+00, +3.7272213911500702E+00, - +9.2591719422049701E-01, +3.7403891707352228E+00, +3.7389943427981649E+00, - +9.2590197299758326E-01, +9.2599894955767703E-01, +9.2611826399971675E-01, - +9.2576484754843247E-01, +3.7170298017526591E+00, +3.7298933986021074E+00, - +3.7649883287219090E+00, +9.2586986693780782E-01, +3.7529905532199597E+00, - +3.7621775881176829E+00, +9.2587667173548660E-01, +9.2604216237640058E-01, - +9.2598322975544489E-01, +9.2588963692121407E-01, +8.1580076383956079E+00, +static const double rost61_m1_ref_cn[22]{ + +3.7743890839860823E+00, +3.7533525980157583E+00, +3.7272213911500702E+00, + +9.2591719422049701E-01, +3.7403891707352228E+00, +3.7389943427981649E+00, + +9.2590197299758326E-01, +9.2599894955767703E-01, +9.2611826399971675E-01, + +9.2576484754843247E-01, +3.7170298017526591E+00, +3.7298933986021074E+00, + +3.7649883287219090E+00, +9.2586986693780782E-01, +3.7529905532199597E+00, + +3.7621775881176829E+00, +9.2587667173548660E-01, +9.2604216237640058E-01, + +9.2598322975544489E-01, +9.2588963692121407E-01, +8.1580076383956079E+00, +8.4210289216027034E-01 }; diff --git a/test/test_param.h b/test/test_param.h index 612d143..6459e47 100644 --- a/test/test_param.h +++ b/test/test_param.h @@ -19,8 +19,8 @@ #define TEST_PARAM_H #include -#include #include +#include using namespace dftd4; @@ -29,158 +29,175 @@ extern int test_rational_damping(const double ref[], TCutoff cutoff); static const int nfuncs = 97; -static const char funcs[nfuncs][32] { - "hf", "b-lyp", "bpbe", "b-p", "bpw", "lb94", "mpwlyp", "mpwpw", - "o-lyp", "opbe", "pbe", "rpbe", "revpbe", "pw86pbe", "rpw86pbe", - "pw91", "pwp", "x-lyp", "b97", "tpss", "revtpss", "scan", "rscan", - "r2scan", "b1lyp", "b3-lyp", "bh-lyp", "b1p", "b3p", "b1pw", "b3pw", - "o3-lyp", "revpbe0", "revpbe38", "pbe0", "pwp1", "pw1pw", "mpw1pw", - "mpw1lyp", "pw6b95", "tpssh", "tpss0", "x3-lyp", "m06l", "m06", - "wb97", "wb97x", "cam-b3lyp", "lc-blyp", "lh07tsvwn", "lh07ssvwn", - "lh12ctssirpw92", "lh12ctssifpw92", "lh14tcalpbe", "lh20t", - "b2plyp", "b2gpplyp", "mpw2plyp", "pwpb95", "dsdblyp", "dsdpbe", - "dsdpbeb95", "dsdpbep86", "dsdsvwn", "dodblyp", "dodpbe", "dodpbeb95", - "dodpbep86", "dodsvwn", "pbe02", "pbe0dh", "dftb(3ob)", "dftb(mio)", - "dftb(pbc)", "dftb(matsci)", "dftb(ob2)", "b1b95", "glyp", "revpbe0dh", - "revtpss0", "revdsd-pbep86", "revdsd-pbe", "revdsd-blyp", - "revdod-pbep86", "b97m", "wb97m", "pbesol", "am05", "mn12sx", - "hse03", "hse06", "hse12", "hse12s", "hsesol", "r2scanh", "r2scan0", +static const char funcs[nfuncs][32]{ + "hf", + "b-lyp", + "bpbe", + "b-p", + "bpw", + "lb94", + "mpwlyp", + "mpwpw", + "o-lyp", + "opbe", + "pbe", + "rpbe", + "revpbe", + "pw86pbe", + "rpw86pbe", + "pw91", + "pwp", + "x-lyp", + "b97", + "tpss", + "revtpss", + "scan", + "rscan", + "r2scan", + "b1lyp", + "b3-lyp", + "bh-lyp", + "b1p", + "b3p", + "b1pw", + "b3pw", + "o3-lyp", + "revpbe0", + "revpbe38", + "pbe0", + "pwp1", + "pw1pw", + "mpw1pw", + "mpw1lyp", + "pw6b95", + "tpssh", + "tpss0", + "x3-lyp", + "m06l", + "m06", + "wb97", + "wb97x", + "cam-b3lyp", + "lc-blyp", + "lh07tsvwn", + "lh07ssvwn", + "lh12ctssirpw92", + "lh12ctssifpw92", + "lh14tcalpbe", + "lh20t", + "b2plyp", + "b2gpplyp", + "mpw2plyp", + "pwpb95", + "dsdblyp", + "dsdpbe", + "dsdpbeb95", + "dsdpbep86", + "dsdsvwn", + "dodblyp", + "dodpbe", + "dodpbeb95", + "dodpbep86", + "dodsvwn", + "pbe02", + "pbe0dh", + "dftb(3ob)", + "dftb(mio)", + "dftb(pbc)", + "dftb(matsci)", + "dftb(ob2)", + "b1b95", + "glyp", + "revpbe0dh", + "revtpss0", + "revdsd-pbep86", + "revdsd-pbe", + "revdsd-blyp", + "revdod-pbep86", + "b97m", + "wb97m", + "pbesol", + "am05", + "mn12sx", + "hse03", + "hse06", + "hse12", + "hse12s", + "hsesol", + "r2scanh", + "r2scan0", "r2scan50" }; -static const double ref[nfuncs] { - -2.82477942248756E-1,-1.83931931447376E-1,-1.67883731806544E-1, - -1.19260344286630E-1,-1.73553643190541E-1,-4.64187008858948E-1, - -1.33350631907425E-1,-1.29120540723649E-1,-4.23162977083251E-1, - -4.19510900741786E-1,-8.78277686099131E-2,-2.88574075471549E-1, - -2.59731166101835E-1,-9.87868654968457E-2,-9.80588487018907E-2, - -7.25055144755650E-2,-3.38536835705987E-2,-1.99101153940959E-1, - -1.51912675830051E-1,-1.17214547334950E-1,-9.43596369853342E-2, - -1.84970291718116E-2,-3.25233366958794E-2,-2.79219862478753E-2, - -1.41170942097673E-1,-1.36036600333643E-1,-9.86753151279879E-2, - -8.99492838831500E-2,-9.44139521417514E-2,-1.35742320110483E-1, - -1.32216887418510E-1,-1.14737686310415E-1,-1.86535361600070E-1, - -1.46887478462845E-1,-7.64961638449203E-2,-3.19815917382690E-2, - -6.61290129604851E-2,-1.05230600646517E-1,-1.02669099979520E-1, - -7.17939219769430E-2,-1.08455784745289E-1,-1.01080258057282E-1, - -1.12546242873133E-1,-1.23576092386848E-2,-2.00031672761610E-2, - -8.37504026938291E-3,-1.58515232511485E-2,-7.77668954242432E-2, - -1.29932494116969E-2,-2.07804101611958E-1,-2.62377725485470E-1, - -3.39614225692246E-1,-3.74698166577280E-1,-8.84516853613503E-2, - -5.21275383612525E-2,-6.24504147646583E-2,-3.94277424205698E-2, - -4.14986337497841E-2,-6.11949137677711E-2,-3.85836295170037E-2, - -4.77329012547881E-2,-4.36022695279866E-2,-1.90151394463137E-2, - -2.31910926279841E-2,-8.51189313158795E-2,-6.56921629858802E-2, - -5.96170924208417E-2,-5.17434304314999E-2,-3.15736132524860E-2, - -8.31688591017636E-3,-4.85454964254385E-2,-6.17262146525898E-2, - -4.55955214285284E-2,-5.85633740272026E-2,-7.13638868029425E-2, - -4.28867615756140E-2,-1.03264576995975E-1,-2.71073699938782E-1, - -1.02702213161577E-1,-8.21356031484763E-2,-5.65853020266326E-2, - -8.56296234267241E-2,-8.91479948664164E-2,-6.13524147568788E-2, - -1.24018192547142E-1,-1.05459213220861E-1,-3.62056543364222E-2, - -3.62056559678757E-2,-2.40058302367699E-2,-6.96544052724497E-2, - -7.14846527169669E-2,-6.94907081073225E-2,-6.20062456876537E-2, - -5.03611950945288E-2,-2.92623318488036E-2,-3.16651212523792E-2, +static const double ref[nfuncs]{ + -2.82477942248756E-1, -1.83931931447376E-1, -1.67883731806544E-1, + -1.19260344286630E-1, -1.73553643190541E-1, -4.64187008858948E-1, + -1.33350631907425E-1, -1.29120540723649E-1, -4.23162977083251E-1, + -4.19510900741786E-1, -8.78277686099131E-2, -2.88574075471549E-1, + -2.59731166101835E-1, -9.87868654968457E-2, -9.80588487018907E-2, + -7.25055144755650E-2, -3.38536835705987E-2, -1.99101153940959E-1, + -1.51912675830051E-1, -1.17214547334950E-1, -9.43596369853342E-2, + -1.84970291718116E-2, -3.25233366958794E-2, -2.79219862478753E-2, + -1.41170942097673E-1, -1.36036600333643E-1, -9.86753151279879E-2, + -8.99492838831500E-2, -9.44139521417514E-2, -1.35742320110483E-1, + -1.32216887418510E-1, -1.14737686310415E-1, -1.86535361600070E-1, + -1.46887478462845E-1, -7.64961638449203E-2, -3.19815917382690E-2, + -6.61290129604851E-2, -1.05230600646517E-1, -1.02669099979520E-1, + -7.17939219769430E-2, -1.08455784745289E-1, -1.01080258057282E-1, + -1.12546242873133E-1, -1.23576092386848E-2, -2.00031672761610E-2, + -8.37504026938291E-3, -1.58515232511485E-2, -7.77668954242432E-2, + -1.29932494116969E-2, -2.07804101611958E-1, -2.62377725485470E-1, + -3.39614225692246E-1, -3.74698166577280E-1, -8.84516853613503E-2, + -5.21275383612525E-2, -6.24504147646583E-2, -3.94277424205698E-2, + -4.14986337497841E-2, -6.11949137677711E-2, -3.85836295170037E-2, + -4.77329012547881E-2, -4.36022695279866E-2, -1.90151394463137E-2, + -2.31910926279841E-2, -8.51189313158795E-2, -6.56921629858802E-2, + -5.96170924208417E-2, -5.17434304314999E-2, -3.15736132524860E-2, + -8.31688591017636E-3, -4.85454964254385E-2, -6.17262146525898E-2, + -4.55955214285284E-2, -5.85633740272026E-2, -7.13638868029425E-2, + -4.28867615756140E-2, -1.03264576995975E-1, -2.71073699938782E-1, + -1.02702213161577E-1, -8.21356031484763E-2, -5.65853020266326E-2, + -8.56296234267241E-2, -8.91479948664164E-2, -6.13524147568788E-2, + -1.24018192547142E-1, -1.05459213220861E-1, -3.62056543364222E-2, + -3.62056559678757E-2, -2.40058302367699E-2, -6.96544052724497E-2, + -7.14846527169669E-2, -6.94907081073225E-2, -6.20062456876537E-2, + -5.03611950945288E-2, -2.92623318488036E-2, -3.16651212523792E-2, -3.45136899744303E-2 }; -static const double ref_no_cutoff[nfuncs] { - -2.8259959590382366E-01, - -1.8404742721234266E-01, - -1.6798937024195304E-01, - -1.1935355217292869E-01, - -1.7366283799329513E-01, - -4.6430923163432714E-01, - -1.3346580384141335E-01, - -1.2923057292683821E-01, - -4.2328482522988720E-01, - -4.1963240189378620E-01, - -8.7933519082332240E-02, - -2.8869613111622339E-01, - -2.5985221342803455E-01, - -9.8893933408153903E-02, - -9.8166674429422632E-02, - -7.2605591514865300E-02, - -3.3915936792186512E-02, - -1.9922033814494208E-01, - -1.5203207292642029E-01, - -1.1732196176499123E-01, - -9.4460953697244290E-02, - -1.8508178547666513E-02, - -3.2572027285843817E-02, - -2.7965656650905928E-02, - -1.4128237865557192E-01, - -1.3614663010332351E-01, - -9.8776916779818535E-02, - -9.0027422431692011E-02, - -9.4497497790429527E-02, - -1.3584429225597902E-01, - -1.3231891834258791E-01, - -1.1484369432188699E-01, - -1.8665428267713263E-01, - -1.4700208605535828E-01, - -7.6591931831600193E-02, - -3.2034113023795184E-02, - -6.6220768631983964E-02, - -1.0533339468079540E-01, - -1.0277805122805916E-01, - -7.1912075732403241E-02, - -1.0855929194374825E-01, - -1.0118342099540853E-01, - -1.1265405513604643E-01, - -1.2361163873034551E-02, - -2.0033492687049504E-02, - -8.3715876932202084E-03, - -1.5871986581992004E-02, - -7.7856722692803276E-02, - -1.2993089916236900E-02, - -2.0792249021524978E-01, - -2.6249584973626233E-01, - -3.3973560197631031E-01, - -3.7481974145720043E-01, - -8.8553324335207614E-02, - -5.2223863418124516E-02, - -6.2549827081499740E-02, - -3.9511828658811962E-02, - -4.1587075619746859E-02, - -6.1315118738896508E-02, - -3.8676661542281436E-02, - -4.7838962519922201E-02, - -4.3721931993326975E-02, - -1.9079922103516497E-02, - -2.3261054412440016E-02, - -8.5231997813884242E-02, - -6.5803424344231073E-02, - -5.9737404157675297E-02, - -5.1851768992444328E-02, - -3.1654165615043431E-02, - -8.3257937599140481E-03, - -4.8624385696080269E-02, - -6.1819725698154936E-02, - -4.5658803921709075E-02, - -5.8629726042068023E-02, - -7.1429204290326231E-02, - -4.2926220854554435E-02, - -1.0337253252219039E-01, - -2.7118928551625038E-01, - -1.0281292290726875E-01, - -8.2229990650757728E-02, - -5.6706096227474322E-02, - -8.5751286380770605E-02, - -8.9270144738356177E-02, - -6.1473208957720508E-02, - -1.2413617396073341E-01, - -1.0557238186347086E-01, - -3.6248687079787421E-02, - -3.6248688696273144E-02, - -2.4036121627224290E-02, - -6.9745139211587009E-02, - -7.1576872338838152E-02, - -6.9580711257046970E-02, - -6.2089574187874155E-02, - -5.0422920264847002E-02, - -2.9305044225664403E-02, - -3.1711738228039008E-02, +static const double ref_no_cutoff[nfuncs]{ + -2.8259959590382366E-01, -1.8404742721234266E-01, -1.6798937024195304E-01, + -1.1935355217292869E-01, -1.7366283799329513E-01, -4.6430923163432714E-01, + -1.3346580384141335E-01, -1.2923057292683821E-01, -4.2328482522988720E-01, + -4.1963240189378620E-01, -8.7933519082332240E-02, -2.8869613111622339E-01, + -2.5985221342803455E-01, -9.8893933408153903E-02, -9.8166674429422632E-02, + -7.2605591514865300E-02, -3.3915936792186512E-02, -1.9922033814494208E-01, + -1.5203207292642029E-01, -1.1732196176499123E-01, -9.4460953697244290E-02, + -1.8508178547666513E-02, -3.2572027285843817E-02, -2.7965656650905928E-02, + -1.4128237865557192E-01, -1.3614663010332351E-01, -9.8776916779818535E-02, + -9.0027422431692011E-02, -9.4497497790429527E-02, -1.3584429225597902E-01, + -1.3231891834258791E-01, -1.1484369432188699E-01, -1.8665428267713263E-01, + -1.4700208605535828E-01, -7.6591931831600193E-02, -3.2034113023795184E-02, + -6.6220768631983964E-02, -1.0533339468079540E-01, -1.0277805122805916E-01, + -7.1912075732403241E-02, -1.0855929194374825E-01, -1.0118342099540853E-01, + -1.1265405513604643E-01, -1.2361163873034551E-02, -2.0033492687049504E-02, + -8.3715876932202084E-03, -1.5871986581992004E-02, -7.7856722692803276E-02, + -1.2993089916236900E-02, -2.0792249021524978E-01, -2.6249584973626233E-01, + -3.3973560197631031E-01, -3.7481974145720043E-01, -8.8553324335207614E-02, + -5.2223863418124516E-02, -6.2549827081499740E-02, -3.9511828658811962E-02, + -4.1587075619746859E-02, -6.1315118738896508E-02, -3.8676661542281436E-02, + -4.7838962519922201E-02, -4.3721931993326975E-02, -1.9079922103516497E-02, + -2.3261054412440016E-02, -8.5231997813884242E-02, -6.5803424344231073E-02, + -5.9737404157675297E-02, -5.1851768992444328E-02, -3.1654165615043431E-02, + -8.3257937599140481E-03, -4.8624385696080269E-02, -6.1819725698154936E-02, + -4.5658803921709075E-02, -5.8629726042068023E-02, -7.1429204290326231E-02, + -4.2926220854554435E-02, -1.0337253252219039E-01, -2.7118928551625038E-01, + -1.0281292290726875E-01, -8.2229990650757728E-02, -5.6706096227474322E-02, + -8.5751286380770605E-02, -8.9270144738356177E-02, -6.1473208957720508E-02, + -1.2413617396073341E-01, -1.0557238186347086E-01, -3.6248687079787421E-02, + -3.6248688696273144E-02, -2.4036121627224290E-02, -6.9745139211587009E-02, + -7.1576872338838152E-02, -6.9580711257046970E-02, -6.2089574187874155E-02, + -5.0422920264847002E-02, -2.9305044225664403E-02, -3.1711738228039008E-02, -3.4563360149465587E-02, }; diff --git a/test/util.cpp b/test/util.cpp index 5843622..4b46d03 100644 --- a/test/util.cpp +++ b/test/util.cpp @@ -27,14 +27,18 @@ using namespace dftd4; - -int get_molecule(int n, const char atoms[][4], const double coord[], TMolecule& mol) { +int get_molecule( + int n, + const char atoms[][4], + const double coord[], + TMolecule &mol +) { ; mol.GetMemory(n); for (int i = 0; i != n; i++) { - mol.xyz(i, 0) = coord[3*i]; - mol.xyz(i, 1) = coord[3*i+1]; - mol.xyz(i, 2) = coord[3*i+2]; + mol.xyz(i, 0) = coord[3 * i]; + mol.xyz(i, 1) = coord[3 * i + 1]; + mol.xyz(i, 2) = coord[3 * i + 2]; mol.at(i) = element(atoms[i]); } @@ -55,9 +59,7 @@ bool check( diff = std::fabs(actual - expected); } - if (diff > epsilon) { - return EXIT_FAILURE; - } + if (diff > epsilon) { return EXIT_FAILURE; } return EXIT_SUCCESS; }; @@ -75,13 +77,10 @@ bool check( diff = std::fabs(actual - expected); } - if (diff > epsilon) { - return EXIT_FAILURE; - } + if (diff > epsilon) { return EXIT_FAILURE; } return EXIT_SUCCESS; }; - void print_fail(const char specifier[32], double obtained, double expected) { printf("Failed for: '%s'\n", specifier); printf(" - Expected %.16f\n", expected); @@ -89,20 +88,19 @@ void print_fail(const char specifier[32], double obtained, double expected) { printf(" - Diff (abs) %.16f\n", fabs(expected - obtained)); } - -int element(const std::string& sym) { +int element(const std::string &sym) { char elem[3]{" "}; char pse[118][3]{ - "h ", "he", "li", "be", "b ", "c ", "n ", "o ", "f ", "ne", "na", "mg", - "al", "si", "p ", "s ", "cl", "ar", "k ", "ca", "sc", "ti", "v ", "cr", - "mn", "fe", "co", "ni", "cu", "zn", "ga", "ge", "as", "se", "br", "kr", - "rb", "sr", "y ", "zr", "nb", "mo", "tc", "ru", "rh", "pd", "ag", "cd", - "in", "sn", "sb", "te", "i ", "xe", "cs", "ba", "la", "ce", "pr", "nd", - "pm", "sm", "eu", "gd", "tb", "dy", "ho", "er", "tm", "yb", "lu", "hf", - "ta", "w ", "re", "os", "ir", "pt", "au", "hg", "tl", "pb", "bi", "po", - "at", "rn", "fr", "ra", "ac", "th", "pa", "u ", "np", "pu", "am", "cm", - "bk", "cf", "es", "fm", "md", "no", "lr", "rf", "db", "sg", "bh", "hs", - "mt", "ds", "rg", "cn", "nh", "fl", "mc", "lv", "ts", "og", + "h ", "he", "li", "be", "b ", "c ", "n ", "o ", "f ", "ne", "na", "mg", + "al", "si", "p ", "s ", "cl", "ar", "k ", "ca", "sc", "ti", "v ", "cr", + "mn", "fe", "co", "ni", "cu", "zn", "ga", "ge", "as", "se", "br", "kr", + "rb", "sr", "y ", "zr", "nb", "mo", "tc", "ru", "rh", "pd", "ag", "cd", + "in", "sn", "sb", "te", "i ", "xe", "cs", "ba", "la", "ce", "pr", "nd", + "pm", "sm", "eu", "gd", "tb", "dy", "ho", "er", "tm", "yb", "lu", "hf", + "ta", "w ", "re", "os", "ir", "pt", "au", "hg", "tl", "pb", "bi", "po", + "at", "rn", "fr", "ra", "ac", "th", "pa", "u ", "np", "pu", "am", "cm", + "bk", "cf", "es", "fm", "md", "no", "lr", "rf", "db", "sg", "bh", "hs", + "mt", "ds", "rg", "cn", "nh", "fl", "mc", "lv", "ts", "og", }; std::transform(sym.begin(), sym.end(), elem, ::tolower); diff --git a/test/util.h b/test/util.h index 8d476c7..4f996dc 100644 --- a/test/util.h +++ b/test/util.h @@ -22,14 +22,25 @@ #include - -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); - -extern void print_fail(const char specifier[32], double obtained, double expected); - -extern int element(const std::string& sym); +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); + +extern void + print_fail(const char specifier[32], double obtained, double expected); + +extern int element(const std::string &sym); #endif // UTIL_H