diff --git a/include/dftd_cblas.h b/include/dftd_cblas.h index 1513de8..e71008a 100644 --- a/include/dftd_cblas.h +++ b/include/dftd_cblas.h @@ -35,8 +35,8 @@ inline int BLAS_Add_Mat_x_Vec( bool Transpose, double alpha ) { - if (A.rows == C.N && A.cols == V.N) { if (Transpose) { + if (A.cols == C.N && A.rows == V.N) { cblas_dgemv( CblasRowMajor, CblasTrans, @@ -52,7 +52,9 @@ inline int BLAS_Add_Mat_x_Vec( 1 ); return EXIT_SUCCESS; - } else { + }; + } else { + if (A.rows == C.N && A.cols == V.N) { cblas_dgemv( CblasRowMajor, CblasNoTrans, diff --git a/include/dftd_ncoord.h b/include/dftd_ncoord.h index 326382c..c174b16 100644 --- a/include/dftd_ncoord.h +++ b/include/dftd_ncoord.h @@ -28,6 +28,143 @@ namespace dftd4 { + +class NCoordBase +{ + public: + TVector cn; + TMatrix dcndr; + static const double rad[]; + double kcn; // Steepness of counting function + double norm_exp; + double cutoff; + /** + * Wrapper for coordination number calculation. + * + * @param mol Molecule object. + * @param dist Distance matrix. + * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}). + * @param lgrad Flag for gradient computation. + * @return Exit status. + */ + int get_ncoord( // with ghost atoms + const TMolecule &mol, + const TMatrix &dist, + const double cutoff, + bool lgrad); + /** + * Wrapper for coordination number calculation. + * + * @param mol Molecule object. + * @param realIdx List for real atoms excluding ghost/non atoms. + * @param dist Distance matrix. + * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}). + * @param cn Vector of coordination numbers. + * @return Exit status. + */ + int get_ncoord( // without ghost atoms + const TMolecule &mol, + const TIVector &realIdx, + const TMatrix &dist, + bool lgrad); + /** + * Calculate the coordination number. + * + * @param mol Molecule object. + * @param realIdx List for real atoms excluding ghost/non atoms. + * @param dist Distance matrix. + * @return Exit status. + */ + int ncoord_base( + const TMolecule &mol, + const TIVector &realIdx, + const TMatrix &dist); + /** + * Calculate error function coordination number and derivative + * w.r.t. nuclear coordinates. + * + * @param mol Molecule object. + * @param realIdx List for real atoms excluding ghost/non atoms. + * @param dist Distance matrix. + * @return Exit status. + */ + int dncoord_base( + const TMolecule &mol, + const TIVector &realIdx, + const TMatrix &dist); + /** + * Electronegativity factor. + * + * @param i atom index + * @param j atom index + * @return Value of the electronegativity factor. + */ + virtual double get_en_factor(int i, int j) const; + /** + * base class function for coordination number contributions. + * + * @param rr atomic distance over covalent radii + * @return Value of the counting function. + */ + virtual double count_fct(double rr) const = 0; + /** + * Derivative of the counting function w.r.t. the distance. + * + * @param rr atomic distance over covalent radii + * @return Derivative of the counting function. + */ + virtual double dcount_fct(double rr) const = 0; + /** + * TCutoff function for large coordination numbers + * + * @param cn_max Maximum CN (not strictly obeyed). + * @param cn On input coordination number, on output modified CN. + * @param dcndr On input derivative of CN w.r.t. cartesian coordinates, + * on output derivative of modified CN. + * @param lgrad Flag for gradient calculation. + */ + virtual int cut_coordination_number(const double cn_max, TVector &cn, TMatrix &dcndr, bool lgrad) = 0; + // Constructor + NCoordBase(double optional_kcn = 7.5, double optional_norm_exp = 1.0, double optional_cutoff = 25.0); + // Virtual destructor + virtual ~NCoordBase() { + cn.DelVec(); + dcndr.DelMat(); + } +}; + +class NCoordErf : public NCoordBase { + public: + // erf() based counting function + double count_fct(double) const override; + // derivative of the erf() based counting function + double dcount_fct(double) const override; + // Soft maximum/cutoff for coordination number + int cut_coordination_number(const double, TVector&, TMatrix&, bool) override; + // Constructor + NCoordErf(double optional_kcn = 7.5, double optional_norm_exp = 1.0, double optional_cutoff = 25.0) + : NCoordBase(optional_kcn, optional_norm_exp, optional_cutoff){} + // Use default destructor; base class handles cleanup + ~NCoordErf() override = default; +}; + +class NCoordErfD4 : public NCoordBase { + public: + // erf() based counting function + double count_fct(double) const override; + // derivative of the erf() based counting function + double dcount_fct(double) const override; + // coordination number scaling factor based on electronegativity difference + double get_en_factor(int, int) const override; + // Soft maximum/cutoff for coordination number + int cut_coordination_number(const double, TVector&, TMatrix&, bool) override; + // Constructor + NCoordErfD4(double optional_kcn = 7.5, double optional_norm_exp = 1.0, double optional_cutoff = 25.0) + : NCoordBase(optional_kcn, optional_norm_exp, optional_cutoff){} + // Use default destructor; base class handles cleanup + ~NCoordErfD4() override = default; +}; + /** * Calculate all distance pairs and store in matrix. * @@ -51,207 +188,6 @@ extern int calc_distances( */ extern void initializeRealIdx(int nat, TVector &realIdx); -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// - -/** - * Wrapper for error function coordination number. - * - * @param mol Molecule object. - * @param realIdx List for real atoms excluding ghost/non atoms. - * @param dist Distance matrix. - * @param cn Vector of coordination numbers. - * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}). - * @param dcndr Derivative of coordination number. - * @param lgrad Flag for gradient computation. - * @return Exit status. - */ -extern int get_ncoord_erf( - const TMolecule &mol, - const TIVector &realIdx, - const TMatrix &dist, - double cutoff, - TVector &cn, - TMatrix &dcndr, - bool lgrad = false -); - -/** - * Wrapper for error function coordination number. - * - * @param mol Molecule object. - * @param dist Distance matrix. - * @param cn Vector of coordination numbers. - * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}). - * @param dcndr Derivative of coordination number. - * @param lgrad Flag for gradient computation. - * @return Exit status. - */ -extern int get_ncoord_erf( - const TMolecule &mol, - const TMatrix &dist, - double cutoff, - TVector &cn, - TMatrix &dcndr, - bool lgrad = false -); - -/** - * Calculate error function coordination number. - * - * @param mol Molecule object. - * @param realIdx List for real atoms excluding ghost/non atoms. - * @param dist Distance matrix. - * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}). - * @param cn Vector of coordination numbers. - * @return Exit status. - */ -extern int ncoord_erf( - const TMolecule &mol, - const TIVector &realIdx, - const TMatrix &dist, - double cutoff, - TVector &cn -); - -/** - * Calculate error function coordination number and derivative - * w.r.t. nuclear coordinates. - * - * @param mol Molecule object. - * @param realIdx List for real atoms excluding ghost/non atoms. - * @param dist Distance matrix. - * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}). - * @param cn Vector of coordination numbers. - * @param dcndr Derivative of coordination number. - * @return Exit status. - */ -extern int dncoord_erf( - const TMolecule &mol, - const TIVector &realIdx, - const TMatrix &dist, - double cutoff, - TVector &cn, - TMatrix &dcndr -); - -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// - -/** - * Wrapper for error function coordination number for DFT-D4. - * - * @param mol Molecule object. - * @param realIdx List for real atoms excluding ghost/non atoms. - * @param dist Distance matrix. - * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}). - * @param cn Vector of coordination numbers. - * @return Exit status. - */ -extern int get_ncoord_d4( - const TMolecule &mol, - const TIVector &realIdx, - const TMatrix &dist, - double cutoff, - TVector &cn, - TMatrix &dcndr, - bool lgrad = false -); - -/** - * Wrapper for error function coordination number for DFT-D4. - * - * @param mol Molecule object. - * @param dist Distance matrix. - * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}). - * @param cn Vector of coordination numbers. - * @return Exit status. - */ -extern int get_ncoord_d4( - const TMolecule &mol, - const TMatrix &dist, - double cutoff, - TVector &cn, - TMatrix &dcndr, - bool lgrad = false -); - -/** - * Calculate covalent coordination number for DFT-D4. - * - * @param mol Molecule object. - * @param realIdx List for real atoms excluding ghost/non atoms. - * @param dist Distance matrix. - * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}). - * @param cn Vector of coordination numbers. - * @return Exit status. - */ -extern int ncoord_d4( - const TMolecule &mol, - const TIVector &realIdx, - const TMatrix &dist, - double cutoff, - TVector &cn -); - -/** - * Calculate covalent coordination number and derivative - * w.r.t. nuclear coordinates - * - * @param mol Molecule object. - * @param realIdx List for real atoms excluding ghost/non atoms. - * @param dist Distance matrix. - * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}). - * @param cn Vector of coordination numbers. - * @param dcndr Derivative of coordination number. - * @return Exit status. - */ -extern int dncoord_d4( - const TMolecule &mol, - const TIVector &realIdx, - const TMatrix &dist, - double cutoff, - TVector &cn, - TMatrix &dcndr -); - -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// - -/** - * Error function counting function for coordination number contributions. - * - * @param k Steepness of the counting function. - * @param rr Current distance over cutoff radius. - * @return Value of the counting function. - */ -extern double erf_count(double k, double rr); - -/** - * Derivative of the counting function w.r.t. the distance. - * - * @param k Steepness of the counting function. - * @param rr TCutoff radius. - * @return Derivative of the counting function. - */ -extern double derf_count(double k, double rr); - -/** - * TCutoff function for large coordination numbers - * - * @param cn_max Maximum CN (not strictly obeyed). - * @param cn On input coordination number, on output modified CN. - * @param dcndr On input derivative of CN w.r.t. cartesian coordinates, - * on output derivative of modified CN. - * @param lgrad Flag for gradient. - */ -extern int cut_coordination_number( - double cn_max, - TVector &cn, - TMatrix &dcndr, - bool lgrad -); - extern inline double log_cn_cut(double cn_max, double cn); extern inline double dlog_cn_cut(double cn_max, double cn); diff --git a/src/dftd_dispersion.cpp b/src/dftd_dispersion.cpp index bfbea23..7ad0d47 100644 --- a/src/dftd_dispersion.cpp +++ b/src/dftd_dispersion.cpp @@ -102,17 +102,14 @@ int get_dispersion( dist.NewMatrix(nat, nat); calc_distances(mol, realIdx, dist); - TVector cn; // D4 coordination number TVector q; // partial charges from EEQ model - TMatrix dcndr; // derivative of D4-CN TMatrix dqdr; // derivative of partial charges TVector gradient; // derivative of dispersion energy + NCoordErfD4 ncoord_erf; multicharge::EEQModel chrg_model; // Charge model - cn.NewVector(nat); q.NewVector(nat); if (lgrad) { - dcndr.NewMatrix(3 * nat, nat); dqdr.NewMatrix(3 * nat, nat); gradient.NewVector(3 * nat); } @@ -122,7 +119,7 @@ int get_dispersion( if (info != EXIT_SUCCESS) return info; // get the D4 coordination number - info = get_ncoord_d4(mol, realIdx, dist, cutoff.cn, cn, dcndr, lgrad); + info = ncoord_erf.get_ncoord(mol, realIdx, dist, lgrad); if (info != EXIT_SUCCESS) return info; // maximum number of reference systems @@ -145,7 +142,7 @@ int get_dispersion( dgwdq.NewMatrix(mref, nat); } info = d4.weight_references( - mol, realIdx, cn, q, refq, gwvec, dgwdcn, dgwdq, lgrad + mol, realIdx, ncoord_erf.cn, q, refq, gwvec, dgwdcn, dgwdq, lgrad ); if (info != EXIT_SUCCESS) return info; @@ -212,11 +209,10 @@ int get_dispersion( dgwdq.NewMatrix(mref, nat); } info = d4.weight_references( - mol, realIdx, cn, q, refq, gwvec, dgwdcn, dgwdq, lgrad + mol, realIdx, ncoord_erf.cn, q, refq, gwvec, dgwdcn, dgwdq, lgrad ); if (info != EXIT_SUCCESS) return info; - cn.Delete(); q.Delete(); refq.Delete(); @@ -253,7 +249,6 @@ int get_dispersion( ); if (info != EXIT_SUCCESS) return info; } else { - cn.Delete(); q.Delete(); refq.Delete(); gwvec.Delete(); @@ -266,9 +261,8 @@ int get_dispersion( dc6dcn.DelMat(); dc6dq.DelMat(); - if (lgrad) { BLAS_Add_Mat_x_Vec(gradient, dcndr, dEdcn, false, 1.0); } + if (lgrad) BLAS_Add_Mat_x_Vec(gradient, ncoord_erf.dcndr, dEdcn, true, 1.0); - dcndr.DelMat(); dEdcn.DelVec(); dEdq.DelVec(); diff --git a/src/dftd_eeq.cpp b/src/dftd_eeq.cpp index 2f6c435..22636eb 100644 --- a/src/dftd_eeq.cpp +++ b/src/dftd_eeq.cpp @@ -72,24 +72,17 @@ int ChargeModel::get_charges( bool lverbose{false}; int nat = realIdx.Max() + 1; - TVector cn; // EEQ cordination number - TMatrix dcndr; // Derivative of EEQ-CN - - cn.NewVec(nat); - if (lgrad) dcndr.NewMat(nat, 3 * nat); + dftd4::NCoordErf ncoord_erf; // get the EEQ coordination number - info = get_ncoord_erf(mol, realIdx, dist, cutoff, cn, dcndr, lgrad); + info = ncoord_erf.get_ncoord(mol, realIdx, dist, lgrad); if (info != EXIT_SUCCESS) return info; // corresponds to model%solve in Fortran info = - eeq_chrgeq(mol, realIdx, dist, charge, cn, q, dcndr, dqdr, lgrad, lverbose); + eeq_chrgeq(mol, realIdx, dist, charge, ncoord_erf.cn, q, ncoord_erf.dcndr, dqdr, lgrad, lverbose); if (info != EXIT_SUCCESS) return info; - dcndr.DelMat(); - cn.DelVec(); - return EXIT_SUCCESS; }; @@ -194,9 +187,9 @@ int ChargeModel::eeq_chrgeq( jj = realIdx(j); if (jj < 0) continue; - dAmat(3 * jj, ii) -= dcndr(jj, 3 * ii) * dxdcn(ii); - dAmat(3 * jj + 1, ii) -= dcndr(jj, 3 * ii + 1) * dxdcn(ii); - dAmat(3 * jj + 2, ii) -= dcndr(jj, 3 * ii + 2) * dxdcn(ii); + dAmat(3 * jj, ii) -= dcndr(ii, 3 * jj ) * dxdcn(ii); + dAmat(3 * jj + 1, ii) -= dcndr(ii, 3 * jj + 1) * dxdcn(ii); + dAmat(3 * jj + 2, ii) -= dcndr(ii, 3 * jj + 2) * dxdcn(ii); } } @@ -364,10 +357,10 @@ int EEQModel::get_damat_0d( atrace(jj, 1) -= dgy * q(ii); atrace(jj, 2) -= dgz * q(ii); - dAmat(3 * ii, jj) = dgx * q(ii); - dAmat(3 * ii + 1, jj) = dgy * q(ii); - dAmat(3 * ii + 2, jj) = dgz * q(ii); - dAmat(3 * jj, ii) = -dgx * q(jj); + dAmat(3 * ii, jj) = dgx * q(ii); + dAmat(3 * ii + 1, jj) = dgy * q(ii); + dAmat(3 * ii + 2, jj) = dgz * q(ii); + dAmat(3 * jj, ii) = -dgx * q(jj); dAmat(3 * jj + 1, ii) = -dgy * q(jj); dAmat(3 * jj + 2, ii) = -dgz * q(jj); } diff --git a/src/dftd_ncoord.cpp b/src/dftd_ncoord.cpp index a6cf431..8970717 100644 --- a/src/dftd_ncoord.cpp +++ b/src/dftd_ncoord.cpp @@ -32,6 +32,14 @@ namespace dftd4 { +NCoordBase::NCoordBase( + double optional_kcn, // defaults for D4 EEQ + double optional_norm_exp, // defaults for D4 EEQ + double optional_cutoff // defaults for D4 EEQ +) : kcn(optional_kcn), norm_exp(optional_norm_exp), cutoff(optional_cutoff) +{} + + /** * Covalent radii (taken from Pyykko and Atsumi, Chem. Eur. J. 15, 2009, * 188-197), values for metals decreased by 10 %. @@ -69,7 +77,7 @@ static const double covalent_rad_d3[119]{ * static const double autoaa = 0.52917721090449243; * static const double aatoau = 1.0 / autoaa; */ -static double rad[119]{ +const double NCoordBase::rad[119]{ 0.00000000000000, 0.80628314650472, 1.15903202310054, 3.02356179939270, 2.36845674285762, 1.94011882127699, 1.88972612462044, 1.78894073130735, 1.58736994468117, 1.61256629300944, 1.68815533799426, 3.52748876595816, @@ -274,35 +282,30 @@ void initializeRealIdx(int nat, TVector &realIdx) { /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// -int get_ncoord_erf( +int NCoordBase::get_ncoord( const TMolecule &mol, const TMatrix &dist, const double cutoff, - TVector &cn, - TMatrix &dcndr, bool lgrad ) { TIVector realIdx; initializeRealIdx(mol.NAtoms, realIdx); - return get_ncoord_erf(mol, realIdx, dist, cutoff, cn, dcndr, lgrad); + return get_ncoord(mol, realIdx, dist, lgrad); }; -int get_ncoord_erf( +int NCoordBase::get_ncoord( const TMolecule &mol, const TIVector &realIdx, const TMatrix &dist, - const double cutoff, - TVector &cn, - TMatrix &dcndr, bool lgrad ) { int info; if (lgrad) { - info = dncoord_erf(mol, realIdx, dist, cutoff, cn, dcndr); + info = dncoord_base(mol, realIdx, dist); } else { - info = ncoord_erf(mol, realIdx, dist, cutoff, cn); + info = ncoord_base(mol, realIdx, dist); } if (info != EXIT_SUCCESS) return info; @@ -312,15 +315,17 @@ int get_ncoord_erf( return EXIT_SUCCESS; }; -int ncoord_erf( +int NCoordBase::ncoord_base( const TMolecule &mol, const TIVector &realIdx, - const TMatrix &dist, - const double cutoff, - TVector &cn + const TMatrix &dist ) { double r = 0.0, rcovij = 0.0, rr = 0.0; double countf = 0.0; + double f_en = 0.0; + // initialize cn to zero + int nat = realIdx.Max() + 1; + cn.NewVector(nat); for (int i = 0, ii = 0; i != mol.NAtoms; i++) { ii = realIdx(i); @@ -335,7 +340,8 @@ int ncoord_erf( rcovij = rad[mol.ATNO(i)] + rad[mol.ATNO(j)]; rr = r / rcovij; - countf = erf_count(kn, rr); + f_en = get_en_factor(mol.ATNO(i), mol.ATNO(j)); + countf = f_en * count_fct(rr); cn(ii) += countf; cn(jj) += countf; } @@ -344,17 +350,20 @@ int ncoord_erf( return EXIT_SUCCESS; } -int dncoord_erf( + +int NCoordBase::dncoord_base( const TMolecule &mol, const TIVector &realIdx, - const TMatrix &dist, - const double cutoff, - TVector &cn, - TMatrix &dcndr + const TMatrix &dist ) { double r = 0.0, rcovij = 0.0, rr = 0.0; double rx = 0.0, ry = 0.0, rz = 0.0; double countf = 0.0, dcountf = 0.0; + double f_en = 0.0; + // initialize cn and dcndr to zero + int nat = realIdx.Max() + 1; + cn.NewVector(nat); + dcndr.NewMatrix(nat, 3 * nat); for (int i = 0, ii = 0; i != mol.NAtoms; i++) { ii = realIdx(i); @@ -367,33 +376,34 @@ int dncoord_erf( r = dist(ii, jj); if (r > cutoff) continue; - rx = (mol.CC(j, 0) - mol.CC(i, 0)) / r; - ry = (mol.CC(j, 1) - mol.CC(i, 1)) / r; - rz = (mol.CC(j, 2) - mol.CC(i, 2)) / r; + rx = (mol.CC(i, 0) - mol.CC(j, 0)) / r; + ry = (mol.CC(i, 1) - mol.CC(j, 1)) / r; + rz = (mol.CC(i, 2) - mol.CC(j, 2)) / r; rcovij = rad[mol.ATNO(i)] + rad[mol.ATNO(j)]; rr = r / rcovij; - countf = erf_count(kn, rr); + f_en = get_en_factor(mol.ATNO(i), mol.ATNO(j)); + countf = f_en * count_fct(rr); cn(ii) += countf; cn(jj) += countf; - dcountf = derf_count(kn, rr) / rcovij; - dcndr(jj, 3 * jj) += dcountf * rx; - dcndr(jj, 3 * jj + 1) += dcountf * ry; - dcndr(jj, 3 * jj + 2) += dcountf * rz; + dcountf = f_en * dcount_fct(rr) / rcovij; + dcndr(jj, 3 * jj ) -= dcountf * rx; + dcndr(jj, 3 * jj + 1) -= dcountf * ry; + dcndr(jj, 3 * jj + 2) -= dcountf * rz; - dcndr(jj, 3 * ii) += dcountf * rx; + dcndr(jj, 3 * ii ) += dcountf * rx; dcndr(jj, 3 * ii + 1) += dcountf * ry; dcndr(jj, 3 * ii + 2) += dcountf * rz; - dcndr(ii, 3 * jj) -= dcountf * rx; + dcndr(ii, 3 * jj ) -= dcountf * rx; dcndr(ii, 3 * jj + 1) -= dcountf * ry; dcndr(ii, 3 * jj + 2) -= dcountf * rz; - dcndr(ii, 3 * ii) -= dcountf * rx; - dcndr(ii, 3 * ii + 1) -= dcountf * ry; - dcndr(ii, 3 * ii + 2) -= dcountf * rz; + dcndr(ii, 3 * ii ) += dcountf * rx; + dcndr(ii, 3 * ii + 1) += dcountf * ry; + dcndr(ii, 3 * ii + 2) += dcountf * rz; } } @@ -403,150 +413,48 @@ int dncoord_erf( /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// -int get_ncoord_d4( - const TMolecule &mol, - const TMatrix &dist, - const double cutoff, - TVector &cn, - TMatrix &dcndr, - bool lgrad -) { - TIVector realIdx; - initializeRealIdx(mol.NAtoms, realIdx); - - return get_ncoord_d4(mol, realIdx, dist, cutoff, cn, dcndr, lgrad); -}; - -int get_ncoord_d4( - const TMolecule &mol, - const TIVector &realIdx, - const TMatrix &dist, - const double cutoff, - TVector &cn, - TMatrix &dcndr, - bool lgrad -) { - if (lgrad) return dncoord_d4(mol, realIdx, dist, cutoff, cn, dcndr); - return ncoord_d4(mol, realIdx, dist, cutoff, cn); -}; - -int ncoord_d4( - const TMolecule &mol, - const TIVector &realIdx, - const TMatrix &dist, - const double cutoff, - TVector &cn -) { - double r = 0.0, rcovij = 0.0, rr = 0.0; - double den = 0.0; - double countf = 0.0; - int izp, jzp; - - for (int i = 0, ii = 0; i != mol.NAtoms; i++) { - ii = realIdx(i); - if (ii < 0) continue; - - izp = mol.ATNO(i); - for (int j = 0, jj = 0; j != i; j++) { - jj = realIdx(j); - if (jj < 0) continue; - - r = dist(ii, jj); - if (r > cutoff) continue; - - jzp = mol.ATNO(j); - rcovij = rad[izp] + rad[jzp]; - rr = r / rcovij; - den = k4 * exp(-pow((fabs(en[izp] - en[jzp]) + k5), 2) / k6); - countf = den * erf_count(kn, rr); - - cn(ii) += countf; - cn(jj) += countf; - } - } - return EXIT_SUCCESS; +double NCoordBase::get_en_factor(int i, int j) const { + return 1.0; } -int dncoord_d4( - const TMolecule &mol, - const TIVector &realIdx, - const TMatrix &dist, - const double cutoff, - TVector &cn, - TMatrix &dcndr -) { - double r = 0.0, rcovij = 0.0, rr = 0.0; - double rx = 0.0, ry = 0.0, rz = 0.0; - double countf = 0.0, dcountf = 0.0, den = 0.0; - int izp, jzp; - - for (int i = 0, ii = 0; i != mol.NAtoms; i++) { - ii = realIdx(i); - if (ii < 0) continue; - - izp = mol.ATNO(i); - for (int j = 0, jj = 0; j != i; j++) { - jj = realIdx(j); - if (jj < 0) continue; - - r = dist(ii, jj); - if (r > cutoff) continue; - - jzp = mol.ATNO(j); - rx = (mol.CC(j, 0) - mol.CC(i, 0)) / r; - ry = (mol.CC(j, 1) - mol.CC(i, 1)) / r; - rz = (mol.CC(j, 2) - mol.CC(i, 2)) / r; - - rcovij = rad[izp] + rad[jzp]; - rr = r / rcovij; - den = k4 * exp(-pow((fabs(en[izp] - en[jzp]) + k5), 2) / k6); - countf = den * erf_count(kn, rr); - cn(ii) += countf; - cn(jj) += countf; - - dcountf = den * derf_count(kn, rr) / rcovij; - dcndr(3 * jj, jj) += dcountf * rx; - dcndr(3 * jj + 1, jj) += dcountf * ry; - dcndr(3 * jj + 2, jj) += dcountf * rz; - dcndr(3 * jj, ii) = dcountf * rx; - dcndr(3 * jj + 1, ii) = dcountf * ry; - dcndr(3 * jj + 2, ii) = dcountf * rz; - dcndr(3 * ii, jj) = -dcountf * rx; - dcndr(3 * ii + 1, jj) = -dcountf * ry; - dcndr(3 * ii + 2, jj) = -dcountf * rz; - dcndr(3 * ii, ii) += -dcountf * rx; - dcndr(3 * ii + 1, ii) += -dcountf * ry; - dcndr(3 * ii + 2, ii) += -dcountf * rz; - } - } - return EXIT_SUCCESS; +double NCoordErfD4::get_en_factor(int i, int j) const { + return k4 * exp(-pow((fabs(en[i] - en[j]) + k5), 2) / k6); } /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// -double erf_count(double k, double rr) { - return 0.5 * (1.0 + erf(-k * (rr - 1.0))); +double NCoordErf::count_fct(double rr) const { + return 0.5 * (1.0 + erf(-kcn * (rr - 1.0))); } -double derf_count(double k, double rr) { - return -k * hlfosqrtpi * exp(-pow(k * (rr - 1.0), 2)); +double NCoordErf::dcount_fct(double rr) const { + return -kcn * hlfosqrtpi * exp(-pow(kcn * (rr - 1.0), 2)); } -int cut_coordination_number( +double NCoordErfD4::count_fct(double rr) const { + return 0.5 * (1.0 + erf(-kcn * (rr - 1.0))); +} + +double NCoordErfD4::dcount_fct(double rr) const { + return -kcn * hlfosqrtpi * exp(-pow(kcn * (rr - 1.0), 2)); +} + +int NCoordErf::cut_coordination_number( const double cn_max, TVector &cn, TMatrix &dcndr, bool lgrad ) { if (lgrad) { + // cutting the cn is not (anti)symmetric, so dcndr is not antisymmetric anymore double dcnpdcn; for (int i = 0; i != cn.N; i++) { dcnpdcn = dlog_cn_cut(cn_max, cn(i)); for (int j = 0; j != cn.N; j++) { - dcndr(j, 3 * i) *= dcnpdcn; - dcndr(j, 3 * i + 1) *= dcnpdcn; - dcndr(j, 3 * i + 2) *= dcnpdcn; + dcndr(i, 3 * j ) *= dcnpdcn; + dcndr(i, 3 * j + 1) *= dcnpdcn; + dcndr(i, 3 * j + 2) *= dcnpdcn; } } } @@ -566,4 +474,13 @@ inline double dlog_cn_cut(const double cn_max, const double cn) { return exp(cn_max) / (exp(cn_max) + exp(cn)); } +int NCoordErfD4::cut_coordination_number( + const double cn_max, + TVector &cn, + TMatrix &dcndr, + bool lgrad +) { + return EXIT_SUCCESS; +} + } // namespace dftd4 diff --git a/test/test_ghost.cpp b/test/test_ghost.cpp index 6381122..401010c 100644 --- a/test/test_ghost.cpp +++ b/test/test_ghost.cpp @@ -59,16 +59,14 @@ int test_water( // COORDINATION NUMBER CHECK // erf-CN without cutoff - TVector cn; - TMatrix dcndr; // empty because no gradient needed - cn.New(nat); - info = get_ncoord_d4(mol, realIdx, dist, 9999.9, cn, dcndr, false); + NCoordErfD4 ncoord_erf(7.5, 1.0, 9999.9); + info = ncoord_erf.get_ncoord(mol, realIdx, dist, false); if (info != EXIT_SUCCESS) return info; // compare to ref for (int i = 0; i != nat; i++) { - if (check(cn(i), water_dimer_ref_cn[i]) == EXIT_FAILURE) { - print_fail("GHOST: CN_D4", cn(i), water_dimer_ref_cn[i]); + if (check(ncoord_erf.cn(i), water_dimer_ref_cn[i]) == EXIT_FAILURE) { + print_fail("GHOST: CN_D4", ncoord_erf.cn(i), water_dimer_ref_cn[i]); return EXIT_FAILURE; } } @@ -148,8 +146,7 @@ int test_water( d4grad[i] = 0.0; } - dcndr.NewMatrix(3 * nat, nat); - info = get_ncoord_d4(mol, realIdx, dist, cutoff.cn, cn, dcndr, lgrad); + info = ncoord_erf.get_ncoord(mol, realIdx, dist, lgrad); if (info != EXIT_SUCCESS) return info; dqdr.NewMatrix(3 * nat, nat); @@ -162,6 +159,7 @@ int test_water( for (int i = 0; i < 3 * mol.NAtoms; i++) { if (check(d4grad[i], ref_grad[i], 1e-8) == EXIT_FAILURE) { print_fail("GHOST: Gradient", d4grad[i], ref_grad[i]); + delete[] d4grad; return EXIT_FAILURE; } } diff --git a/test/test_grad.cpp b/test/test_grad.cpp index ca82b97..d2ecdd5 100644 --- a/test/test_grad.cpp +++ b/test/test_grad.cpp @@ -20,6 +20,8 @@ #include #include #include +#include +#include #include "molecules.h" #include "test_grad.h" @@ -81,6 +83,7 @@ int test_numgrad(TMolecule &mol, const int charge, const dparam &par) { for (int c = 0; c < 3; c++) { if (check(d4grad[3 * i + c], numgrad(i, c), thr) != EXIT_SUCCESS) { print_fail("Gradient mismatch", d4grad[3 * i + c], numgrad(i, c)); + delete[] d4grad; return EXIT_FAILURE; } } @@ -126,6 +129,7 @@ int test_pbed4_mb01() { par.s8 = 0.95948085; par.a1 = 0.38574991; par.a2 = 4.80688534; + par.s10 = 0.0; // assemble molecule int charge = mb16_43_01_charge; @@ -165,9 +169,98 @@ int test_tpss0d4mbd_rost61m1() { return test_numgrad(mol, charge, par); } + +int test_numgrad_dqdr( + int n, + const char atoms[][3], + const double coord[] +) { + // assemble molecule + int info; + TMolecule mol; + info = get_molecule(n, atoms, coord, mol); + if (info != EXIT_SUCCESS) return info; + TVector q_r, q_l; + double step{1.0e-6}; // accurate up to 1.0E-8 + double thr{1.2e-8}; + + TMatrix dist; + dist.NewMat(mol.NAtoms, mol.NAtoms); + TMatrix num_dqdr; // numerical gradient of the partial charges + num_dqdr.NewMat(3 * mol.NAtoms, mol.NAtoms); + TMatrix analytic_dqdr; // analytical gradient of the partial charges + analytic_dqdr.NewMat(3 * mol.NAtoms, mol.NAtoms); + multicharge::EEQModel eeq_model; + TVector q; + TMatrix dqdr; + q.NewVec(mol.NAtoms); + dqdr.NewMat(3 * mol.NAtoms, mol.NAtoms); + + + TCutoff cutoff; + + // masking (nothing excluded) + TVector realIdx; + realIdx.NewVec(mol.NAtoms); + int nat = 0; + for (int i = 0; i != mol.NAtoms; i++) { + realIdx(i) = nat; + nat++; + } + + // analytical gradient + calc_distances(mol, realIdx, dist); + info = + eeq_model.get_charges(mol, realIdx, dist, 0, cutoff.cn_eeq, q, analytic_dqdr, true); + if (info != EXIT_SUCCESS) return info; + + // calculate numerical gradient via finite difference method + for (int i = 0; i < mol.NAtoms; i++) { + for (int c = 0; c < 3; c++) { + // calculate forward point + q_r.NewVec(n); + mol.CC(i, c) += step; + calc_distances(mol, realIdx, dist); + eeq_model.get_charges(mol, realIdx, dist, 0, cutoff.cn_eeq, q_r, dqdr, false); + + // calculate backward point + q_l.NewVec(n); + mol.CC(i, c) = mol.CC(i, c) - 2 * step; + calc_distances(mol, realIdx, dist); + eeq_model.get_charges(mol, realIdx, dist, 0, cutoff.cn_eeq, q_l, dqdr, false); + + // calculate numerical gradient as finite difference + mol.CC(i, c) = mol.CC(i, c) + step; + for (int j = 0; j < mol.NAtoms; j++) { + num_dqdr(3 * i + c, j) = 0.5 * (q_r(j) - q_l(j)) / step; + + } + } + } + + // compare against numerical gradient + for (int i = 0; i < mol.NAtoms; i++) { + for (int c = 0; c < 3; c++) { + for (int j = 0; j < mol.NAtoms; j++) { + if (check(analytic_dqdr(3 * i + c, j), num_dqdr(3 * i + c, j), thr) != EXIT_SUCCESS) { + print_fail("Gradient mismatch for dqdr\n", analytic_dqdr(3 * i + c, j), num_dqdr(3 * i + c, j)); + return EXIT_FAILURE; + } + } + } + } + + return EXIT_SUCCESS; +} + + int test_grad() { int info{0}; + info = test_numgrad_dqdr( + mb16_43_01_n, mb16_43_01_atoms, mb16_43_01_coord); + if (info != EXIT_SUCCESS) return info; + info = test_pbed4_mb01(); if (info != EXIT_SUCCESS) return info; diff --git a/test/test_ncoord.cpp b/test/test_ncoord.cpp index f9c1612..21cd503 100644 --- a/test/test_ncoord.cpp +++ b/test/test_ncoord.cpp @@ -61,16 +61,14 @@ int test_cn( calc_distances(mol, realIdx, dist); // erf-CN without cutoff - TVector cn; - TMatrix dcndr; // empty because no gradient needed - cn.New(n); - info = get_ncoord_d4(mol, realIdx, dist, 9999.9, cn, dcndr, false); + NCoordErfD4 ncoord_erf(7.5, 1.0, 9999.9); + info = ncoord_erf.get_ncoord(mol, realIdx, dist, false); if (info != EXIT_SUCCESS) return info; // compare to ref for (int i = 0; i != n; i++) { - if (check(cn(i), ref(i)) == EXIT_FAILURE) { - print_fail("CN_D4", cn(i), ref(i)); + if (check(ncoord_erf.cn(i), ref(i)) == EXIT_FAILURE) { + print_fail("CN_D4", ncoord_erf.cn(i), ref(i)); return EXIT_FAILURE; } } @@ -78,6 +76,167 @@ int test_cn( return EXIT_SUCCESS; } + +int test_numgrad_d4( + int n, + const char atoms[][3], + const double coord[] +) { + // assemble molecule + int info; + TMolecule mol; + info = get_molecule(n, atoms, coord, mol); + if (info != EXIT_SUCCESS) return info; + double step{1.0e-6}; + double thr{1e-8}; + + TMatrix dist; + dist.NewMat(mol.NAtoms, mol.NAtoms); + TMatrix num_dcndr; // numerical gradient of the coordination number + num_dcndr.NewMat(mol.NAtoms, 3 * mol.NAtoms); + TMatrix analytic_dcndr; // analytical gradient of the coordination number + analytic_dcndr.NewMat(mol.NAtoms, 3 * mol.NAtoms); + NCoordErfD4 ncoord_erf_d4; + NCoordErfD4 ncoord_erf_d4_r; + NCoordErfD4 ncoord_erf_d4_l; + + TCutoff cutoff; + + // masking (nothing excluded) + TVector realIdx; + realIdx.NewVec(mol.NAtoms); + int nat = 0; + for (int i = 0; i != mol.NAtoms; i++) { + realIdx(i) = nat; + nat++; + } + + // analytical gradient + calc_distances(mol, realIdx, dist); + ncoord_erf_d4.get_ncoord(mol, realIdx, dist, true); + analytic_dcndr.CopyMat(ncoord_erf_d4.dcndr); + + // check if analytical gradient is antisymmetric + for (int c = 0; c < 3; c++) { + for (int i = 0; i < mol.NAtoms; i++) { + for (int k = 0; k < i; k++) { + if (abs(analytic_dcndr(k, 3 * i + c) + analytic_dcndr(i, 3 * k + c)) > 1.0e-9 ) { + print_fail("Analytical CN-gradient is not antisymmetric for NCoordErfD4",analytic_dcndr(k, 3 * i + c) + analytic_dcndr(i, 3 * k + c) , 0.0); + return EXIT_FAILURE; + } + } + } + } + + // numerical gradient + for (int i = 0; i < mol.NAtoms; i++) { + for (int c = 0; c < 3; c++) { + mol.CC(i, c) += step; + calc_distances(mol, realIdx, dist); + ncoord_erf_d4_r.get_ncoord(mol, realIdx, dist, false); + + mol.CC(i, c) = mol.CC(i, c) - 2 * step; + calc_distances(mol, realIdx, dist); + ncoord_erf_d4_l.get_ncoord(mol, realIdx, dist, false); + + mol.CC(i, c) = mol.CC(i, c) + step; + for (int j = 0; j < mol.NAtoms; j++) { + // numerical CN gradient: dCN(j)/ dr(i)^c with c = x,y, or z + num_dcndr(j, 3 * i + c) = 0.5 * (ncoord_erf_d4_r.cn(j) - ncoord_erf_d4_l.cn(j)) / step; + } + } + } + + // compare against numerical gradient + for (int i = 0; i < mol.NAtoms; i++) { + for (int c = 0; c < 3; c++) { + for (int j = 0; j < mol.NAtoms; j++) { + if (check(analytic_dcndr(j, 3 * i + c), num_dcndr(j, 3 * i + c), thr) != EXIT_SUCCESS) { + print_fail("Gradient mismatch for NCoordErfD4 dcndr", analytic_dcndr(j, 3 * i + c), num_dcndr(j, 3 * i + c)); + return EXIT_FAILURE; + } + } + } + } + + return EXIT_SUCCESS; +} + + +int test_numgrad( + int n, + const char atoms[][3], + const double coord[] +) { + // assemble molecule + int info; + TMolecule mol; + info = get_molecule(n, atoms, coord, mol); + if (info != EXIT_SUCCESS) return info; + double step{1.0e-6}; + double thr{1e-8}; + + TMatrix dist; + dist.NewMat(mol.NAtoms, mol.NAtoms); + TMatrix num_dcndr; // numerical gradient of the coordination number + num_dcndr.NewMat(mol.NAtoms, 3 * mol.NAtoms); + TMatrix analytic_dcndr; // analytical gradient of the coordination number + analytic_dcndr.NewMat(mol.NAtoms, 3 * mol.NAtoms); + NCoordErf ncoord_erf; + NCoordErf ncoord_erf_r; + NCoordErf ncoord_erf_l; + + TCutoff cutoff; + + // masking (nothing excluded) + TVector realIdx; + realIdx.NewVec(mol.NAtoms); + int nat = 0; + for (int i = 0; i != mol.NAtoms; i++) { + realIdx(i) = nat; + nat++; + } + + // analytical gradient + calc_distances(mol, realIdx, dist); + ncoord_erf.get_ncoord(mol, realIdx, dist, true); + analytic_dcndr.CopyMat(ncoord_erf.dcndr); + + // numerical gradient + for (int i = 0; i < mol.NAtoms; i++) { + for (int c = 0; c < 3; c++) { + mol.CC(i, c) += step; + calc_distances(mol, realIdx, dist); + ncoord_erf_r.get_ncoord(mol, realIdx, dist, false); + + mol.CC(i, c) = mol.CC(i, c) - 2 * step; + calc_distances(mol, realIdx, dist); + ncoord_erf_l.get_ncoord(mol, realIdx, dist, false); + + mol.CC(i, c) = mol.CC(i, c) + step; + for (int j = 0; j < mol.NAtoms; j++) { + // numerical CN gradient: dCN(j)/ dr(i)^c with c = x,y, or z + num_dcndr(j, 3 * i + c) = 0.5 * (ncoord_erf_r.cn(j) - ncoord_erf_l.cn(j)) / step; + } + } + } + + // compare against numerical gradient + for (int i = 0; i < mol.NAtoms; i++) { + for (int c = 0; c < 3; c++) { + for (int j = 0; j < mol.NAtoms; j++) { + if (check(analytic_dcndr(j, 3 * i + c), num_dcndr(j, 3 * i + c), thr) != EXIT_SUCCESS) { + print_fail("Gradient mismatch for NCoordErf dcndr", analytic_dcndr(j, 3 * i + c), num_dcndr(j, 3 * i + c)); + return EXIT_FAILURE; + } + } + } + } + + return EXIT_SUCCESS; +} + + int test_ncoord() { int info; @@ -93,5 +252,24 @@ int test_ncoord() { test_cn(rost61_m1_n, rost61_m1_atoms, rost61_m1_coord, rost61_m1_ref_cn); if (info != EXIT_SUCCESS) return info; + info = test_numgrad_d4( + mb16_43_01_n, mb16_43_01_atoms, mb16_43_01_coord + ); + if (info != EXIT_SUCCESS) return info; + + info = test_numgrad( + mb16_43_01_n, mb16_43_01_atoms, mb16_43_01_coord + ); + if (info != EXIT_SUCCESS) return info; + + info = test_numgrad_d4( + water_n, water_atoms, water_coord + ); + if (info != EXIT_SUCCESS) return info; + + info = test_numgrad( + water_n, water_atoms, water_coord + ); + if (info != EXIT_SUCCESS) return info; return EXIT_SUCCESS; }