tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'github/develop' into collected-small-changes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2025-01-22 18:01:43 -05:00
parent 19d64fd0f9 e548c656ce
commit 94eee049b4
3 changed files with 188 additions and 31 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
doc/src/pair_hybrid.rst
View File

@ -70,6 +70,12 @@ Examples
pair_coeff 1 1 lj/cut 1.0 1.0 2.5
pair_coeff 1 1 morse 1.0 1.0 1.0 2.5
variable peratom1 atom 1/(1+exp(-$k*vx^2)
variable peratom2 atom 1-v_peratom1
pair_style hybrid/scaled v_peratom1 lj/cut 2.5 v_peratom2 morse 2.5
pair_coeff 1 1 lj/cut 1.0 1.0 2.5
pair_coeff 1 1 morse 1.0 1.0 1.0 2.5
Description
"""""""""""
@ -78,7 +84,7 @@ styles enable the use of multiple pair styles in one simulation. With
the *hybrid* style, exactly one pair style is assigned to each pair of
atom types. With the *hybrid/overlay* and *hybrid/scaled* styles, one
or more pair styles can be assigned to each pair of atom types. With
the hybrid/molecular style, pair styles are assigned to either intra-
the *hybrid/molecular* style, pair styles are assigned to either intra-
or inter-molecular interactions.
The assignment of pair styles to type pairs is made via the
@ -114,16 +120,26 @@ restrictions discussed below.
If the *hybrid/scaled* style is used instead of *hybrid/overlay*,
contributions from sub-styles are weighted by their scale factors, which
may be fractional or even negative. Furthermore the scale factors may
be variables that may change during a simulation. This enables
may be fractional or even negative. Furthermore the scale factor for
each sub-style may a constant, an *equal* style variable, or an *atom*
style variable. Variable scale factors may change during the simulation.
Different sub-styles may use different scale factor styles.
In the case of a sub-style scale factor that is an *atom* style variable,
the force contribution to each atom from that sub-style is weighted
by the value of the variable for that atom, while the contribution
from that sub-style to the global potential energy is zero.
All other contributions to the per-atom energy, per-atom
virial, and global virial (if not obtained from forces)
from that sub-style are zero.
This enables
switching smoothly between two different pair styles or two different
parameter sets during a run in a similar fashion as could be done
with :doc:`fix adapt <fix_adapt>` or :doc:`fix alchemy <fix_alchemy>`.
All pair styles that will be used are listed as "sub-styles" following
the *hybrid* or *hybrid/overlay* keyword, in any order. In case of the
*hybrid/scaled* pair style, each sub-style is prefixed with a scale
factor. The scale factor is either a floating point number or an equal
factor. The scale factor is either a floating point number or an
*equal* or *atom*
style (or equivalent) variable. Each sub-style's name is followed by
its usual arguments, as illustrated in the examples above. See the doc
pages of the individual pair styles for a listing and explanation of the
@ -374,7 +390,7 @@ between all atoms of types 1,3,4 will be computed by that potential.
Pair_style hybrid allows interactions between type pairs 2-2, 1-2,
2-3, 2-4 to be specified for computation by other pair styles. You
could even add a second interaction for 1-1 to be computed by another
pair style, assuming pair_style hybrid/overlay is used.
pair style, assuming pair_style *hybrid/overlay* is used.
But you should not, as a general rule, attempt to exclude the many-body
interactions for some subset of the type pairs within the set of 1,3,4
@ -414,7 +430,7 @@ passed to the Tersoff potential, which means it would compute no
3-body interactions containing both type 1 and 2 atoms.
Here is another example to use 2 many-body potentials together in an
overlapping manner using hybrid/overlay. Imagine you have CNT (C atoms)
overlapping manner using *hybrid/overlay*. Imagine you have CNT (C atoms)
on a Si surface. You want to use Tersoff for Si/Si and Si/C
interactions, and AIREBO for C/C interactions. Si atoms are type 1; C
atoms are type 2. Something like this will work:

184
src/pair_hybrid_scaled.cpp
View File

@ -32,9 +32,14 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairHybridScaled::PairHybridScaled(LAMMPS *lmp) :
PairHybrid(lmp), fsum(nullptr), tsum(nullptr), scaleval(nullptr), scaleidx(nullptr)
PairHybrid(lmp), fsum(nullptr), tsum(nullptr), scaleval(nullptr), scaleidx(nullptr),
atomvar(nullptr), atomscale(nullptr)
{
nmaxfsum = -1;
// set comm size needed by this Pair (if atomscaleflag)
comm_forward = 1;
}
/* ---------------------------------------------------------------------- */
@ -45,6 +50,8 @@ PairHybridScaled::~PairHybridScaled()
memory->destroy(tsum);
delete[] scaleval;
delete[] scaleidx;
delete[] atomvar;
memory->destroy(atomscale);
}
/* ----------------------------------------------------------------------
@ -66,18 +73,35 @@ void PairHybridScaled::compute(int eflag, int vflag)
// update scale values from variables where needed
const int nvars = scalevars.size();
int atomscaleflag = 0;
if (nvars > 0) {
auto vals = new double[nvars];
auto vars = new int[nvars];
for (int k = 0; k < nvars; ++k) {
int m = input->variable->find(scalevars[k].c_str());
if (m < 0)
error->all(FLERR, "Variable '{}' not found when updating scale factors", scalevars[k]);
vals[k] = input->variable->compute_equal(m);
// for equal-style, compute variable, set variable index to -1
if (input->variable->equalstyle(m)) {
vals[k] = input->variable->compute_equal(m);
vars[k] = -1;
// for atom-style, store variable index, set variable to 0.0, set atomscaleflag
} else if (input->variable->atomstyle(m)) {
vals[k] = 0.0;
vars[k] = m;
atomscaleflag = 1;
} else
error->all(FLERR, "Variable '{}' has incompatible style", scalevars[k]);
}
for (int k = 0; k < nstyles; ++k) {
if (scaleidx[k] >= 0) scaleval[k] = vals[scaleidx[k]];
if (scaleidx[k] >= 0) {
scaleval[k] = vals[scaleidx[k]];
atomvar[k] = vars[scaleidx[k]];
}
}
delete[] vals;
delete[] vars;
}
// check if no_virial_fdotr_compute is set and global component of
@ -95,9 +119,11 @@ void PairHybridScaled::compute(int eflag, int vflag)
if (atom->nmax > nmaxfsum) {
memory->destroy(fsum);
if (atom->torque_flag) memory->destroy(tsum);
if (atomscaleflag) memory->destroy(atomscale);
nmaxfsum = atom->nmax;
memory->create(fsum, nmaxfsum, 3, "pair:fsum");
if (atom->torque_flag) memory->create(tsum, nmaxfsum, 3, "pair:tsum");
if (atomscaleflag) memory->create(atomscale, nmaxfsum, "pair:atomscale");
}
const int nall = atom->nlocal + atom->nghost;
auto f = atom->f;
@ -157,14 +183,34 @@ void PairHybridScaled::compute(int eflag, int vflag)
// add scaled forces to global sum
const double scale = scaleval[m];
for (i = 0; i < nall; ++i) {
fsum[i][0] += scale * f[i][0];
fsum[i][1] += scale * f[i][1];
fsum[i][2] += scale * f[i][2];
if (atom->torque_flag) {
tsum[i][0] += scale * t[i][0];
tsum[i][1] += scale * t[i][1];
tsum[i][2] += scale * t[i][2];
// if scale factor is constant or equal-style variable
if (scaleidx[m] < 0 || atomvar[m] < 0) {
for (i = 0; i < nall; ++i) {
fsum[i][0] += scale * f[i][0];
fsum[i][1] += scale * f[i][1];
fsum[i][2] += scale * f[i][2];
if (atom->torque_flag) {
tsum[i][0] += scale * t[i][0];
tsum[i][1] += scale * t[i][1];
tsum[i][2] += scale * t[i][2];
}
}
// if scale factor is atom-style variable
} else {
const int igroupall = 0;
input->variable->compute_atom(atomvar[m], igroupall, atomscale, 1, 0);
comm->forward_comm(this);
for (i = 0; i < nall; ++i) {
const double ascale = atomscale[i];
fsum[i][0] += ascale * f[i][0];
fsum[i][1] += ascale * f[i][1];
fsum[i][2] += ascale * f[i][2];
if (atom->torque_flag) {
tsum[i][0] += ascale * t[i][0];
tsum[i][1] += ascale * t[i][1];
tsum[i][2] += ascale * t[i][2];
}
}
}
@ -288,6 +334,7 @@ void PairHybridScaled::settings(int narg, char **arg)
scaleval = new double[narg];
scaleidx = new int[narg];
atomvar = new int[narg];
scalevars.reserve(narg);
// allocate each sub-style
@ -303,7 +350,8 @@ void PairHybridScaled::settings(int narg, char **arg)
while (iarg < narg - 1) {
// first process scale factor or variable
// idx < 0 indicates constant value otherwise index in variable name list
// scaleidx[k] < 0 indicates constant value, otherwise index in variable name list
// initialize atomvar[k] to -1, indicates not atom-style variable
double val = 0.0;
int idx = -1;
@ -323,6 +371,7 @@ void PairHybridScaled::settings(int narg, char **arg)
}
scaleval[nstyles] = val;
scaleidx[nstyles] = idx;
atomvar[nstyles] = -1;
++iarg;
if (utils::strmatch(arg[iarg], "^hybrid"))
@ -387,18 +436,35 @@ double PairHybridScaled::single(int i, int j, int itype, int jtype, double rsq,
// update scale values from variables where needed
const int nvars = scalevars.size();
int atomscaleflag = 0;
if (nvars > 0) {
auto vals = new double[nvars];
auto vars = new int[nvars];
for (int k = 0; k < nvars; ++k) {
int m = input->variable->find(scalevars[k].c_str());
if (m < 0)
error->all(FLERR, "Variable '{}' not found when updating scale factors", scalevars[k]);
vals[k] = input->variable->compute_equal(m);
// for equal-style, compute variable, set variable index to -1
if (input->variable->equalstyle(m)) {
vals[k] = input->variable->compute_equal(m);
vars[k] = -1;
// for atom-style, store variable index, set variable to 0.0, set atomscaleflag
} else if (input->variable->atomstyle(m)) {
vals[k] = 0.0;
vars[k] = m;
atomscaleflag = 1;
} else
error->all(FLERR, "Variable '{}' has incompatible style", scalevars[k]);
}
for (int k = 0; k < nstyles; ++k) {
if (scaleidx[k] >= 0) scaleval[k] = vals[scaleidx[k]];
if (scaleidx[k] >= 0) {
scaleval[k] = vals[scaleidx[k]];
atomvar[k] = vars[scaleidx[k]];
}
}
delete[] vals;
delete[] vars;
}
double fone;
@ -417,7 +483,18 @@ double PairHybridScaled::single(int i, int j, int itype, int jtype, double rsq,
double scale = scaleval[map[itype][jtype][m]];
esum += scale * pstyle->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fone);
fforce += scale * fone;
// if scale factor is constant or equal-style variable
if (scaleidx[m] < 0 || atomvar[m] < 0) {
fforce += scale * fone;
// if scale factor is atom-style variable, average i and j
} else {
const int igroupall = 0;
input->variable->compute_atom(atomvar[m], igroupall, atomscale, 1, 0);
comm->forward_comm(this);
const double ascale = 0.5 * (atomscale[i] + atomscale[j]);
fforce += ascale * fone;
}
}
}
@ -440,18 +517,35 @@ void PairHybridScaled::born_matrix(int i, int j, int itype, int jtype, double rs
// update scale values from variables where needed
const int nvars = scalevars.size();
int atomscaleflag = 0;
if (nvars > 0) {
double *vals = new double[nvars];
auto vals = new double[nvars];
auto vars = new int[nvars];
for (int k = 0; k < nvars; ++k) {
int m = input->variable->find(scalevars[k].c_str());
if (m < 0)
error->all(FLERR, "Variable '{}' not found when updating scale factors", scalevars[k]);
vals[k] = input->variable->compute_equal(m);
// for equal-style, compute variable, set variable index to -1
if (input->variable->equalstyle(m)) {
vals[k] = input->variable->compute_equal(m);
vars[k] = -1;
// for atom-style, store variable index, set variable to 0.0, set atomscaleflag
} else if (input->variable->atomstyle(m)) {
vals[k] = 0.0;
vars[k] = m;
atomscaleflag = 1;
} else
error->all(FLERR, "Variable '{}' has incompatible style", scalevars[k]);
}
for (int k = 0; k < nstyles; ++k) {
if (scaleidx[k] >= 0) scaleval[k] = vals[scaleidx[k]];
if (scaleidx[k] >= 0) {
scaleval[k] = vals[scaleidx[k]];
atomvar[k] = vars[scaleidx[k]];
}
}
delete[] vals;
delete[] vars;
}
double du, du2, scale;
@ -460,18 +554,30 @@ void PairHybridScaled::born_matrix(int i, int j, int itype, int jtype, double rs
for (int m = 0; m < nmap[itype][jtype]; m++) {
auto pstyle = styles[map[itype][jtype][m]];
if (rsq < pstyle->cutsq[itype][jtype]) {
if (pstyle->born_matrix_enable == 0)
error->one(FLERR, "Pair hybrid sub-style does not support born_matrix call");
if (pstyle->single_enable == 0)
error->one(FLERR, "Pair hybrid sub-style does not support single call");
if ((special_lj[map[itype][jtype][m]] != nullptr) ||
(special_coul[map[itype][jtype][m]] != nullptr))
error->one(FLERR, "Pair hybrid born_matrix() does not support per sub-style special_bond");
error->one(FLERR, "Pair hybrid single() does not support per sub-style special_bond");
du = du2 = 0.0;
scale = scaleval[map[itype][jtype][m]];
double scale = scaleval[map[itype][jtype][m]];
pstyle->born_matrix(i, j, itype, jtype, rsq, factor_coul, factor_lj, du, du2);
dupair += scale * du;
du2pair += scale * du2;
// if scale factor is constant or equal-style variable
if (scaleidx[m] < 0 || atomvar[m] < 0) {
dupair += scale * du;
du2pair += scale * du2;
// if scale factor is atom-style variable, average i and j
} else {
const int igroupall = 0;
input->variable->compute_atom(atomvar[m], igroupall, atomscale, 1, 0);
comm->forward_comm(this);
const double ascale = 0.5 * (atomscale[i] + atomscale[j]);
dupair += ascale * du;
du2pair += ascale * du2;
}
}
}
}
@ -574,6 +680,7 @@ void PairHybridScaled::write_restart(FILE *fp)
fwrite(scaleval, sizeof(double), nstyles, fp);
fwrite(scaleidx, sizeof(int), nstyles, fp);
fwrite(atomvar, sizeof(int), nstyles, fp);
int n = scalevars.size();
fwrite(&n, sizeof(int), 1, fp);
@ -594,17 +701,21 @@ void PairHybridScaled::read_restart(FILE *fp)
delete[] scaleval;
delete[] scaleidx;
delete[] atomvar;
scalevars.clear();
scaleval = new double[nstyles];
scaleidx = new int[nstyles];
atomvar = new int[nstyles];
int n, me = comm->me;
if (me == 0) {
utils::sfread(FLERR, scaleval, sizeof(double), nstyles, fp, nullptr, error);
utils::sfread(FLERR, scaleidx, sizeof(int), nstyles, fp, nullptr, error);
utils::sfread(FLERR, atomvar, sizeof(int), nstyles, fp, nullptr, error);
}
MPI_Bcast(scaleval, nstyles, MPI_DOUBLE, 0, world);
MPI_Bcast(scaleidx, nstyles, MPI_INT, 0, world);
MPI_Bcast(atomvar, nstyles, MPI_INT, 0, world);
char *tmp;
if (me == 0) utils::sfread(FLERR, &n, sizeof(int), 1, fp, nullptr, error);
@ -667,3 +778,28 @@ void PairHybridScaled::copy_svector(int itype, int jtype)
}
}
}
/* ---------------------------------------------------------------------- */
int PairHybridScaled::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = atomscale[j];
}
return m;
}
/* ---------------------------------------------------------------------- */
void PairHybridScaled::unpack_forward_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) atomscale[i] = buf[m++];
}

5
src/pair_hybrid_scaled.h
View File

@ -44,12 +44,17 @@ class PairHybridScaled : public PairHybrid {
void init_svector() override;
void copy_svector(int, int) override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
protected:
double **fsum, **tsum;
double *scaleval;
int *scaleidx;
std::vector<std::string> scalevars;
int nmaxfsum;
int *atomvar; // indices of atom-style variables
double *atomscale; // vector of atom-style variable values
};
} // namespace LAMMPS_NS