From 5f035a2d4bb9cb69dffd41d8b87d1aa8dbea1a2b Mon Sep 17 00:00:00 2001 From: Thomas Rose Date: Fri, 6 Jun 2025 17:56:54 +0200 Subject: [PATCH 01/12] refactor: move charge models to new multicharge namespace and create ChargeModel class --- include/qceeq_param.h | 110 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 110 insertions(+) create mode 100644 include/qceeq_param.h diff --git a/include/qceeq_param.h b/include/qceeq_param.h new file mode 100644 index 0000000..f6e7260 --- /dev/null +++ b/include/qceeq_param.h @@ -0,0 +1,110 @@ +/* +Electronegativity equilibration (EEQ_BC) model + +This file contains the model parameters and getter functions + +Thomas Rose () + +*/ + +#ifndef QCEEQ_PARAM_H +#define QCEEQ_PARAM_H + +namespace EEQ_param { + +constexpr int MAXELEMENT = 104; // 103 + dummy + +// Element-specific electronegativity for the EEQ charges. +static const double xi[MAXELEMENT]{ + +0.00000000, // dummy + +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, +1.06598299, +0.68184178, +1.04581665, +1.09888688, + +1.07206461, +1.09821942, +1.10900303, +1.01039812, +1.00095966, +1.11003303, + +1.16831853, +1.00887482, +1.05928842, +1.07672363, +1.11308426, +1.14340090, + +1.13714110, +}; + +// Element-specific chemical hardnesses for the EEQ charges. +static const double gam[MAXELEMENT]{ + +0.00000000, // dummy + -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.10034731, +0.15801991, -0.00071039, -0.00170887, + -0.00133327, -0.00104386, -0.00094936, -0.00111390, -0.00125257, -0.00095936, + -0.00102814, -0.00104450, -0.00112666, -0.00101529, -0.00059592, -0.00012585, + -0.00140896, +}; + +// Element-specific CN scaling constant for the EEQ charges. +static const double kappa[MAXELEMENT]{ + +0.00000000, // dummy + +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.00013506, -0.00020631, +0.00473118, +0.01590519, + +0.00369763, +0.00417543, +0.00706682, +0.00488679, +0.00505103, +0.00710682, + +0.00463050, +0.00387799, +0.00296795, +0.00400648, +0.00548481, +0.01350400, + +0.00675380, +}; + +static const double alp[MAXELEMENT]{ + +0.00000000, // dummy + +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, +1.04059593, +0.60550051, +1.22262145, +1.28736399, + +1.44431317, +1.29032833, +1.41009404, +1.25501213, +1.15181468, +1.42010424, + +1.43955530, +1.28565237, +1.35017463, +1.33011749, +1.30745135, +1.26526071, + +1.34071499, +}; + +} + +#endif // QCEEQ_PARAM_H \ No newline at end of file From 7fec355261d41efd6115474f1ac3ac9b723cb843 Mon Sep 17 00:00:00 2001 From: Thomas Rose Date: Thu, 12 Jun 2025 13:18:31 +0200 Subject: [PATCH 02/12] refactor: create class for coordination number calculation --- include/dftd_ncoord.h | 73 +++++++++++++++++++++++++++++++++ src/dftd_dispersion.cpp | 21 +++++----- src/dftd_eeq.cpp | 15 ++++--- src/dftd_ncoord.cpp | 91 +++++++++++++++++++---------------------- test/test_ghost.cpp | 15 ++++--- test/test_ncoord.cpp | 11 +++-- 6 files changed, 145 insertions(+), 81 deletions(-) diff --git a/include/dftd_ncoord.h b/include/dftd_ncoord.h index 326382c..d9f7c27 100644 --- a/include/dftd_ncoord.h +++ b/include/dftd_ncoord.h @@ -28,6 +28,79 @@ namespace dftd4 { + +class NCoordBase +{ + public: + TVector cn; + TMatrix dcndr; + static const double rad[]; + double kcn; // Steepness of counting function + double norm_exp; + double cutoff; + // Get the coordination number + int get_ncoord( // with ghost atoms + const TMolecule&, + const TMatrix&, + const double, + TVector&, + TMatrix&, + bool); + int get_ncoord( // without ghost atoms + const TMolecule&, + const TIVector&, + const TMatrix&, + bool); + // Calculate the coordination number using the virtual counting function + int ncoord_base( + const TMolecule&, + const TIVector&, + const TMatrix&); + // Calculate the derivative of the coordination number + int dr_ncoord_base( + const TMolecule&, + const TIVector&, + const TMatrix&); + // Get the DFT-D4 coordination number + int get_ncoord_d4( + const TMolecule&, + const TMatrix&, + bool); + int get_ncoord_d4( + const TMolecule&, + const TIVector&, + const TMatrix&, + bool); + int ncoord_d4( + const TMolecule&, + const TIVector&, + const TMatrix&); + int dncoord_d4( + const TMolecule&, + const TIVector&, + const TMatrix&); + + // Counting function + virtual double count_fct(double) const = 0; + // Derivative of the counting function + virtual double dr_count_fct(double) const = 0; + // Constructor + NCoordBase(double, double, double); + // Virtual destructor + virtual ~NCoordBase() = default; +}; + +class NCoordErf : public NCoordBase { + public: + // erf() based counting function + double count_fct(double) const override; + // derivative of the erf() based counting function + double dr_count_fct(double) const 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){} +}; + /** * Calculate all distance pairs and store in matrix. * diff --git a/src/dftd_dispersion.cpp b/src/dftd_dispersion.cpp index bfbea23..eadfb84 100644 --- a/src/dftd_dispersion.cpp +++ b/src/dftd_dispersion.cpp @@ -102,17 +102,16 @@ 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 + NCoordErf ncoord_erf; multicharge::EEQModel chrg_model; // Charge model - cn.NewVector(nat); + ncoord_erf.cn.NewVector(nat); q.NewVector(nat); if (lgrad) { - dcndr.NewMatrix(3 * nat, nat); + ncoord_erf.dcndr.NewMatrix(3 * nat, nat); dqdr.NewMatrix(3 * nat, nat); gradient.NewVector(3 * nat); } @@ -122,7 +121,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_d4(mol, realIdx, dist, lgrad); if (info != EXIT_SUCCESS) return info; // maximum number of reference systems @@ -145,7 +144,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 +211,11 @@ 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(); + ncoord_erf.cn.Delete(); q.Delete(); refq.Delete(); @@ -253,7 +252,7 @@ int get_dispersion( ); if (info != EXIT_SUCCESS) return info; } else { - cn.Delete(); + ncoord_erf.cn.Delete(); q.Delete(); refq.Delete(); gwvec.Delete(); @@ -266,9 +265,9 @@ 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, false, 1.0); } - dcndr.DelMat(); + ncoord_erf.dcndr.DelMat(); dEdcn.DelVec(); dEdq.DelVec(); diff --git a/src/dftd_eeq.cpp b/src/dftd_eeq.cpp index 2f6c435..6c43e00 100644 --- a/src/dftd_eeq.cpp +++ b/src/dftd_eeq.cpp @@ -72,23 +72,22 @@ int ChargeModel::get_charges( bool lverbose{false}; int nat = realIdx.Max() + 1; - TVector cn; // EEQ cordination number - TMatrix dcndr; // Derivative of EEQ-CN + dftd4::NCoordErf ncoord_erf; - cn.NewVec(nat); - if (lgrad) dcndr.NewMat(nat, 3 * nat); + ncoord_erf.cn.NewVec(nat); + if (lgrad) ncoord_erf.dcndr.NewMat(nat, 3 * nat); // 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(); + ncoord_erf.dcndr.DelMat(); + ncoord_erf.cn.DelVec(); return EXIT_SUCCESS; }; diff --git a/src/dftd_ncoord.cpp b/src/dftd_ncoord.cpp index a6cf431..ec2207b 100644 --- a/src/dftd_ncoord.cpp +++ b/src/dftd_ncoord.cpp @@ -32,6 +32,17 @@ namespace dftd4 { +NCoordBase::NCoordBase( + double optional_kcn = 7.5, // defaults for D4 EEQ + double optional_norm_exp = 1.0, // defaults for D4 EEQ + double optional_cutoff = 25.0 //@TR NOTE check if this is EEQ default +){ + NCoordBase::kcn = optional_kcn; + NCoordBase::norm_exp = optional_norm_exp; + NCoordBase::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 +80,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,7 +285,7 @@ void initializeRealIdx(int nat, TVector &realIdx) { /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// -int get_ncoord_erf( +int NCoordBase::get_ncoord( const TMolecule &mol, const TMatrix &dist, const double cutoff, @@ -285,24 +296,21 @@ int get_ncoord_erf( 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 = dr_ncoord_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,12 +320,10 @@ 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; @@ -335,7 +341,7 @@ int ncoord_erf( rcovij = rad[mol.ATNO(i)] + rad[mol.ATNO(j)]; rr = r / rcovij; - countf = erf_count(kn, rr); + countf = count_fct(rr); cn(ii) += countf; cn(jj) += countf; } @@ -344,13 +350,10 @@ int ncoord_erf( return EXIT_SUCCESS; } -int dncoord_erf( +int NCoordBase::dr_ncoord_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; @@ -374,11 +377,11 @@ int dncoord_erf( rcovij = rad[mol.ATNO(i)] + rad[mol.ATNO(j)]; rr = r / rcovij; - countf = erf_count(kn, rr); + countf = count_fct(rr); cn(ii) += countf; cn(jj) += countf; - dcountf = derf_count(kn, rr) / rcovij; + dcountf = dr_count_fct(rr) / rcovij; dcndr(jj, 3 * jj) += dcountf * rx; dcndr(jj, 3 * jj + 1) += dcountf * ry; dcndr(jj, 3 * jj + 2) += dcountf * rz; @@ -403,39 +406,31 @@ int dncoord_erf( /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// -int get_ncoord_d4( +int NCoordBase::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); + return get_ncoord_d4(mol, realIdx, dist, lgrad); }; -int get_ncoord_d4( +int NCoordBase::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); + if (lgrad) return dncoord_d4(mol, realIdx, dist); + return ncoord_d4(mol, realIdx, dist); }; -int ncoord_d4( +int NCoordBase::ncoord_d4( 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 den = 0.0; @@ -458,7 +453,7 @@ int ncoord_d4( 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); + countf = den * count_fct(rr); cn(ii) += countf; cn(jj) += countf; @@ -467,13 +462,10 @@ int ncoord_d4( return EXIT_SUCCESS; } -int dncoord_d4( +int NCoordBase::dncoord_d4( 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; @@ -500,11 +492,11 @@ int dncoord_d4( 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); + countf = den * count_fct(rr); cn(ii) += countf; cn(jj) += countf; - dcountf = den * derf_count(kn, rr) / rcovij; + dcountf = den * dr_count_fct(rr) / rcovij; dcndr(3 * jj, jj) += dcountf * rx; dcndr(3 * jj + 1, jj) += dcountf * ry; dcndr(3 * jj + 2, jj) += dcountf * rz; @@ -525,12 +517,15 @@ int dncoord_d4( /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// -double erf_count(double k, double rr) { - return 0.5 * (1.0 + erf(-k * (rr - 1.0))); +//NCoordErf::NCoordErf(double optional_kcn, double optional_norm_exp, double optional_cutoff) +// : NCoordBase(optional_kcn, optional_norm_exp, optional_cutoff) {} + +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::dr_count_fct(double rr) const { + return -kcn * hlfosqrtpi * exp(-pow(kcn * (rr - 1.0), 2)); } int cut_coordination_number( diff --git a/test/test_ghost.cpp b/test/test_ghost.cpp index 6381122..f613823 100644 --- a/test/test_ghost.cpp +++ b/test/test_ghost.cpp @@ -59,16 +59,15 @@ 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); + NCoordErf ncoord_erf(7.5, 1.0, 9999.9); + ncoord_erf.cn.New(nat); + info = ncoord_erf.get_ncoord_d4(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 +147,8 @@ int test_water( d4grad[i] = 0.0; } - dcndr.NewMatrix(3 * nat, nat); - info = get_ncoord_d4(mol, realIdx, dist, cutoff.cn, cn, dcndr, lgrad); + ncoord_erf.dcndr.NewMatrix(3 * nat, nat); + info = ncoord_erf.get_ncoord_d4(mol, realIdx, dist, lgrad); if (info != EXIT_SUCCESS) return info; dqdr.NewMatrix(3 * nat, nat); diff --git a/test/test_ncoord.cpp b/test/test_ncoord.cpp index f9c1612..7b77da5 100644 --- a/test/test_ncoord.cpp +++ b/test/test_ncoord.cpp @@ -61,16 +61,15 @@ 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); + NCoordErf ncoord_erf(7.5, 1.0, 9999.9); + ncoord_erf.cn.New(n); + info = ncoord_erf.get_ncoord_d4(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; } } From b85b2c74f0bacbbdea3c8464e11ddea758ef3a30 Mon Sep 17 00:00:00 2001 From: Thomas Rose Date: Fri, 13 Jun 2025 11:18:27 +0200 Subject: [PATCH 03/12] Update comment in dftd_ncoord.cpp --- src/dftd_ncoord.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dftd_ncoord.cpp b/src/dftd_ncoord.cpp index ec2207b..537d08f 100644 --- a/src/dftd_ncoord.cpp +++ b/src/dftd_ncoord.cpp @@ -35,7 +35,7 @@ namespace dftd4 { NCoordBase::NCoordBase( double optional_kcn = 7.5, // defaults for D4 EEQ double optional_norm_exp = 1.0, // defaults for D4 EEQ - double optional_cutoff = 25.0 //@TR NOTE check if this is EEQ default + double optional_cutoff = 25.0 // defaults for D4 EEQ ){ NCoordBase::kcn = optional_kcn; NCoordBase::norm_exp = optional_norm_exp; From 6886aa7d818804c38b7964a93c703f890741da0f Mon Sep 17 00:00:00 2001 From: Thomas Rose Date: Wed, 18 Jun 2025 10:42:19 +0200 Subject: [PATCH 04/12] Move cn and dcndr allocation into ncoord and dncoord functions. --- src/dftd_dispersion.cpp | 2 -- src/dftd_eeq.cpp | 3 --- src/dftd_ncoord.cpp | 14 ++++++++++++++ test/test_ghost.cpp | 2 -- test/test_grad.cpp | 1 + test/test_ncoord.cpp | 1 - 6 files changed, 15 insertions(+), 8 deletions(-) diff --git a/src/dftd_dispersion.cpp b/src/dftd_dispersion.cpp index eadfb84..717fb61 100644 --- a/src/dftd_dispersion.cpp +++ b/src/dftd_dispersion.cpp @@ -108,10 +108,8 @@ int get_dispersion( NCoordErf ncoord_erf; multicharge::EEQModel chrg_model; // Charge model - ncoord_erf.cn.NewVector(nat); q.NewVector(nat); if (lgrad) { - ncoord_erf.dcndr.NewMatrix(3 * nat, nat); dqdr.NewMatrix(3 * nat, nat); gradient.NewVector(3 * nat); } diff --git a/src/dftd_eeq.cpp b/src/dftd_eeq.cpp index 6c43e00..4c5e70e 100644 --- a/src/dftd_eeq.cpp +++ b/src/dftd_eeq.cpp @@ -74,9 +74,6 @@ int ChargeModel::get_charges( dftd4::NCoordErf ncoord_erf; - ncoord_erf.cn.NewVec(nat); - if (lgrad) ncoord_erf.dcndr.NewMat(nat, 3 * nat); - // get the EEQ coordination number info = ncoord_erf.get_ncoord(mol, realIdx, dist, lgrad); if (info != EXIT_SUCCESS) return info; diff --git a/src/dftd_ncoord.cpp b/src/dftd_ncoord.cpp index 537d08f..37ca84c 100644 --- a/src/dftd_ncoord.cpp +++ b/src/dftd_ncoord.cpp @@ -327,6 +327,9 @@ int NCoordBase::ncoord_base( ) { double r = 0.0, rcovij = 0.0, rr = 0.0; double countf = 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); @@ -358,6 +361,10 @@ int NCoordBase::dr_ncoord_base( 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; + // 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); @@ -436,6 +443,9 @@ int NCoordBase::ncoord_d4( double den = 0.0; double countf = 0.0; int izp, jzp; + // 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); @@ -471,6 +481,10 @@ int NCoordBase::dncoord_d4( double rx = 0.0, ry = 0.0, rz = 0.0; double countf = 0.0, dcountf = 0.0, den = 0.0; int izp, jzp; + // initialize cn and dcndr to zero + int nat = realIdx.Max() + 1; + cn.NewVector(nat); + dcndr.NewMatrix(3*nat, nat); for (int i = 0, ii = 0; i != mol.NAtoms; i++) { ii = realIdx(i); diff --git a/test/test_ghost.cpp b/test/test_ghost.cpp index f613823..f402891 100644 --- a/test/test_ghost.cpp +++ b/test/test_ghost.cpp @@ -60,7 +60,6 @@ int test_water( // erf-CN without cutoff NCoordErf ncoord_erf(7.5, 1.0, 9999.9); - ncoord_erf.cn.New(nat); info = ncoord_erf.get_ncoord_d4(mol, realIdx, dist, false); if (info != EXIT_SUCCESS) return info; @@ -147,7 +146,6 @@ int test_water( d4grad[i] = 0.0; } - ncoord_erf.dcndr.NewMatrix(3 * nat, nat); info = ncoord_erf.get_ncoord_d4(mol, realIdx, dist, lgrad); if (info != EXIT_SUCCESS) return info; diff --git a/test/test_grad.cpp b/test/test_grad.cpp index ca82b97..21ef6ce 100644 --- a/test/test_grad.cpp +++ b/test/test_grad.cpp @@ -126,6 +126,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; diff --git a/test/test_ncoord.cpp b/test/test_ncoord.cpp index 7b77da5..7a855f7 100644 --- a/test/test_ncoord.cpp +++ b/test/test_ncoord.cpp @@ -62,7 +62,6 @@ int test_cn( // erf-CN without cutoff NCoordErf ncoord_erf(7.5, 1.0, 9999.9); - ncoord_erf.cn.New(n); info = ncoord_erf.get_ncoord_d4(mol, realIdx, dist, false); if (info != EXIT_SUCCESS) return info; From c65dc277047064144f8cdfc423eaeb2bf4cc15a2 Mon Sep 17 00:00:00 2001 From: Thomas Rose Date: Wed, 18 Jun 2025 14:54:08 +0200 Subject: [PATCH 05/12] NCoordBase constructor initializer list and code cleanup --- include/dftd_ncoord.h | 2 +- include/qceeq_param.h | 110 ------------------------------------------ src/dftd_ncoord.cpp | 16 ++---- 3 files changed, 6 insertions(+), 122 deletions(-) delete mode 100644 include/qceeq_param.h diff --git a/include/dftd_ncoord.h b/include/dftd_ncoord.h index d9f7c27..4ef9e6d 100644 --- a/include/dftd_ncoord.h +++ b/include/dftd_ncoord.h @@ -85,7 +85,7 @@ class NCoordBase // Derivative of the counting function virtual double dr_count_fct(double) const = 0; // Constructor - NCoordBase(double, double, double); + NCoordBase(double optional_kcn = 7.5, double optional_norm_exp = 1.0, double optional_cutoff = 25.0); // Virtual destructor virtual ~NCoordBase() = default; }; diff --git a/include/qceeq_param.h b/include/qceeq_param.h deleted file mode 100644 index f6e7260..0000000 --- a/include/qceeq_param.h +++ /dev/null @@ -1,110 +0,0 @@ -/* -Electronegativity equilibration (EEQ_BC) model - -This file contains the model parameters and getter functions - -Thomas Rose () - -*/ - -#ifndef QCEEQ_PARAM_H -#define QCEEQ_PARAM_H - -namespace EEQ_param { - -constexpr int MAXELEMENT = 104; // 103 + dummy - -// Element-specific electronegativity for the EEQ charges. -static const double xi[MAXELEMENT]{ - +0.00000000, // dummy - +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, +1.06598299, +0.68184178, +1.04581665, +1.09888688, - +1.07206461, +1.09821942, +1.10900303, +1.01039812, +1.00095966, +1.11003303, - +1.16831853, +1.00887482, +1.05928842, +1.07672363, +1.11308426, +1.14340090, - +1.13714110, -}; - -// Element-specific chemical hardnesses for the EEQ charges. -static const double gam[MAXELEMENT]{ - +0.00000000, // dummy - -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.10034731, +0.15801991, -0.00071039, -0.00170887, - -0.00133327, -0.00104386, -0.00094936, -0.00111390, -0.00125257, -0.00095936, - -0.00102814, -0.00104450, -0.00112666, -0.00101529, -0.00059592, -0.00012585, - -0.00140896, -}; - -// Element-specific CN scaling constant for the EEQ charges. -static const double kappa[MAXELEMENT]{ - +0.00000000, // dummy - +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.00013506, -0.00020631, +0.00473118, +0.01590519, - +0.00369763, +0.00417543, +0.00706682, +0.00488679, +0.00505103, +0.00710682, - +0.00463050, +0.00387799, +0.00296795, +0.00400648, +0.00548481, +0.01350400, - +0.00675380, -}; - -static const double alp[MAXELEMENT]{ - +0.00000000, // dummy - +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, +1.04059593, +0.60550051, +1.22262145, +1.28736399, - +1.44431317, +1.29032833, +1.41009404, +1.25501213, +1.15181468, +1.42010424, - +1.43955530, +1.28565237, +1.35017463, +1.33011749, +1.30745135, +1.26526071, - +1.34071499, -}; - -} - -#endif // QCEEQ_PARAM_H \ No newline at end of file diff --git a/src/dftd_ncoord.cpp b/src/dftd_ncoord.cpp index 37ca84c..adcb321 100644 --- a/src/dftd_ncoord.cpp +++ b/src/dftd_ncoord.cpp @@ -33,14 +33,11 @@ namespace dftd4 { NCoordBase::NCoordBase( - double optional_kcn = 7.5, // defaults for D4 EEQ - double optional_norm_exp = 1.0, // defaults for D4 EEQ - double optional_cutoff = 25.0 // defaults for D4 EEQ -){ - NCoordBase::kcn = optional_kcn; - NCoordBase::norm_exp = optional_norm_exp; - NCoordBase::cutoff = optional_cutoff; -} + 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) +{} /** @@ -531,9 +528,6 @@ int NCoordBase::dncoord_d4( /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// -//NCoordErf::NCoordErf(double optional_kcn, double optional_norm_exp, double optional_cutoff) -// : NCoordBase(optional_kcn, optional_norm_exp, optional_cutoff) {} - double NCoordErf::count_fct(double rr) const { return 0.5 * (1.0 + erf(-kcn * (rr - 1.0))); } From b2a6fc0bd5290c6bde60b9f77269606f24f6a126 Mon Sep 17 00:00:00 2001 From: Thomas Rose Date: Wed, 18 Jun 2025 14:59:42 +0200 Subject: [PATCH 06/12] rename functions for consistency --- include/dftd_ncoord.h | 6 +++--- src/dftd_ncoord.cpp | 10 +++++----- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/include/dftd_ncoord.h b/include/dftd_ncoord.h index 4ef9e6d..717ba5f 100644 --- a/include/dftd_ncoord.h +++ b/include/dftd_ncoord.h @@ -57,7 +57,7 @@ class NCoordBase const TIVector&, const TMatrix&); // Calculate the derivative of the coordination number - int dr_ncoord_base( + int dncoord_base( const TMolecule&, const TIVector&, const TMatrix&); @@ -83,7 +83,7 @@ class NCoordBase // Counting function virtual double count_fct(double) const = 0; // Derivative of the counting function - virtual double dr_count_fct(double) const = 0; + virtual double dcount_fct(double) const = 0; // Constructor NCoordBase(double optional_kcn = 7.5, double optional_norm_exp = 1.0, double optional_cutoff = 25.0); // Virtual destructor @@ -95,7 +95,7 @@ class NCoordErf : public NCoordBase { // erf() based counting function double count_fct(double) const override; // derivative of the erf() based counting function - double dr_count_fct(double) const override; + double dcount_fct(double) const 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){} diff --git a/src/dftd_ncoord.cpp b/src/dftd_ncoord.cpp index adcb321..3a2bc85 100644 --- a/src/dftd_ncoord.cpp +++ b/src/dftd_ncoord.cpp @@ -305,7 +305,7 @@ int NCoordBase::get_ncoord( int info; if (lgrad) { - info = dr_ncoord_base(mol, realIdx, dist); + info = dncoord_base(mol, realIdx, dist); } else { info = ncoord_base(mol, realIdx, dist); } @@ -350,7 +350,7 @@ int NCoordBase::ncoord_base( return EXIT_SUCCESS; } -int NCoordBase::dr_ncoord_base( +int NCoordBase::dncoord_base( const TMolecule &mol, const TIVector &realIdx, const TMatrix &dist @@ -385,7 +385,7 @@ int NCoordBase::dr_ncoord_base( cn(ii) += countf; cn(jj) += countf; - dcountf = dr_count_fct(rr) / rcovij; + dcountf = dcount_fct(rr) / rcovij; dcndr(jj, 3 * jj) += dcountf * rx; dcndr(jj, 3 * jj + 1) += dcountf * ry; dcndr(jj, 3 * jj + 2) += dcountf * rz; @@ -507,7 +507,7 @@ int NCoordBase::dncoord_d4( cn(ii) += countf; cn(jj) += countf; - dcountf = den * dr_count_fct(rr) / rcovij; + dcountf = den * dcount_fct(rr) / rcovij; dcndr(3 * jj, jj) += dcountf * rx; dcndr(3 * jj + 1, jj) += dcountf * ry; dcndr(3 * jj + 2, jj) += dcountf * rz; @@ -532,7 +532,7 @@ double NCoordErf::count_fct(double rr) const { return 0.5 * (1.0 + erf(-kcn * (rr - 1.0))); } -double NCoordErf::dr_count_fct(double rr) const { +double NCoordErf::dcount_fct(double rr) const { return -kcn * hlfosqrtpi * exp(-pow(kcn * (rr - 1.0), 2)); } From fc8baed7225ff5f9f84797530b0fb528d300c5de Mon Sep 17 00:00:00 2001 From: Thomas Rose Date: Wed, 18 Jun 2025 15:09:41 +0200 Subject: [PATCH 07/12] Leave deletion of arrays to object destructor --- src/dftd_dispersion.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/dftd_dispersion.cpp b/src/dftd_dispersion.cpp index 717fb61..3294cd5 100644 --- a/src/dftd_dispersion.cpp +++ b/src/dftd_dispersion.cpp @@ -213,7 +213,6 @@ int get_dispersion( ); if (info != EXIT_SUCCESS) return info; - ncoord_erf.cn.Delete(); q.Delete(); refq.Delete(); @@ -250,7 +249,6 @@ int get_dispersion( ); if (info != EXIT_SUCCESS) return info; } else { - ncoord_erf.cn.Delete(); q.Delete(); refq.Delete(); gwvec.Delete(); @@ -265,7 +263,6 @@ int get_dispersion( if (lgrad) { BLAS_Add_Mat_x_Vec(gradient, ncoord_erf.dcndr, dEdcn, false, 1.0); } - ncoord_erf.dcndr.DelMat(); dEdcn.DelVec(); dEdq.DelVec(); From 5b2233968c7bf498fe46fc4de97013c81f389b65 Mon Sep 17 00:00:00 2001 From: Thomas Rose Date: Fri, 18 Jul 2025 11:35:00 +0200 Subject: [PATCH 08/12] Clean up NCoordBase member variables in destructor --- include/dftd_ncoord.h | 7 ++++++- src/dftd_eeq.cpp | 3 --- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/include/dftd_ncoord.h b/include/dftd_ncoord.h index 717ba5f..9095045 100644 --- a/include/dftd_ncoord.h +++ b/include/dftd_ncoord.h @@ -87,7 +87,10 @@ class NCoordBase // Constructor NCoordBase(double optional_kcn = 7.5, double optional_norm_exp = 1.0, double optional_cutoff = 25.0); // Virtual destructor - virtual ~NCoordBase() = default; + virtual ~NCoordBase() { + cn.DelVec(); + dcndr.DelMat(); + } }; class NCoordErf : public NCoordBase { @@ -99,6 +102,8 @@ class NCoordErf : public NCoordBase { // 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; }; /** diff --git a/src/dftd_eeq.cpp b/src/dftd_eeq.cpp index 4c5e70e..1ac4e98 100644 --- a/src/dftd_eeq.cpp +++ b/src/dftd_eeq.cpp @@ -83,9 +83,6 @@ int ChargeModel::get_charges( eeq_chrgeq(mol, realIdx, dist, charge, ncoord_erf.cn, q, ncoord_erf.dcndr, dqdr, lgrad, lverbose); if (info != EXIT_SUCCESS) return info; - ncoord_erf.dcndr.DelMat(); - ncoord_erf.cn.DelVec(); - return EXIT_SUCCESS; }; From 2afd44ea9fbffef549afae2dc7a04b5680351009 Mon Sep 17 00:00:00 2001 From: Thomas Rose Date: Tue, 22 Jul 2025 10:40:14 +0200 Subject: [PATCH 09/12] NCoord refactor and bug fixes - Merged ncoord_d4 and ncoord_base - Fixed dcndr calculation --- include/dftd_cblas.h | 6 +- include/dftd_ncoord.h | 103 +++++----------------- src/dftd_dispersion.cpp | 6 +- src/dftd_eeq.cpp | 14 +-- src/dftd_ncoord.cpp | 174 ++++++++++--------------------------- test/test_ghost.cpp | 7 +- test/test_grad.cpp | 107 +++++++++++++++++++++++ test/test_ncoord.cpp | 184 +++++++++++++++++++++++++++++++++++++++- 8 files changed, 375 insertions(+), 226 deletions(-) diff --git a/include/dftd_cblas.h b/include/dftd_cblas.h index a9ed3a3..8601e4d 100644 --- a/include/dftd_cblas.h +++ b/include/dftd_cblas.h @@ -41,8 +41,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, @@ -58,7 +58,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 9095045..d354793 100644 --- a/include/dftd_ncoord.h +++ b/include/dftd_ncoord.h @@ -80,10 +80,14 @@ class NCoordBase const TIVector&, const TMatrix&); + // Electronegativity factor + virtual double get_en_factor(int, int) const; // Counting function virtual double count_fct(double) const = 0; // Derivative of the counting function virtual double dcount_fct(double) const = 0; + // Soft maximum/cutoff for coordination number + virtual int cut_coordination_number(const double, TVector&, TMatrix&, bool) = 0; // Constructor NCoordBase(double optional_kcn = 7.5, double optional_norm_exp = 1.0, double optional_cutoff = 25.0); // Virtual destructor @@ -99,6 +103,8 @@ class NCoordErf : public NCoordBase { 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){} @@ -106,6 +112,23 @@ class NCoordErf : public NCoordBase { ~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. * @@ -216,86 +239,6 @@ extern int dncoord_erf( /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// -/** - * 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. * diff --git a/src/dftd_dispersion.cpp b/src/dftd_dispersion.cpp index 3294cd5..7ad0d47 100644 --- a/src/dftd_dispersion.cpp +++ b/src/dftd_dispersion.cpp @@ -105,7 +105,7 @@ int get_dispersion( TVector q; // partial charges from EEQ model TMatrix dqdr; // derivative of partial charges TVector gradient; // derivative of dispersion energy - NCoordErf ncoord_erf; + NCoordErfD4 ncoord_erf; multicharge::EEQModel chrg_model; // Charge model q.NewVector(nat); @@ -119,7 +119,7 @@ int get_dispersion( if (info != EXIT_SUCCESS) return info; // get the D4 coordination number - info = ncoord_erf.get_ncoord_d4(mol, realIdx, dist, lgrad); + info = ncoord_erf.get_ncoord(mol, realIdx, dist, lgrad); if (info != EXIT_SUCCESS) return info; // maximum number of reference systems @@ -261,7 +261,7 @@ int get_dispersion( dc6dcn.DelMat(); dc6dq.DelMat(); - if (lgrad) { BLAS_Add_Mat_x_Vec(gradient, ncoord_erf.dcndr, dEdcn, false, 1.0); } + if (lgrad) BLAS_Add_Mat_x_Vec(gradient, ncoord_erf.dcndr, dEdcn, true, 1.0); dEdcn.DelVec(); dEdq.DelVec(); diff --git a/src/dftd_eeq.cpp b/src/dftd_eeq.cpp index 1ac4e98..22636eb 100644 --- a/src/dftd_eeq.cpp +++ b/src/dftd_eeq.cpp @@ -187,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); } } @@ -357,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 3a2bc85..cd4975c 100644 --- a/src/dftd_ncoord.cpp +++ b/src/dftd_ncoord.cpp @@ -324,6 +324,7 @@ int NCoordBase::ncoord_base( ) { double r = 0.0, rcovij = 0.0, rr = 0.0; double countf = 0.0; + double f_en; // initialize cn to zero int nat = realIdx.Max() + 1; cn.NewVector(nat); @@ -341,7 +342,8 @@ int NCoordBase::ncoord_base( rcovij = rad[mol.ATNO(i)] + rad[mol.ATNO(j)]; rr = r / rcovij; - countf = count_fct(rr); + f_en = get_en_factor(mol.ATNO(i), mol.ATNO(j)); + countf = f_en * count_fct(rr); cn(ii) += countf; cn(jj) += countf; } @@ -350,6 +352,7 @@ int NCoordBase::ncoord_base( return EXIT_SUCCESS; } + int NCoordBase::dncoord_base( const TMolecule &mol, const TIVector &realIdx, @@ -358,6 +361,7 @@ int NCoordBase::dncoord_base( 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); @@ -374,33 +378,34 @@ int NCoordBase::dncoord_base( 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 = count_fct(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 = dcount_fct(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; } } @@ -410,119 +415,12 @@ int NCoordBase::dncoord_base( /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// -int NCoordBase::get_ncoord_d4( - const TMolecule &mol, - const TMatrix &dist, - bool lgrad -) { - TIVector realIdx; - initializeRealIdx(mol.NAtoms, realIdx); - - return get_ncoord_d4(mol, realIdx, dist, lgrad); -}; - -int NCoordBase::get_ncoord_d4( - const TMolecule &mol, - const TIVector &realIdx, - const TMatrix &dist, - bool lgrad -) { - if (lgrad) return dncoord_d4(mol, realIdx, dist); - return ncoord_d4(mol, realIdx, dist); -}; - -int NCoordBase::ncoord_d4( - const TMolecule &mol, - const TIVector &realIdx, - const TMatrix &dist -) { - double r = 0.0, rcovij = 0.0, rr = 0.0; - double den = 0.0; - double countf = 0.0; - int izp, jzp; - // 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); - 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 * count_fct(rr); - - cn(ii) += countf; - cn(jj) += countf; - } - } - return EXIT_SUCCESS; +double NCoordBase::get_en_factor(int i, int j) const { + return 1.0; } -int NCoordBase::dncoord_d4( - const TMolecule &mol, - const TIVector &realIdx, - 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, den = 0.0; - int izp, jzp; - // initialize cn and dcndr to zero - int nat = realIdx.Max() + 1; - cn.NewVector(nat); - dcndr.NewMatrix(3*nat, nat); - - 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 * count_fct(rr); - cn(ii) += countf; - cn(jj) += countf; - - dcountf = den * dcount_fct(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); } /////////////////////////////////////////////////////////////////////////////// @@ -536,20 +434,29 @@ 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; } } } @@ -569,4 +476,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 f402891..401010c 100644 --- a/test/test_ghost.cpp +++ b/test/test_ghost.cpp @@ -59,8 +59,8 @@ int test_water( // COORDINATION NUMBER CHECK // erf-CN without cutoff - NCoordErf ncoord_erf(7.5, 1.0, 9999.9); - info = ncoord_erf.get_ncoord_d4(mol, realIdx, dist, 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 @@ -146,7 +146,7 @@ int test_water( d4grad[i] = 0.0; } - info = ncoord_erf.get_ncoord_d4(mol, realIdx, dist, lgrad); + info = ncoord_erf.get_ncoord(mol, realIdx, dist, lgrad); if (info != EXIT_SUCCESS) return info; dqdr.NewMatrix(3 * nat, nat); @@ -159,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 21ef6ce..9674874 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; } } @@ -166,9 +169,113 @@ 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; + TVector q_dum; + 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, dqdr, true); + if (info != EXIT_SUCCESS) return info; + for (int i = 0; i < mol.NAtoms; i++) { + for (int c = 0; c < 3; c++) { + for (int k = 0; k < mol.NAtoms; k++) { + analytic_dqdr(3 * i + c, k) = dqdr(3 * i + c, k); + } + } + } + + // numerical gradient + for (int i = 0; i < mol.NAtoms; i++) { + for (int c = 0; c < 3; c++) { + multicharge::EEQModel eeq_model_l; + multicharge::EEQModel eeq_model_r; + q_r.NewVec(n); + q_l.NewVec(n); + q_dum.NewVec(mol.NAtoms); + mol.CC(i, c) += step; + calc_distances(mol, realIdx, dist); + eeq_model_r.get_charges(mol, realIdx, dist, 0, cutoff.cn_eeq, q_dum, dqdr, false); + + for (int k = 0; k < mol.NAtoms; k++) { + q_r(k) = q_dum(k); + } + + mol.CC(i, c) = mol.CC(i, c) - 2 * step; + calc_distances(mol, realIdx, dist); + eeq_model_l.get_charges(mol, realIdx, dist, 0, cutoff.cn_eeq, q_dum, dqdr, false); + for (int k = 0; k < mol.NAtoms; k++) { + q_l(k) = q_dum(k); + } + + 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 7a855f7..21cd503 100644 --- a/test/test_ncoord.cpp +++ b/test/test_ncoord.cpp @@ -61,8 +61,8 @@ int test_cn( calc_distances(mol, realIdx, dist); // erf-CN without cutoff - NCoordErf ncoord_erf(7.5, 1.0, 9999.9); - info = ncoord_erf.get_ncoord_d4(mol, realIdx, dist, 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 @@ -76,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; @@ -91,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; } From 3c7eed979506375f7d04bb98ff5f23323dac3562 Mon Sep 17 00:00:00 2001 From: Thomas3R Date: Wed, 13 Aug 2025 09:28:44 +0200 Subject: [PATCH 10/12] Update src/dftd_ncoord.cpp Co-authored-by: Marvin Friede <51965259+marvinfriede@users.noreply.github.com> --- src/dftd_ncoord.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dftd_ncoord.cpp b/src/dftd_ncoord.cpp index cd4975c..6dfca69 100644 --- a/src/dftd_ncoord.cpp +++ b/src/dftd_ncoord.cpp @@ -324,7 +324,7 @@ int NCoordBase::ncoord_base( ) { double r = 0.0, rcovij = 0.0, rr = 0.0; double countf = 0.0; - double f_en; + double f_en = 0.0; // initialize cn to zero int nat = realIdx.Max() + 1; cn.NewVector(nat); From f6a37e7f24d6494185e5ee8bde2b0c2fcb29c312 Mon Sep 17 00:00:00 2001 From: Thomas Rose Date: Wed, 13 Aug 2025 10:41:01 +0200 Subject: [PATCH 11/12] Address PR review comments --- include/dftd_ncoord.h | 249 ++++++++++++++---------------------------- src/dftd_ncoord.cpp | 2 - test/test_grad.cpp | 32 ++---- 3 files changed, 91 insertions(+), 192 deletions(-) diff --git a/include/dftd_ncoord.h b/include/dftd_ncoord.h index d354793..c174b16 100644 --- a/include/dftd_ncoord.h +++ b/include/dftd_ncoord.h @@ -38,56 +38,92 @@ class NCoordBase double kcn; // Steepness of counting function double norm_exp; double cutoff; - // Get the coordination number + /** + * 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&, - const TMatrix&, - const double, - TVector&, - TMatrix&, - bool); + 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&, - const TIVector&, - const TMatrix&, - bool); - // Calculate the coordination number using the virtual counting function + 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&, - const TIVector&, - const TMatrix&); - // Calculate the derivative of the coordination number + 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&, - const TIVector&, - const TMatrix&); - // Get the DFT-D4 coordination number - int get_ncoord_d4( - const TMolecule&, - const TMatrix&, - bool); - int get_ncoord_d4( - const TMolecule&, - const TIVector&, - const TMatrix&, - bool); - int ncoord_d4( - const TMolecule&, - const TIVector&, - const TMatrix&); - int dncoord_d4( - const TMolecule&, - const TIVector&, - const TMatrix&); - - // Electronegativity factor - virtual double get_en_factor(int, int) const; - // Counting function - virtual double count_fct(double) const = 0; - // Derivative of the counting function - virtual double dcount_fct(double) const = 0; - // Soft maximum/cutoff for coordination number - virtual int cut_coordination_number(const double, TVector&, TMatrix&, bool) = 0; + 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 @@ -152,127 +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 -); - -/////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////// - -/** - * 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_ncoord.cpp b/src/dftd_ncoord.cpp index cd4975c..cd84812 100644 --- a/src/dftd_ncoord.cpp +++ b/src/dftd_ncoord.cpp @@ -286,8 +286,6 @@ int NCoordBase::get_ncoord( const TMolecule &mol, const TMatrix &dist, const double cutoff, - TVector &cn, - TMatrix &dcndr, bool lgrad ) { TIVector realIdx; diff --git a/test/test_grad.cpp b/test/test_grad.cpp index 9674874..fb04a9f 100644 --- a/test/test_grad.cpp +++ b/test/test_grad.cpp @@ -192,7 +192,6 @@ int test_numgrad_dqdr( analytic_dqdr.NewMat(3 * mol.NAtoms, mol.NAtoms); multicharge::EEQModel eeq_model; TVector q; - TVector q_dum; TMatrix dqdr; q.NewVec(mol.NAtoms); dqdr.NewMat(3 * mol.NAtoms, mol.NAtoms); @@ -212,39 +211,26 @@ int test_numgrad_dqdr( // analytical gradient calc_distances(mol, realIdx, dist); info = - eeq_model.get_charges(mol, realIdx, dist, 0, cutoff.cn_eeq, q, dqdr, true); + eeq_model.get_charges(mol, realIdx, dist, 0, cutoff.cn_eeq, q, analytic_dqdr, true); if (info != EXIT_SUCCESS) return info; - for (int i = 0; i < mol.NAtoms; i++) { - for (int c = 0; c < 3; c++) { - for (int k = 0; k < mol.NAtoms; k++) { - analytic_dqdr(3 * i + c, k) = dqdr(3 * i + c, k); - } - } - } - // numerical gradient + // calculate numerical gradient via finite difference method + multicharge::EEQModel eeq_model_num; // charge model instance for numerical gradient for (int i = 0; i < mol.NAtoms; i++) { for (int c = 0; c < 3; c++) { - multicharge::EEQModel eeq_model_l; - multicharge::EEQModel eeq_model_r; + // calculate forward point q_r.NewVec(n); - q_l.NewVec(n); - q_dum.NewVec(mol.NAtoms); mol.CC(i, c) += step; calc_distances(mol, realIdx, dist); - eeq_model_r.get_charges(mol, realIdx, dist, 0, cutoff.cn_eeq, q_dum, dqdr, false); - - for (int k = 0; k < mol.NAtoms; k++) { - q_r(k) = q_dum(k); - } + eeq_model_num.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_l.get_charges(mol, realIdx, dist, 0, cutoff.cn_eeq, q_dum, dqdr, false); - for (int k = 0; k < mol.NAtoms; k++) { - q_l(k) = q_dum(k); - } + eeq_model_num.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; From 3f8cb39a2f7e81875c8abd0361f8b8de457b154d Mon Sep 17 00:00:00 2001 From: Thomas Rose Date: Wed, 13 Aug 2025 13:06:47 +0200 Subject: [PATCH 12/12] remove reduntant charge model intance from numerical gradient test --- test/test_grad.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/test/test_grad.cpp b/test/test_grad.cpp index fb04a9f..d2ecdd5 100644 --- a/test/test_grad.cpp +++ b/test/test_grad.cpp @@ -215,20 +215,19 @@ int test_numgrad_dqdr( if (info != EXIT_SUCCESS) return info; // calculate numerical gradient via finite difference method - multicharge::EEQModel eeq_model_num; // charge model instance for numerical gradient 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_num.get_charges(mol, realIdx, dist, 0, cutoff.cn_eeq, q_r, dqdr, false); + 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_num.get_charges(mol, realIdx, dist, 0, cutoff.cn_eeq, q_l, dqdr, false); + 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;