jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

make computes rdf and adf error out multi cutoff neighbor lists if needed

This commit is contained in:
Axel Kohlmeyer 2024-03-03 12:27:39 -05:00
parent a776d8425f
commit 2809428fe2
No known key found for this signature in database
GPG Key ID: D9B44E93BF0C375A
6 changed files with 108 additions and 91 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
doc/src/compute_adf.rst
View File

@ -204,8 +204,23 @@ angles per atom satisfying the ADF criteria.
Restrictions
""""""""""""
This compute is part of the EXTRA-COMPUTE package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info.
This compute is part of the EXTRA-COMPUTE package. It is only enabled
if LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
By default, the ADF is not computed for distances longer than the
largest force cutoff, since the neighbor list creation will only contain
pairs up to that distance (plus neighbor list skin). If you use outer
cutoffs larger than that, you must use :doc:`neighbor style 'bin' or
'nsq' <neighbor>`.
If you want an ADF for a larger outer cutoff, you can also use the
:doc:`rerun <rerun>` command to post-process a dump file, use :doc:`pair
style zero <pair_zero>` and set the force cutoff to be larger in the
rerun script. Note that in the rerun context, the force cutoff is
arbitrary and with pair style zero you are not computing any forces, and
since you are not running dynamics you are not changing the model that
generated the trajectory.
The ADF is not computed for neighbors outside the force cutoff,
since processors (in parallel) don't know about atom coordinates for

37
doc/src/compute_rdf.rst
View File

@ -176,22 +176,29 @@ also numbers :math:`\ge 0.0`.
Restrictions
""""""""""""
The RDF is not computed for distances longer than the force cutoff,
since processors (in parallel) do not know about atom coordinates for
atoms further away than that distance. If you want an RDF for larger
distances, you can use the :doc:`rerun <rerun>` command to post-process
a dump file and set the cutoff for the potential to be longer in the
By default, the RDF is not computed for distances longer than the
largest force cutoff, since the neighbor list creation will only contain
pairs up to that distance (plus neighbor list skin). This distance can
be increased using the *cutoff* keyword but this keyword is only valid
with :doc:`neighbor styles 'bin' and 'nsq' <neighbor>`.
If you want an RDF for larger distances, you can also use the
:doc:`rerun <rerun>` command to post-process a dump file, use :doc:`pair
style zero <pair_zero>` and set the force cutoff to be longer in the
rerun script. Note that in the rerun context, the force cutoff is
arbitrary, since you are not running dynamics and thus are not changing
your model. The definition of :math:`g(r)` used by LAMMPS is only appropriate
for characterizing atoms that are uniformly distributed throughout the
simulation cell. In such cases, the coordination number is still
correct and meaningful. As an example, if a large simulation cell
contains only one atom of type *itypeN* and one of *jtypeN*, then :math:`g(r)`
will register an arbitrarily large spike at whatever distance they
happen to be at, and zero everywhere else.
The function :math:`\text{coord}(r)` will show a step
change from zero to one at the location of the spike in :math:`g(r)`.
arbitrary and with pair style zero you are not computing any forces, and
you are not running dynamics you are not changing the model that
generated the trajectory.
The definition of :math:`g(r)` used by LAMMPS is only appropriate for
characterizing atoms that are uniformly distributed throughout the
simulation cell. In such cases, the coordination number is still correct
and meaningful. As an example, if a large simulation cell contains only
one atom of type *itypeN* and one of *jtypeN*, then :math:`g(r)` will
register an arbitrarily large spike at whatever distance they happen to
be at, and zero everywhere else. The function :math:`\text{coord}(r)`
will show a step change from zero to one at the location of the spike in
:math:`g(r)`.
.. note::

24
examples/rdf-adf/in.spce
View File

@ -1,22 +1,22 @@
# Liquid water RDFs and ADFs (~12 O-O-O/atom, ~1 O-H...O/atom)
units real
atom_style full
units real
atom_style full
read_data data.spce
read_data data.spce
pair_style lj/cut/coul/long 12.0 12.0
pair_coeff * * 0.0 1.0
pair_coeff 1 1 0.15535 3.166
kspace_style pppm 1.0e-6
kspace_style pppm 1.0e-6
bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none
bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none
bond_coeff 1 1000.00 1.000
angle_coeff 1 100.0 109.47
bond_coeff 1 1000.00 1.000
angle_coeff 1 100.0 109.47
# need to set bond/angle inclusion to > 0.0
# so that intramolecular pairs are included in neighbor lists (required for second ADF)
@ -26,8 +26,8 @@ neighbor 2.0 bin
timestep 2.0
neigh_modify every 1 delay 2 check yes
fix 1 all shake 0.0001 20 0 b 1 a 1
fix 2 all nvt temp 300.0 300.0 100.0
fix 1 all shake 0.0001 20 0 b 1 a 1
fix 2 all nvt temp 300.0 300.0 100.0
velocity all create 300.0 6244325

110
src/EXTRA-COMPUTE/compute_adf.cpp
View File

@ -34,32 +34,27 @@
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
using MathConst::MY_PI;
using MathConst::RAD2DEG;
enum{DEGREE, RADIAN, COSINE};
enum { DEGREE, RADIAN, COSINE };
/* ----------------------------------------------------------------------
compute angular distribution functions for I, J, K atoms
---------------------------------------------------------------------- */
ComputeADF::ComputeADF(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
ilo(nullptr), ihi(nullptr), jlo(nullptr), jhi(nullptr), klo(nullptr), khi(nullptr),
hist(nullptr), histall(nullptr),
rcutinnerj(nullptr), rcutinnerk(nullptr),
rcutouterj(nullptr), rcutouterk(nullptr),
list(nullptr),
iatomcount(nullptr), iatomcountall(nullptr), iatomflag(nullptr),
maxjatom(nullptr), maxkatom(nullptr),
numjatom(nullptr), numkatom(nullptr),
neighjatom(nullptr),neighkatom(nullptr),
jatomflag(nullptr), katomflag(nullptr),
maxjkatom(nullptr), numjkatom(nullptr),
neighjkatom(nullptr), bothjkatom(nullptr), delrjkatom(nullptr)
Compute(lmp, narg, arg), ilo(nullptr), ihi(nullptr), jlo(nullptr), jhi(nullptr), klo(nullptr),
khi(nullptr), hist(nullptr), histall(nullptr), rcutinnerj(nullptr), rcutinnerk(nullptr),
rcutouterj(nullptr), rcutouterk(nullptr), list(nullptr), iatomcount(nullptr),
iatomcountall(nullptr), iatomflag(nullptr), maxjatom(nullptr), maxkatom(nullptr),
numjatom(nullptr), numkatom(nullptr), neighjatom(nullptr), neighkatom(nullptr),
jatomflag(nullptr), katomflag(nullptr), maxjkatom(nullptr), numjkatom(nullptr),
neighjkatom(nullptr), bothjkatom(nullptr), delrjkatom(nullptr)
{
int nargsperadf = 7;
if (narg < 4 ) error->all(FLERR,"Illegal compute adf command");
if (narg < 4) utils::missing_cmd_args(FLERR, "compute adf", error);
array_flag = 1;
extarray = 0;
@ -89,17 +84,16 @@ ComputeADF::ComputeADF(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[iarg+1],"degree") == 0) ordinate_style = DEGREE;
else if (strcmp(arg[iarg+1],"radian") == 0) ordinate_style = RADIAN;
else if (strcmp(arg[iarg+1],"cosine") == 0) ordinate_style = COSINE;
else error->all(FLERR,"Illegal compute adf command");
else error->all(FLERR,"Unknown compute adf ordinate flag {}",arg[iarg+1]);
iarg += 2;
} else error->all(FLERR,"Illegal compute adf command");
} else error->all(FLERR,"Unknown compute adf keyword {}", arg[iarg]);
}
// triplewise args
if (!nargtriple) ntriples = 1;
else {
if (nargtriple % nargsperadf)
error->all(FLERR,"Illegal compute adf command");
if (nargtriple % nargsperadf) error->all(FLERR,"Illegal compute adf command");
ntriples = nargtriple/nargsperadf;
}
@ -107,12 +101,9 @@ ComputeADF::ComputeADF(LAMMPS *lmp, int narg, char **arg) :
size_array_cols = 1 + 2*ntriples;
int ntypes = atom->ntypes;
memory->create(iatomflag,ntriples,ntypes+1,
"adf:iatomflag");
memory->create(jatomflag,ntriples,ntypes+1,
"adf:jatomflag");
memory->create(katomflag,ntriples,ntypes+1,
"adf:katomflag");
memory->create(iatomflag,ntriples,ntypes+1,"adf:iatomflag");
memory->create(jatomflag,ntriples,ntypes+1,"adf:jatomflag");
memory->create(katomflag,ntriples,ntypes+1,"adf:katomflag");
ilo = new int[ntriples];
ihi = new int[ntriples];
@ -134,14 +125,14 @@ ComputeADF::ComputeADF(LAMMPS *lmp, int narg, char **arg) :
klo[0] = 1; khi[0] = ntypes;
} else {
cutflag = 1;
if ((neighbor->style == Neighbor::MULTI) || (neighbor->style == Neighbor::MULTI_OLD))
error->all(FLERR, "Compute adf with custom cutoffs requires neighbor style 'bin' or 'nsq'");
iarg = 4;
for (int m = 0; m < ntriples; m++) {
utils::bounds(FLERR,arg[iarg],1,atom->ntypes,ilo[m],ihi[m],error);
utils::bounds(FLERR,arg[iarg+1],1,atom->ntypes,jlo[m],jhi[m],error);
utils::bounds(FLERR,arg[iarg+2],1,atom->ntypes,klo[m],khi[m],error);
if (ilo[m] > ihi[m] ||
jlo[m] > jhi[m] ||
klo[m] > khi[m])
if ((ilo[m] > ihi[m]) || (jlo[m] > jhi[m]) || (klo[m] > khi[m]))
error->all(FLERR,"Illegal compute adf command");
rcutinnerj[m] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
rcutouterj[m] = utils::numeric(FLERR,arg[iarg+4],false,lmp);
@ -221,8 +212,6 @@ ComputeADF::ComputeADF(LAMMPS *lmp, int narg, char **arg) :
memory->create(bothjkatom[m],maxjkatom[m],"adf:bothjkatom");
memory->create(delrjkatom[m],maxjkatom[m],4,"adf:delrjkatom");
}
rad2deg = 180.0 / MY_PI;
}
/* ---------------------------------------------------------------------- */
@ -230,47 +219,47 @@ ComputeADF::ComputeADF(LAMMPS *lmp, int narg, char **arg) :
ComputeADF::~ComputeADF()
{
memory->destroy(iatomflag);
delete [] ilo;
delete [] ihi;
delete [] jlo;
delete [] jhi;
delete [] klo;
delete [] khi;
delete [] iatomcount;
delete [] iatomcountall;
delete[] ilo;
delete[] ihi;
delete[] jlo;
delete[] jhi;
delete[] klo;
delete[] khi;
delete[] iatomcount;
delete[] iatomcountall;
memory->destroy(hist);
memory->destroy(histall);
memory->destroy(array);
memory->destroy(jatomflag);
delete [] rcutinnerj;
delete [] rcutouterj;
delete [] maxjatom;
delete [] numjatom;
delete[] rcutinnerj;
delete[] rcutouterj;
delete[] maxjatom;
delete[] numjatom;
for (int m = 0; m < ntriples; m++)
memory->destroy(neighjatom[m]);
delete [] neighjatom;
delete[] neighjatom;
memory->destroy(katomflag);
delete [] rcutinnerk;
delete [] rcutouterk;
delete [] maxkatom;
delete [] numkatom;
delete[] rcutinnerk;
delete[] rcutouterk;
delete[] maxkatom;
delete[] numkatom;
for (int m = 0; m < ntriples; m++)
memory->destroy(neighkatom[m]);
delete [] neighkatom;
delete[] neighkatom;
delete [] maxjkatom;
delete [] numjkatom;
delete[] maxjkatom;
delete[] numjkatom;
for (int m = 0; m < ntriples; m++)
memory->destroy(neighjkatom[m]);
delete [] neighjkatom;
delete[] neighjkatom;
for (int m = 0; m < ntriples; m++)
memory->destroy(bothjkatom[m]);
delete [] bothjkatom;
delete[] bothjkatom;
for (int m = 0; m < ntriples; m++)
memory->destroy(delrjkatom[m]);
delete [] delrjkatom;
delete[] delrjkatom;
}
/* ---------------------------------------------------------------------- */
@ -282,8 +271,7 @@ void ComputeADF::init()
if (!cutflag) {
if (!force->pair)
error->all(FLERR,"Compute adf requires a pair style be defined "
"or an outer cutoff specified");
error->all(FLERR,"Compute adf requires a pair style be defined or an outer cutoff specified");
rcutinnerj[0] = 0.0;
rcutinnerk[0] = 0.0;
rcutouterj[0] = force->pair->cutforce;
@ -298,7 +286,7 @@ void ComputeADF::init()
// specify mycutneigh if force cutoff too small or non-existent
if (!(force->pair) || maxouter > force->pair->cutforce) {
if (!(force->pair) || (maxouter > force->pair->cutforce)) {
double skin = neighbor->skin;
mycutneigh = maxouter + skin;
if (mycutneigh > comm->cutghostuser)
@ -310,7 +298,7 @@ void ComputeADF::init()
int x0;
if (ordinate_style == DEGREE) {
deltax = MY_PI / nbin * rad2deg;
deltax = MY_PI / nbin * RAD2DEG;
deltaxinv = nbin / MY_PI;
x0 = 0.0;
@ -337,7 +325,11 @@ void ComputeADF::init()
// than maxouter apart, just like a normal neighbor list does
auto req = neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL);
if (mycutneigh > 0.0) req->set_cutoff(mycutneigh);
if (mycutneigh > 0.0) {
if ((neighbor->style == Neighbor::MULTI) || (neighbor->style == Neighbor::MULTI_OLD))
error->all(FLERR, "Compute adf with custom cutoffs requires neighbor style 'bin' or 'nsq'");
req->set_cutoff(mycutneigh);
}
}
/* ---------------------------------------------------------------------- */

1
src/EXTRA-COMPUTE/compute_adf.h
View File

@ -59,7 +59,6 @@ class ComputeADF : public Compute {
int **bothjkatom; // 1 if atom is in both jatom and katom lists
double ***delrjkatom; // list of 4-vectors: delx, dely, delx, and 1/r
double rad2deg; // conversion factor from radians to degrees
int ordinate_style; // DEGREE, RADIAN, or COSINE
int cutflag; // 1 if at least one outer cutoff specified
};

8
src/compute_rdf.cpp
View File

@ -96,7 +96,7 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) :
jlo = new int[npairs];
jhi = new int[npairs];
if (nargpair == 0) {
if (!nargpair) {
ilo[0] = 1; ihi[0] = ntypes;
jlo[0] = 1; jhi[0] = ntypes;
} else {
@ -206,7 +206,11 @@ void ComputeRDF::init()
// than cutoff_user apart, just like a normal neighbor list does
auto req = neighbor->add_request(this, NeighConst::REQ_OCCASIONAL);
if (cutflag) req->set_cutoff(mycutneigh);
if (cutflag) {
if ((neighbor->style == Neighbor::MULTI) || (neighbor->style == Neighbor::MULTI_OLD))
error->all(FLERR, "Compute rdf with custom cutoff requires neighbor style 'bin' or 'nsq'");
req->set_cutoff(mycutneigh);
}
}
/* ---------------------------------------------------------------------- */