From 86c2ae6dabec5d2528a88f5f60229acdd0b399ae Mon Sep 17 00:00:00 2001 From: Trung Nguyen Date: Fri, 9 Dec 2022 23:34:49 -0600 Subject: [PATCH] Switched to using q_scaled, keeping q as the real, unscaled charges --- src/DIELECTRIC/atom_vec_dielectric.cpp | 43 ++++++------------- src/DIELECTRIC/atom_vec_dielectric.h | 4 +- src/DIELECTRIC/pair_coul_cut_dielectric.cpp | 7 +-- src/DIELECTRIC/pair_coul_long_dielectric.cpp | 2 +- .../pair_lj_cut_coul_cut_dielectric.cpp | 7 +-- .../pair_lj_cut_coul_debye_dielectric.cpp | 7 +-- .../pair_lj_cut_coul_long_dielectric.cpp | 14 +++--- .../pair_lj_cut_coul_msm_dielectric.cpp | 24 ++++++++--- .../pair_lj_long_coul_long_dielectric.cpp | 5 ++- src/atom.cpp | 6 +-- src/atom.h | 2 +- src/set.cpp | 7 +-- 12 files changed, 59 insertions(+), 69 deletions(-) diff --git a/src/DIELECTRIC/atom_vec_dielectric.cpp b/src/DIELECTRIC/atom_vec_dielectric.cpp index 23dee2f38d..18fcac8d1c 100644 --- a/src/DIELECTRIC/atom_vec_dielectric.cpp +++ b/src/DIELECTRIC/atom_vec_dielectric.cpp @@ -60,28 +60,28 @@ AtomVecDielectric::AtomVecDielectric(LAMMPS *_lmp) : AtomVec(_lmp) "angle_atom1", "angle_atom2", "angle_atom3", "num_dihedral", "dihedral_type", "dihedral_atom1", "dihedral_atom2", "dihedral_atom3", "dihedral_atom4", "num_improper", "improper_type", "improper_atom1", "improper_atom2", "improper_atom3", "improper_atom4", "nspecial", "special", - "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"}; + "mu", "area", "ed", "em", "epsilon", "curvature", "q_scaled"}; fields_copy = {"q", "molecule", "num_bond", "bond_type", "bond_atom", "num_angle", "angle_type", "angle_atom1", "angle_atom2", "angle_atom3", "num_dihedral", "dihedral_type", "dihedral_atom1", "dihedral_atom2", "dihedral_atom3", "dihedral_atom4", "num_improper", "improper_type", "improper_atom1", "improper_atom2", "improper_atom3", "improper_atom4", "nspecial", "special", - "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"}; - fields_comm = {"q", "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"}; - fields_border = {"q", "molecule", "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"}; + "mu", "area", "ed", "em", "epsilon", "curvature", "q_scaled"}; + fields_comm = {"q", "mu", "area", "ed", "em", "epsilon", "curvature", "q_scaled"}; + fields_border = {"q", "molecule", "mu", "area", "ed", "em", "epsilon", "curvature", "q_scaled"}; fields_border_vel = {"q", "molecule", "mu", "area", "ed", "em", "epsilon", "curvature", - "q_unscaled"}; + "q_scaled"}; fields_exchange = {"q", "molecule", "num_bond", "bond_type", "bond_atom", "num_angle", "angle_type", "angle_atom1", "angle_atom2", "angle_atom3", "num_dihedral", "dihedral_type", "dihedral_atom1", "dihedral_atom2", "dihedral_atom3", "dihedral_atom4", "num_improper", "improper_type", "improper_atom1", "improper_atom2", "improper_atom3", "improper_atom4", - "nspecial", "special", "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"}; + "nspecial", "special", "mu", "area", "ed", "em", "epsilon", "curvature", "q_scaled"}; fields_restart = {"q", "molecule", "num_bond", "bond_type", "bond_atom", "num_angle", "angle_type", "angle_atom1", "angle_atom2", "angle_atom3", "num_dihedral", "dihedral_type", "dihedral_atom1", "dihedral_atom2", "dihedral_atom3", "dihedral_atom4", "num_improper", "improper_type", "improper_atom1", "improper_atom2", "improper_atom3", "improper_atom4", - "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"}; + "mu", "area", "ed", "em", "epsilon", "curvature", "q_scaled"}; fields_create = {"q", "molecule", "num_bond", "num_angle", "num_dihedral", "num_improper", - "nspecial", "mu", "area", "ed", "em", "epsilon", "curvature", "q_unscaled"}; + "nspecial", "mu", "area", "ed", "em", "epsilon", "curvature", "q_scaled"}; fields_data_atom = {"id", "molecule", "type", "q", "x", "mu3", "area", "ed", "em", "epsilon", "curvature"}; fields_data_vel = {"id", "v"}; @@ -151,7 +151,7 @@ void AtomVecDielectric::grow_pointers() em = atom->em; epsilon = atom->epsilon; curvature = atom->curvature; - q_unscaled = atom->q_unscaled; + q_scaled = atom->q_scaled; } /* ---------------------------------------------------------------------- @@ -181,31 +181,12 @@ void AtomVecDielectric::data_atom_post(int ilocal) nspecial[ilocal][2] = 0; double *q = atom->q; - q_unscaled[ilocal] = q[ilocal]; - q[ilocal] /= epsilon[ilocal]; + q_scaled[ilocal] = q[ilocal]/epsilon[ilocal]; double *mu_one = mu[ilocal]; mu_one[3] = sqrt(mu_one[0] * mu_one[0] + mu_one[1] * mu_one[1] + mu_one[2] * mu_one[2]); } -/* ---------------------------------------------------------------------- - restore original data for writing the data file -------------------------------------------------------------------------- */ - -void AtomVecDielectric::pack_data_pre(int ilocal) -{ - atom->q[ilocal] = q_unscaled[ilocal]; -} - -/* ---------------------------------------------------------------------- - undo restore and get back to post read data state -------------------------------------------------------------------------- */ - -void AtomVecDielectric::pack_data_post(int ilocal) -{ - atom->q[ilocal] /= epsilon[ilocal]; -} - /* ---------------------------------------------------------------------- initialize other atom quantities after AtomVec::unpack_restart() ------------------------------------------------------------------------- */ @@ -229,7 +210,7 @@ int AtomVecDielectric::property_atom(const std::string &name) if (name == "em") return 2; if (name == "epsilon") return 3; if (name == "curvature") return 4; - if (name == "q_unscaled") return 5; + if (name == "q_scaled") return 5; return -1; } @@ -287,7 +268,7 @@ void AtomVecDielectric::pack_property_atom(int index, double *buf, int nvalues, } else if (index == 5) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) - buf[n] = q_unscaled[i]; + buf[n] = q_scaled[i]; else buf[n] = 0.0; n += nvalues; diff --git a/src/DIELECTRIC/atom_vec_dielectric.h b/src/DIELECTRIC/atom_vec_dielectric.h index 77b6001896..9e1bff8c4d 100644 --- a/src/DIELECTRIC/atom_vec_dielectric.h +++ b/src/DIELECTRIC/atom_vec_dielectric.h @@ -38,8 +38,6 @@ class AtomVecDielectric : public AtomVec { void unpack_restart_init(int) override; int property_atom(const std::string &) override; void pack_property_atom(int, double *, int, int) override; - void pack_data_pre(int) override; - void pack_data_post(int) override; protected: int *num_bond, *num_angle, *num_dihedral, *num_improper; @@ -49,7 +47,7 @@ class AtomVecDielectric : public AtomVec { int bond_per_atom, angle_per_atom, dihedral_per_atom, improper_per_atom; double **mu; - double *area, *ed, *em, *epsilon, *curvature, *q_unscaled; + double *area, *ed, *em, *epsilon, *curvature, *q_scaled; }; } // namespace LAMMPS_NS diff --git a/src/DIELECTRIC/pair_coul_cut_dielectric.cpp b/src/DIELECTRIC/pair_coul_cut_dielectric.cpp index 912f0fd2e5..e9d7ac3455 100644 --- a/src/DIELECTRIC/pair_coul_cut_dielectric.cpp +++ b/src/DIELECTRIC/pair_coul_cut_dielectric.cpp @@ -69,7 +69,7 @@ void PairCoulCutDielectric::compute(int eflag, int vflag) double **x = atom->x; double **f = atom->f; - double *q = atom->q; + double *q = atom->q_scaled; double *eps = atom->epsilon; double **norm = atom->mu; double *curvature = atom->curvature; @@ -165,10 +165,11 @@ double PairCoulCutDielectric::single(int i, int j, int /*itype*/, int /*jtype*/, double factor_coul, double /*factor_lj*/, double &fforce) { double r2inv, phicoul, ei, ej; + double *q = atom->q_scaled; double *eps = atom->epsilon; r2inv = 1.0 / rsq; - fforce = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv) * eps[i]; + fforce = force->qqrd2e * q[i] * q[j] * sqrt(r2inv) * eps[i]; double eng = 0.0; if (eps[i] == 1) @@ -179,7 +180,7 @@ double PairCoulCutDielectric::single(int i, int j, int /*itype*/, int /*jtype*/, ej = 0; else ej = eps[j]; - phicoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv); + phicoul = force->qqrd2e * q[i] * q[j] * sqrt(r2inv); phicoul *= 0.5 * (ei + ej); eng += factor_coul * phicoul; diff --git a/src/DIELECTRIC/pair_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_coul_long_dielectric.cpp index 991135199e..da9528f8a2 100644 --- a/src/DIELECTRIC/pair_coul_long_dielectric.cpp +++ b/src/DIELECTRIC/pair_coul_long_dielectric.cpp @@ -73,7 +73,7 @@ void PairCoulLongDielectric::compute(int eflag, int vflag) double **x = atom->x; double **f = atom->f; - double *q = atom->q; + double *q = atom->q_scaled; double *eps = atom->epsilon; double **norm = atom->mu; double *curvature = atom->curvature; diff --git a/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp index 2234ef1e83..afe657191e 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp @@ -74,7 +74,7 @@ void PairLJCutCoulCutDielectric::compute(int eflag, int vflag) double **x = atom->x; double **f = atom->f; - double *q = atom->q; + double *q = atom->q_scaled; double *eps = atom->epsilon; double **norm = atom->mu; double *curvature = atom->curvature; @@ -195,11 +195,12 @@ double PairLJCutCoulCutDielectric::single(int i, int j, int itype, int jtype, do double factor_coul, double factor_lj, double &fforce) { double r2inv, r6inv, forcecoul, forcelj, phicoul, ei, ej, philj; + double *q = atom->q_scaled; double *eps = atom->epsilon; r2inv = 1.0 / rsq; if (rsq < cut_coulsq[itype][jtype]) - forcecoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv) * eps[i]; + forcecoul = force->qqrd2e * q[i] * q[j] * sqrt(r2inv) * eps[i]; else forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { @@ -219,7 +220,7 @@ double PairLJCutCoulCutDielectric::single(int i, int j, int itype, int jtype, do else ej = eps[j]; if (rsq < cut_coulsq[itype][jtype]) { - phicoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv); + phicoul = force->qqrd2e * q[i] * q[j] * sqrt(r2inv); phicoul *= 0.5 * (ei + ej); eng += factor_coul * phicoul; } diff --git a/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp index 46b032a369..3f39206251 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp @@ -75,7 +75,7 @@ void PairLJCutCoulDebyeDielectric::compute(int eflag, int vflag) double **x = atom->x; double **f = atom->f; - double *q = atom->q; + double *q = atom->q_scaled; double *eps = atom->epsilon; double **norm = atom->mu; double *curvature = atom->curvature; @@ -197,6 +197,7 @@ double PairLJCutCoulDebyeDielectric::single(int i, int j, int itype, int jtype, { double r2inv, r6inv, forcecoul, forcelj, phicoul, ei, ej, philj; double r, rinv, screening; + double *q = atom->q_scaled; double *eps = atom->epsilon; r2inv = 1.0 / rsq; @@ -204,7 +205,7 @@ double PairLJCutCoulDebyeDielectric::single(int i, int j, int itype, int jtype, r = sqrt(rsq); rinv = 1.0 / r; screening = exp(-kappa * r); - forcecoul = force->qqrd2e * atom->q[i] * atom->q[j] * screening * (kappa + rinv) * eps[i]; + forcecoul = force->qqrd2e * q[i] * q[j] * screening * (kappa + rinv) * eps[i]; } else forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { @@ -224,7 +225,7 @@ double PairLJCutCoulDebyeDielectric::single(int i, int j, int itype, int jtype, else ej = eps[j]; if (rsq < cut_coulsq[itype][jtype]) { - phicoul = force->qqrd2e * atom->q[i] * atom->q[j] * rinv * screening; + phicoul = force->qqrd2e * q[i] * q[j] * rinv * screening; phicoul *= 0.5 * (ei + ej); eng += factor_coul * phicoul; } diff --git a/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp index 2bd02b500e..a99699883a 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp @@ -81,7 +81,7 @@ void PairLJCutCoulLongDielectric::compute(int eflag, int vflag) double **x = atom->x; double **f = atom->f; - double *q = atom->q; + double *q = atom->q_scaled; double *eps = atom->epsilon; double **norm = atom->mu; double *curvature = atom->curvature; @@ -252,6 +252,7 @@ double PairLJCutCoulLongDielectric::single(int i, int j, int itype, int jtype, d double r2inv, r6inv, r, grij, expm2, t, erfc, ei, ej, prefactor; double fraction, table, forcecoul, forcelj, phicoul, philj; int itable; + double *q = atom->q_scaled; double *eps = atom->epsilon; r2inv = 1.0 / rsq; @@ -262,7 +263,7 @@ double PairLJCutCoulLongDielectric::single(int i, int j, int itype, int jtype, d expm2 = exp(-grij * grij); t = 1.0 / (1.0 + EWALD_P * grij); erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2; - prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r; + prefactor = force->qqrd2e * q[i] * q[j] / r; forcecoul = prefactor * (erfc + EWALD_F * grij * expm2); if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; } else { @@ -272,10 +273,10 @@ double PairLJCutCoulLongDielectric::single(int i, int j, int itype, int jtype, d itable >>= ncoulshiftbits; fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable]; table = ftable[itable] + fraction * dftable[itable]; - forcecoul = atom->q[i] * atom->q[j] * table; + forcecoul = q[i] * q[j] * table; if (factor_coul < 1.0) { table = ctable[itable] + fraction * dctable[itable]; - prefactor = atom->q[i] * atom->q[j] * table; + prefactor = q[i] * q[j] * table; forcecoul -= (1.0 - factor_coul) * prefactor; } } @@ -301,12 +302,11 @@ double PairLJCutCoulLongDielectric::single(int i, int j, int itype, int jtype, d ej = eps[j]; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) - phicoul = prefactor * (ei + ej) * erfc; + phicoul = prefactor * 0.5 * (ei + ej) * erfc; else { table = etable[itable] + fraction * detable[itable]; - phicoul = atom->q[i] * atom->q[j] * (ei + ej) * table; + phicoul = q[i] * q[j] * 0.5 * (ei + ej) * table; } - phicoul *= 0.5; if (factor_coul < 1.0) phicoul -= (1.0 - factor_coul) * prefactor; eng += phicoul; } diff --git a/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp b/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp index d5cedd7713..ade05ddb04 100644 --- a/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp @@ -101,7 +101,7 @@ void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag) double **x = atom->x; double **f = atom->f; - double *q = atom->q; + double *q = atom->q_scaled; double *eps = atom->epsilon; double **norm = atom->mu; double *curvature = atom->curvature; @@ -269,15 +269,17 @@ void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag) double PairLJCutCoulMSMDielectric::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_lj, double &fforce) { - double r2inv, r6inv, r, egamma, fgamma, prefactor; + double r2inv, r6inv, r, egamma, fgamma, ei, ej, prefactor; double fraction, table, forcecoul, forcelj, phicoul, philj; int itable; + double *q = atom->q_scaled; + double *eps = atom->epsilon; r2inv = 1.0 / rsq; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) { r = sqrt(rsq); - prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r; + prefactor = force->qqrd2e * q[i] * q[j] / r; egamma = 1.0 - (r / cut_coul) * force->kspace->gamma(r / cut_coul); fgamma = 1.0 + (rsq / cut_coulsq) * force->kspace->dgamma(r / cut_coul); forcecoul = prefactor * fgamma; @@ -289,10 +291,10 @@ double PairLJCutCoulMSMDielectric::single(int i, int j, int itype, int jtype, do itable >>= ncoulshiftbits; fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable]; table = ftable[itable] + fraction * dftable[itable]; - forcecoul = atom->q[i] * atom->q[j] * table; + forcecoul = q[i] * q[j] * table; if (factor_coul < 1.0) { table = ctable[itable] + fraction * dctable[itable]; - prefactor = atom->q[i] * atom->q[j] * table; + prefactor = q[i] * q[j] * table; forcecoul -= (1.0 - factor_coul) * prefactor; } } @@ -308,12 +310,20 @@ double PairLJCutCoulMSMDielectric::single(int i, int j, int itype, int jtype, do fforce = (forcecoul + factor_lj * forcelj) * r2inv; double eng = 0.0; + if (eps[i] == 1) + ei = 0; + else + ei = eps[i]; + if (eps[j] == 1) + ej = 0; + else + ej = eps[j]; if (rsq < cut_coulsq) { if (!ncoultablebits || rsq <= tabinnersq) - phicoul = prefactor * egamma; + phicoul = prefactor * 0.5 * (ei + ej) * egamma; else { table = etable[itable] + fraction * detable[itable]; - phicoul = atom->q[i] * atom->q[j] * table; + phicoul = q[i] * q[j] * 0.5 * (ei + ej) * table; } if (factor_coul < 1.0) phicoul -= (1.0 - factor_coul) * prefactor; eng += phicoul; diff --git a/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp b/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp index 189bc78f20..1de5e8afcf 100644 --- a/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp +++ b/src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp @@ -92,7 +92,7 @@ void PairLJLongCoulLongDielectric::compute(int eflag, int vflag) double **x = atom->x, *x0 = x[0]; double **f = atom->f; - double *q = atom->q; + double *q = atom->q_scaled; double *eps = avec->epsilon; double **norm = avec->mu; double *curvature = avec->curvature; @@ -292,7 +292,8 @@ double PairLJLongCoulLongDielectric::single(int i, int j, int itype, int jtype, double factor_coul, double factor_lj, double &fforce) { double r2inv, r6inv, force_coul, force_lj, ei, ej; - double g2 = g_ewald_6 * g_ewald_6, g6 = g2 * g2 * g2, g8 = g6 * g2, *q = atom->q; + double g2 = g_ewald_6 * g_ewald_6, g6 = g2 * g2 * g2, g8 = g6 * g2; + double *q = atom->q_scaled; double *eps = avec->epsilon; double eng = 0.0; diff --git a/src/atom.cpp b/src/atom.cpp index 3d843cbbbc..0fde1d3f8b 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -214,7 +214,7 @@ Atom::Atom(LAMMPS *_lmp) : Pointers(_lmp) // DIELECTRIC package - area = ed = em = epsilon = curvature = q_unscaled = nullptr; + area = ed = em = epsilon = curvature = q_scaled = nullptr; // end of customization section // -------------------------------------------------------------------- @@ -563,7 +563,7 @@ void Atom::peratom_create() add_peratom("em",&em,DOUBLE,0); add_peratom("epsilon",&epsilon,DOUBLE,0); add_peratom("curvature",&curvature,DOUBLE,0); - add_peratom("q_unscaled",&q_unscaled,DOUBLE,0); + add_peratom("q_scaled",&q_scaled,DOUBLE,0); // end of customization section // -------------------------------------------------------------------- @@ -2950,7 +2950,7 @@ void *Atom::extract(const char *name) if (strcmp(name,"em") == 0) return (void *) em; if (strcmp(name,"epsilon") == 0) return (void *) epsilon; if (strcmp(name,"curvature") == 0) return (void *) curvature; - if (strcmp(name,"q_unscaled") == 0) return (void *) q_unscaled; + if (strcmp(name,"q_scaled") == 0) return (void *) q_scaled; // end of customization section // -------------------------------------------------------------------- diff --git a/src/atom.h b/src/atom.h index 70e064c0ff..61ae0d047b 100644 --- a/src/atom.h +++ b/src/atom.h @@ -172,7 +172,7 @@ class Atom : protected Pointers { // DIELECTRIC package - double *area, *ed, *em, *epsilon, *curvature, *q_unscaled; + double *area, *ed, *em, *epsilon, *curvature, *q_scaled; // end of customization section // -------------------------------------------------------------------- diff --git a/src/set.cpp b/src/set.cpp index 0131b1e16c..b45b6e5798 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -1087,14 +1087,11 @@ void Set::set(int keyword) else if (keyword == EPSILON) { if (dvalue >= 0.0) { - // compute the unscaled charge value (i.e. atom->q_unscaled) // assign the new local dielectric constant - // update both the scaled and unscaled charge values + // update both the scaled charge value - double q_unscaled = atom->q[i] * atom->epsilon[i]; atom->epsilon[i] = dvalue; - atom->q[i] = q_unscaled / dvalue; - atom->q_unscaled[i] = q_unscaled; + atom->q_scaled[i] = atom->q[i] / dvalue; } }