From 2809428fe20cbfd08d9e108fa88c72b27ef76c31 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 3 Mar 2024 12:27:39 -0500 Subject: [PATCH] make computes rdf and adf error out multi cutoff neighbor lists if needed --- doc/src/compute_adf.rst | 19 +++++- doc/src/compute_rdf.rst | 37 ++++++---- examples/rdf-adf/in.spce | 24 +++---- src/EXTRA-COMPUTE/compute_adf.cpp | 110 ++++++++++++++---------------- src/EXTRA-COMPUTE/compute_adf.h | 1 - src/compute_rdf.cpp | 8 ++- 6 files changed, 108 insertions(+), 91 deletions(-) diff --git a/doc/src/compute_adf.rst b/doc/src/compute_adf.rst index fc1ad1ae0a..a43a10207c 100644 --- a/doc/src/compute_adf.rst +++ b/doc/src/compute_adf.rst @@ -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 ` 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 +` 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' `. + +If you want an ADF for a larger outer cutoff, you can also use the +:doc:`rerun ` command to post-process a dump file, use :doc:`pair +style 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 diff --git a/doc/src/compute_rdf.rst b/doc/src/compute_rdf.rst index ed73800f82..85e758016e 100644 --- a/doc/src/compute_rdf.rst +++ b/doc/src/compute_rdf.rst @@ -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 ` 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' `. + +If you want an RDF for larger distances, you can also use the +:doc:`rerun ` command to post-process a dump file, use :doc:`pair +style 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:: diff --git a/examples/rdf-adf/in.spce b/examples/rdf-adf/in.spce index 9a9d99fd42..6627924adc 100644 --- a/examples/rdf-adf/in.spce +++ b/examples/rdf-adf/in.spce @@ -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 diff --git a/src/EXTRA-COMPUTE/compute_adf.cpp b/src/EXTRA-COMPUTE/compute_adf.cpp index 35ff8bfd33..20b1749fa9 100644 --- a/src/EXTRA-COMPUTE/compute_adf.cpp +++ b/src/EXTRA-COMPUTE/compute_adf.cpp @@ -34,32 +34,27 @@ #include 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); + } } /* ---------------------------------------------------------------------- */ diff --git a/src/EXTRA-COMPUTE/compute_adf.h b/src/EXTRA-COMPUTE/compute_adf.h index 5f30995aa2..f1f95d325e 100644 --- a/src/EXTRA-COMPUTE/compute_adf.h +++ b/src/EXTRA-COMPUTE/compute_adf.h @@ -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 }; diff --git a/src/compute_rdf.cpp b/src/compute_rdf.cpp index abc72ff202..89f3c91017 100644 --- a/src/compute_rdf.cpp +++ b/src/compute_rdf.cpp @@ -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); + } } /* ---------------------------------------------------------------------- */