ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Switched to using q_scaled, keeping q as the real, unscaled charges

Browse Source
This commit is contained in:
Trung Nguyen
2022-12-09 23:34:49 -06:00
parent 9d962363a6
commit 86c2ae6dab
12 changed files with 59 additions and 69 deletions
Download Patch File Download Diff File Expand all files Collapse all files

43
src/DIELECTRIC/atom_vec_dielectric.cpp
View File

@ -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;

4
src/DIELECTRIC/atom_vec_dielectric.h
View File

@ -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

7
src/DIELECTRIC/pair_coul_cut_dielectric.cpp
View File

@ -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;

2
src/DIELECTRIC/pair_coul_long_dielectric.cpp
View File

@ -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;

7
src/DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp
View File

@ -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;
}

7
src/DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp
View File

@ -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;
}

14
src/DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp
View File

@ -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;
}

24
src/DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp
View File

@ -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;

5
src/DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp
View File

@ -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;

6
src/atom.cpp
View File

@ -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
// --------------------------------------------------------------------

2
src/atom.h
View File

@ -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
// --------------------------------------------------------------------

7
src/set.cpp
View File

@ -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;
}
}