From c660dcefd2c6bcf7d6811c158a65b68a85b000d6 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 8 Apr 2021 17:42:32 -0400 Subject: [PATCH 01/14] add minimal framework for pair style hybrid/scaled --- src/pair_hybrid_scaled.cpp | 201 +++++++++++++++++++++++++++++++++++++ src/pair_hybrid_scaled.h | 63 ++++++++++++ 2 files changed, 264 insertions(+) create mode 100644 src/pair_hybrid_scaled.cpp create mode 100644 src/pair_hybrid_scaled.h diff --git a/src/pair_hybrid_scaled.cpp b/src/pair_hybrid_scaled.cpp new file mode 100644 index 0000000000..9f3a5d44be --- /dev/null +++ b/src/pair_hybrid_scaled.cpp @@ -0,0 +1,201 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://lammps.sandia.gov/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "pair_hybrid_scaled.h" + +#include "atom.h" +#include "error.h" +#include "memory.h" + +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairHybridScaled::PairHybridScaled(LAMMPS *lmp) : + PairHybrid(lmp), fsum(nullptr), + scaleval(nullptr), scaleidx(nullptr), scalevar(nullptr) +{ + nscalevar = 0; + nmaxfsum = -1; +} + +/* ---------------------------------------------------------------------- */ + +PairHybridScaled::~PairHybridScaled() +{ + for (int i=0; i < nscalevar; ++i) + delete[] scalevar[i]; + delete[] scalevar; + + if (nmaxfsum > 0) + memory->destroy(fsum); + + if (allocated) { + memory->destroy(scaleval); + memory->destroy(scaleidx); + } +} +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairHybridScaled::allocate() +{ + PairHybrid::allocate(); + + const int n = atom->ntypes; + + memory->create(scaleval,n+1,n+1,"pair:scaleval"); + memory->create(scaleidx,n+1,n+1,"pair:scaleidx"); + for (int i = 1; i <= n; i++) { + for (int j = i; j <= n; j++) { + scaleval[i][j] = 0.0; + scaleidx[i][j] = -1; + } + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairHybridScaled::coeff(int narg, char **arg) +{ + if (narg < 3) error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); + utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + + // 3rd arg = scale factor for sub-style, must be either + // a constant or equal stye or compatible variable + + double scale = utils::numeric(FLERR,arg[2],false,lmp); + + // 4th arg = pair sub-style name + // 5th arg = pair sub-style index if name used multiple times + // + // allow for "none" as valid sub-style name + + int multflag; + int m; + + for (m = 0; m < nstyles; m++) { + multflag = 0; + if (strcmp(arg[3],keywords[m]) == 0) { + if (multiple[m]) { + multflag = 1; + if (narg < 5) error->all(FLERR,"Incorrect args for pair coefficients"); + if (multiple[m] == utils::inumeric(FLERR,arg[4],false,lmp)) break; + else continue; + } else break; + } + } + + int none = 0; + if (m == nstyles) { + if (strcmp(arg[3],"none") == 0) none = 1; + else error->all(FLERR,"Pair coeff for hybrid has invalid style"); + } + + // move 1st/2nd args to 3rd/4th args + // if multflag: move 1st/2nd args to 4th/5th args + // just copy ptrs, since arg[] points into original input line + + arg[3+multflag] = arg[1]; + arg[2+multflag] = arg[0]; + + // invoke sub-style coeff() starting with 1st remaining arg + + if (!none) styles[m]->coeff(narg-2-multflag,arg+2+multflag); + + // set setflag and which type pairs map to which sub-style + // if sub-style is none: set hybrid subflag, wipe out map + // else: set hybrid setflag & map only if substyle setflag is set + // if sub-style is new for type pair, add as multiple mapping + // if sub-style exists for type pair, don't add, just update coeffs + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + if (none) { + setflag[i][j] = 1; + nmap[i][j] = 0; + count++; + } else if (styles[m]->setflag[i][j]) { + int k; + for (k = 0; k < nmap[i][j]; k++) + if (map[i][j][k] == m) break; + if (k == nmap[i][j]) map[i][j][nmap[i][j]++] = m; + setflag[i][j] = 1; + scaleval[i][j] = scale; + scaleidx[i][j] = -1; + count++; + } + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + + +/* ---------------------------------------------------------------------- + we need to handle Pair::svector special for hybrid/scaled +------------------------------------------------------------------------- */ + +void PairHybridScaled::init_svector() +{ + // single_extra = list all sub-style single_extra + // allocate svector + + single_extra = 0; + for (int m = 0; m < nstyles; m++) + single_extra += styles[m]->single_extra; + + if (single_extra) { + delete [] svector; + svector = new double[single_extra]; + } +} + +/* ---------------------------------------------------------------------- + we need to handle Pair::svector special for hybrid/scaled +------------------------------------------------------------------------- */ + +void PairHybridScaled::copy_svector(int itype, int jtype) +{ + int n=0; + Pair *this_style; + + // fill svector array. + // copy data from active styles and use 0.0 for inactive ones + for (int m = 0; m < nstyles; m++) { + for (int k = 0; k < nmap[itype][jtype]; ++k) { + if (m == map[itype][jtype][k]) { + this_style = styles[m]; + } else { + this_style = nullptr; + } + } + for (int l = 0; l < styles[m]->single_extra; ++l) { + if (this_style) { + svector[n++] = this_style->svector[l]; + } else { + svector[n++] = 0.0; + } + } + } +} diff --git a/src/pair_hybrid_scaled.h b/src/pair_hybrid_scaled.h new file mode 100644 index 0000000000..d0d2d8838a --- /dev/null +++ b/src/pair_hybrid_scaled.h @@ -0,0 +1,63 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(hybrid/scaled,PairHybridScaled) + +#else + +#ifndef LMP_PAIR_HYBRID_SCALED_H +#define LMP_PAIR_HYBRID_SCALED_H + +#include "pair_hybrid.h" + +namespace LAMMPS_NS { + +class PairHybridScaled : public PairHybrid { + public: + PairHybridScaled(class LAMMPS *); + virtual ~PairHybridScaled(); + void coeff(int, char **); + //void compute(int, int); + + void init_svector(); + void copy_svector(int,int); + + protected: + double **fsum; + double **scaleval; + int **scaleidx; + char **scalevar; + int nscalevar; + int nmaxfsum; + + void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair coeff for hybrid has invalid style + +Style in pair coeff must have been listed in pair_style command. + +*/ From 8ed6d80b851a4dbe02aa61bcd50ca12e5d3c15b1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 8 Apr 2021 21:02:28 -0400 Subject: [PATCH 02/14] first working version, forces only, no restart --- doc/src/pair_hybrid.rst | 98 ++++++---- src/pair_hybrid.h | 2 +- src/pair_hybrid_scaled.cpp | 358 ++++++++++++++++++++++++++++++++----- src/pair_hybrid_scaled.h | 17 +- 4 files changed, 386 insertions(+), 89 deletions(-) diff --git a/doc/src/pair_hybrid.rst b/doc/src/pair_hybrid.rst index 558c326175..c851423029 100644 --- a/doc/src/pair_hybrid.rst +++ b/doc/src/pair_hybrid.rst @@ -2,6 +2,7 @@ .. index:: pair_style hybrid/kk .. index:: pair_style hybrid/overlay .. index:: pair_style hybrid/overlay/kk +.. index:: pair_style hybrid/scaled pair_style hybrid command ========================= @@ -13,6 +14,8 @@ pair_style hybrid/overlay command Accelerator Variants: *hybrid/overlay/kk* +pair_style hybrid/scale command + Syntax """""" @@ -20,8 +23,10 @@ Syntax pair_style hybrid style1 args style2 args ... pair_style hybrid/overlay style1 args style2 args ... + pair_style hybrid/scaled factor1 style1 args factor2 style 2 args ... * style1,style2 = list of one or more pair styles and their arguments +* factor1,factor2 = scale factors for pair styles, may be a variable Examples """""""" @@ -37,15 +42,24 @@ Examples pair_coeff * * lj/cut 1.0 1.0 pair_coeff * * coul/long + variable one equal ramp(1.0,0.0) + variable two equal 1.0-v_one + pair_style hybrid/scaled v_one lj/cut 2.5 v_two 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 """"""""""" -The *hybrid* and *hybrid/overlay* 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* style, one or more pair styles can be assigned to -each pair of atom types. The assignment of pair styles to type pairs -is made via the :doc:`pair_coeff ` command. +The *hybrid*, *hybrid/overlay*, and *hybrid/scaled* 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. The assignment of pair styles +to type pairs is made via the :doc:`pair_coeff ` command. +The *hybrid/scaled* style differs from the *hybrid/overlay* style by +requiring a factor for each pair style that is used to scale all +forces and energies computed by the pair style. Here are two examples of hybrid simulations. The *hybrid* style could be used for a simulation of a metal droplet on a LJ surface. The @@ -61,12 +75,19 @@ it would be more efficient to use the single combined potential, but in general any combination of pair potentials can be used together in to produce an interaction that is not encoded in any single pair_style file, e.g. adding Coulombic forces between granular particles. +The *hybrid/scaled* style enables more complex combinations of pair +styles than a simple sum as *hybrid/overlay* does. Furthermore, since +the scale factors can be variables, they can change during a simulation +which would allow to smoothly switch between two different pair styles +or two different parameter sets. All pair styles that will be used are listed as "sub-styles" following -the *hybrid* or *hybrid/overlay* keyword, in any order. Each -sub-style's name is followed by its usual arguments, as illustrated in -the example above. See the doc pages of individual pair styles for a -listing and explanation of the appropriate arguments. +the *hybrid* or *hybrid/overlay* keyword, in any order. In case of the +*hybrid/scaled* pair style each sub-style is prefixed with its scale +factor. The scale factor may be an equal style (or equivalent) +variable. Each sub-style's name is followed by its usual arguments, as +illustrated in the example above. See the doc pages of individual pair +styles for a listing and explanation of the appropriate arguments. Note that an individual pair style can be used multiple times as a sub-style. For efficiency this should only be done if your model @@ -143,16 +164,16 @@ one sub-style. Just as with a simulation using a single pair style, if you specify the same atom type pair in a second pair_coeff command, the previous assignment will be overwritten. -For the *hybrid/overlay* style, each atom type pair I,J can be -assigned to one or more sub-styles. If you specify the same atom type -pair in a second pair_coeff command with a new sub-style, then the -second sub-style is added to the list of potentials that will be -calculated for two interacting atoms of those types. If you specify -the same atom type pair in a second pair_coeff command with a -sub-style that has already been defined for that pair of atoms, then -the new pair coefficients simply override the previous ones, as in the -normal usage of the pair_coeff command. E.g. these two sets of -commands are the same: +For the *hybrid/overlay* and *hybrid/scaled* style, each atom type pair +I,J can be assigned to one or more sub-styles. If you specify the same +atom type pair in a second pair_coeff command with a new sub-style, then +the second sub-style is added to the list of potentials that will be +calculated for two interacting atoms of those types. If you specify the +same atom type pair in a second pair_coeff command with a sub-style that +has already been defined for that pair of atoms, then the new pair +coefficients simply override the previous ones, as in the normal usage +of the pair_coeff command. E.g. these two sets of commands are the +same: .. code-block:: LAMMPS @@ -170,19 +191,20 @@ data file or restart files read by the :doc:`read_data ` or :doc:`read_restart ` commands, or by mixing as described below. -For both the *hybrid* and *hybrid/overlay* styles, every atom type -pair I,J (where I <= J) must be assigned to at least one sub-style via -the :doc:`pair_coeff ` command as in the examples above, or -in the data file read by the :doc:`read_data `, or by mixing -as described below. +For all of the *hybrid*, *hybrid/overlay*, and *hybrid/scaled* styles, +every atom type pair I,J (where I <= J) must be assigned to at least one +sub-style via the :doc:`pair_coeff ` command as in the +examples above, or in the data file read by the :doc:`read_data +`, or by mixing as described below. If you want there to be no interactions between a particular pair of atom types, you have 3 choices. You can assign the type pair to some sub-style and use the :doc:`neigh_modify exclude type ` command. You can assign it to some sub-style and set the coefficients -so that there is effectively no interaction (e.g. epsilon = 0.0 in a -LJ potential). Or, for *hybrid* and *hybrid/overlay* simulations, you -can use this form of the pair_coeff command in your input script: +so that there is effectively no interaction (e.g. epsilon = 0.0 in a LJ +potential). Or, for *hybrid*, *hybrid/overlay*, or *hybrid/scaled* +simulations, you can use this form of the pair_coeff command in your +input script: .. code-block:: LAMMPS @@ -260,7 +282,10 @@ the specific syntax, requirements and restrictions. ---------- The potential energy contribution to the overall system due to an -individual sub-style can be accessed and output via the :doc:`compute pair ` command. +individual sub-style can be accessed and output via the :doc:`compute +pair ` command. Note that in the case of pair style +*hybrid/scaled* this is the **unscaled** potential energy of the +selected sub-style. ---------- @@ -388,12 +413,12 @@ The hybrid pair styles supports the :doc:`pair_modify ` shift, table, and tail options for an I,J pair interaction, if the associated sub-style supports it. -For the hybrid pair styles, the list of sub-styles and their -respective settings are written to :doc:`binary restart files `, so a :doc:`pair_style ` command does -not need to specified in an input script that reads a restart file. -However, the coefficient information is not stored in the restart -file. Thus, pair_coeff commands need to be re-specified in the -restart input script. +For the hybrid pair styles, the list of sub-styles and their respective +settings are written to :doc:`binary restart files `, so a +:doc:`pair_style ` command does not need to specified in an +input script that reads a restart file. However, the coefficient +information is not stored in the restart file. Thus, pair_coeff +commands need to be re-specified in the restart input script. These pair styles support the use of the *inner*\ , *middle*\ , and *outer* keywords of the :doc:`run_style respa ` command, if @@ -409,6 +434,9 @@ e.g. *lj/cut/coul/long* or *buck/coul/long*\ . You must insure that the short-range Coulombic cutoff used by each of these long pair styles is the same or else LAMMPS will generate an error. +Pair style *hybrid/scaled* currently only works for non-accelerated +pair styles. + Related commands """""""""""""""" diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h index 7717d1fd51..3bde620514 100644 --- a/src/pair_hybrid.h +++ b/src/pair_hybrid.h @@ -37,7 +37,7 @@ class PairHybrid : public Pair { PairHybrid(class LAMMPS *); virtual ~PairHybrid(); virtual void compute(int, int); - void settings(int, char **); + virtual void settings(int, char **); virtual void coeff(int, char **); void init_style(); double init_one(int, int); diff --git a/src/pair_hybrid_scaled.cpp b/src/pair_hybrid_scaled.cpp index 9f3a5d44be..7e1d634627 100644 --- a/src/pair_hybrid_scaled.cpp +++ b/src/pair_hybrid_scaled.cpp @@ -15,7 +15,12 @@ #include "atom.h" #include "error.h" +#include "force.h" +#include "input.h" #include "memory.h" +#include "respa.h" +#include "update.h" +#include "variable.h" #include @@ -23,11 +28,9 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairHybridScaled::PairHybridScaled(LAMMPS *lmp) : - PairHybrid(lmp), fsum(nullptr), - scaleval(nullptr), scaleidx(nullptr), scalevar(nullptr) +PairHybridScaled::PairHybridScaled(LAMMPS *lmp) + : PairHybrid(lmp), fsum(nullptr), scaleval(nullptr), scaleidx(nullptr) { - nscalevar = 0; nmaxfsum = -1; } @@ -35,36 +38,307 @@ PairHybridScaled::PairHybridScaled(LAMMPS *lmp) : PairHybridScaled::~PairHybridScaled() { - for (int i=0; i < nscalevar; ++i) - delete[] scalevar[i]; - delete[] scalevar; - - if (nmaxfsum > 0) - memory->destroy(fsum); - - if (allocated) { - memory->destroy(scaleval); - memory->destroy(scaleidx); - } + memory->destroy(fsum); + memory->destroy(scaleval); } + /* ---------------------------------------------------------------------- - allocate all arrays + call each sub-style's compute() or compute_outer() function + accumulate sub-style global/peratom energy/virial in hybrid + for global vflag = VIRIAL_PAIR: + each sub-style computes own virial[6] + sum sub-style virial[6] to hybrid's virial[6] + for global vflag = VIRIAL_FDOTR: + call sub-style with adjusted vflag to prevent it calling + virial_fdotr_compute() + hybrid calls virial_fdotr_compute() on final accumulated f ------------------------------------------------------------------------- */ -void PairHybridScaled::allocate() +void PairHybridScaled::compute(int eflag, int vflag) { - PairHybrid::allocate(); + int i,j,m,n; - const int n = atom->ntypes; + // update scale values from variables where needed - memory->create(scaleval,n+1,n+1,"pair:scaleval"); - memory->create(scaleidx,n+1,n+1,"pair:scaleidx"); - for (int i = 1; i <= n; i++) { - for (int j = i; j <= n; j++) { - scaleval[i][j] = 0.0; - scaleidx[i][j] = -1; + const int nvars = scalevars.size(); + if (nvars > 0) { + double *vals = new double[nvars]; + for (i = 0; i < nvars; ++i) { + j = input->variable->find(scalevars[i].c_str()); + vals[i] = input->variable->compute_equal(j); + } + for (i = 0; i < nstyles; ++i) { + if (scaleidx[i] >= 0) + scaleval[i] = vals[scaleidx[i]]; + } + delete[] vals; + } + + // check if no_virial_fdotr_compute is set and global component of + // incoming vflag = VIRIAL_FDOTR + // if so, reset vflag as if global component were VIRIAL_PAIR + // necessary since one or more sub-styles cannot compute virial as F dot r + + if (no_virial_fdotr_compute && (vflag & VIRIAL_FDOTR)) + vflag = VIRIAL_PAIR | (vflag & ~VIRIAL_FDOTR); + + ev_init(eflag,vflag); + + // grow fsum array if needed, and copy existing forces (usually 0.0) to it. + + if (atom->nmax > nmaxfsum) { + memory->destroy(fsum); + nmaxfsum = atom->nmax; + memory->create(fsum,nmaxfsum,3,"pair:fsum"); + } + const int nall = atom->nlocal + atom->nghost; + auto f = atom->f; + for (i = 0; i < nall; ++i) { + fsum[i][0] = f[i][0]; + fsum[i][1] = f[i][1]; + fsum[i][2] = f[i][2]; + } + + // check if global component of incoming vflag = VIRIAL_FDOTR + // if so, reset vflag passed to substyle so VIRIAL_FDOTR is turned off + // necessary so substyle will not invoke virial_fdotr_compute() + + int vflag_substyle; + if (vflag & VIRIAL_FDOTR) vflag_substyle = vflag & ~VIRIAL_FDOTR; + else vflag_substyle = vflag; + + double *saved_special = save_special(); + + // check if we are running with r-RESPA using the hybrid keyword + + Respa *respa = nullptr; + respaflag = 0; + if (utils::strmatch(update->integrate_style,"^respa")) { + respa = (Respa *) update->integrate; + if (respa->nhybrid_styles > 0) respaflag = 1; + } + + for (m = 0; m < nstyles; m++) { + + // clear forces + + memset(&f[0][0],0,nall*3*sizeof(double)); + + set_special(m); + + if (!respaflag || (respaflag && respa->hybrid_compute[m])) { + + // invoke compute() unless compute flag is turned off or + // outerflag is set and sub-style has a compute_outer() method + + if (styles[m]->compute_flag == 0) continue; + if (outerflag && styles[m]->respa_enable) + styles[m]->compute_outer(eflag,vflag_substyle); + else styles[m]->compute(eflag,vflag_substyle); + } + + // 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]; + } + + restore_special(saved_special); + + // jump to next sub-style if r-RESPA does not want global accumulated data + + if (respaflag && !respa->tally_global) continue; + + if (eflag_global) { + eng_vdwl += scale * styles[m]->eng_vdwl; + eng_coul += scale * styles[m]->eng_coul; + } + if (vflag_global) { + for (n = 0; n < 6; n++) virial[n] += scale * styles[m]->virial[n]; + } + if (eflag_atom) { + n = atom->nlocal; + if (force->newton_pair) n += atom->nghost; + double *eatom_substyle = styles[m]->eatom; + for (i = 0; i < n; i++) eatom[i] += scale * eatom_substyle[i]; + } + if (vflag_atom) { + n = atom->nlocal; + if (force->newton_pair) n += atom->nghost; + double **vatom_substyle = styles[m]->vatom; + for (i = 0; i < n; i++) + for (j = 0; j < 6; j++) + vatom[i][j] += scale * vatom_substyle[i][j]; + } + + // substyles may be CENTROID_SAME or CENTROID_AVAIL + + if (cvflag_atom) { + n = atom->nlocal; + if (force->newton_pair) n += atom->nghost; + if (styles[m]->centroidstressflag == CENTROID_AVAIL) { + double **cvatom_substyle = styles[m]->cvatom; + for (i = 0; i < n; i++) + for (j = 0; j < 9; j++) + cvatom[i][j] += scale * cvatom_substyle[i][j]; + } else { + double **vatom_substyle = styles[m]->vatom; + for (i = 0; i < n; i++) { + for (j = 0; j < 6; j++) { + cvatom[i][j] += scale * vatom_substyle[i][j]; + } + for (j = 6; j < 9; j++) { + cvatom[i][j] += scale * vatom_substyle[i][j-3]; + } + } + } } } + + // copy accumulated forces to original force array + + for (i = 0; i < nall; ++i) { + f[i][0] = fsum[i][0]; + f[i][1] = fsum[i][1]; + f[i][2] = fsum[i][2]; + } + delete [] saved_special; + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + create one pair style for each arg in list +------------------------------------------------------------------------- */ + +void PairHybridScaled::settings(int narg, char **arg) +{ + if (narg < 1) error->all(FLERR,"Illegal pair_style command"); + if (lmp->kokkos && !utils::strmatch(force->pair_style,"^hybrid.*/kk$")) + error->all(FLERR,fmt::format("Must use pair_style {}/kk with Kokkos", + force->pair_style)); + + // delete old lists, since cannot just change settings + + if (nstyles > 0) { + for (int m = 0; m < nstyles; m++) { + delete styles[m]; + delete [] keywords[m]; + if (special_lj[m]) delete [] special_lj[m]; + if (special_coul[m]) delete [] special_coul[m]; + } + delete[] styles; + delete[] keywords; + delete[] multiple; + delete[] special_lj; + delete[] special_coul; + delete[] compute_tally; + delete[] scaleval; + delete[] scaleidx; + scalevars.clear(); + } + + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + memory->destroy(cutghost); + memory->destroy(nmap); + memory->destroy(map); + } + allocated = 0; + + // allocate list of sub-styles as big as possibly needed if no extra args + + styles = new Pair*[narg]; + keywords = new char*[narg]; + multiple = new int[narg]; + + special_lj = new double*[narg]; + special_coul = new double*[narg]; + compute_tally = new int[narg]; + + scaleval = new double[narg]; + scaleidx = new int[narg]; + scalevars.reserve(narg); + + // allocate each sub-style + // allocate uses suffix, but don't store suffix version in keywords, + // else syntax in coeff() will not match + // call settings() with set of args that are not pair style names + // use force->pair_map to determine which args these are + + int iarg,jarg,dummy; + + iarg = 0; + nstyles = 0; + while (iarg < narg-1) { + + // first process scale factor or variable + // idx < 0 indicates constant value otherwise index in variable name list + + double val = 0.0; + int idx = -1; + if (utils::strmatch(arg[iarg],"^v_")) { + for (std::size_t i=0; i < scalevars.size(); ++i) { + if (scalevars[i] == arg[iarg]+2) { + idx = i; + break; + } + } + if (idx < 0) { + idx = scalevars.size(); + scalevars.push_back(arg[iarg]+2); + } + } else { + val = utils::numeric(FLERR,arg[iarg],false,lmp); + } + scaleval[nstyles] = val; + scaleidx[nstyles] = idx; + ++iarg; + + if (utils::strmatch(arg[iarg],"^hybrid")) + error->all(FLERR,"Pair style hybrid cannot have hybrid as an argument"); + if (strcmp(arg[iarg],"none") == 0) + error->all(FLERR,"Pair style hybrid cannot have none as an argument"); + + styles[nstyles] = force->new_pair(arg[iarg],1,dummy); + force->store_style(keywords[nstyles],arg[iarg],0); + special_lj[nstyles] = special_coul[nstyles] = nullptr; + compute_tally[nstyles] = 1; + + // determine list of arguments for pair style settings + // by looking for the next known pair style name. + + jarg = iarg + 1; + while ((jarg < narg) + && !force->pair_map->count(arg[jarg]) + && !lmp->match_style("pair",arg[jarg])) jarg++; + + // decrement to account for scale factor except when last argument + + if (jarg < narg) --jarg; + + styles[nstyles]->settings(jarg-iarg-1,arg+iarg+1); + iarg = jarg; + nstyles++; + } + + // multiple[i] = 1 to M if sub-style used multiple times, else 0 + + for (int i = 0; i < nstyles; i++) { + int count = 0; + for (int j = 0; j < nstyles; j++) { + if (strcmp(keywords[j],keywords[i]) == 0) count++; + if (j == i) multiple[i] = count; + } + if (count == 1) multiple[i] = 0; + } + + // set pair flags from sub-style flags + + flags(); } /* ---------------------------------------------------------------------- @@ -80,14 +354,8 @@ void PairHybridScaled::coeff(int narg, char **arg) utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); - // 3rd arg = scale factor for sub-style, must be either - // a constant or equal stye or compatible variable - - double scale = utils::numeric(FLERR,arg[2],false,lmp); - - // 4th arg = pair sub-style name - // 5th arg = pair sub-style index if name used multiple times - // + // 3rd arg = pair sub-style name + // 4th arg = pair sub-style index if name used multiple times // allow for "none" as valid sub-style name int multflag; @@ -95,11 +363,14 @@ void PairHybridScaled::coeff(int narg, char **arg) for (m = 0; m < nstyles; m++) { multflag = 0; - if (strcmp(arg[3],keywords[m]) == 0) { + if (strcmp(arg[2],keywords[m]) == 0) { if (multiple[m]) { multflag = 1; - if (narg < 5) error->all(FLERR,"Incorrect args for pair coefficients"); - if (multiple[m] == utils::inumeric(FLERR,arg[4],false,lmp)) break; + if (narg < 4) error->all(FLERR,"Incorrect args for pair coefficients"); + if (!isdigit(arg[3][0])) + error->all(FLERR,"Incorrect args for pair coefficients"); + int index = utils::inumeric(FLERR,arg[3],false,lmp); + if (index == multiple[m]) break; else continue; } else break; } @@ -107,20 +378,20 @@ void PairHybridScaled::coeff(int narg, char **arg) int none = 0; if (m == nstyles) { - if (strcmp(arg[3],"none") == 0) none = 1; + if (strcmp(arg[2],"none") == 0) none = 1; else error->all(FLERR,"Pair coeff for hybrid has invalid style"); } - // move 1st/2nd args to 3rd/4th args - // if multflag: move 1st/2nd args to 4th/5th args + // move 1st/2nd args to 2nd/3rd args + // if multflag: move 1st/2nd args to 3rd/4th args // just copy ptrs, since arg[] points into original input line - arg[3+multflag] = arg[1]; - arg[2+multflag] = arg[0]; + arg[2+multflag] = arg[1]; + arg[1+multflag] = arg[0]; // invoke sub-style coeff() starting with 1st remaining arg - if (!none) styles[m]->coeff(narg-2-multflag,arg+2+multflag); + if (!none) styles[m]->coeff(narg-1-multflag,&arg[1+multflag]); // set setflag and which type pairs map to which sub-style // if sub-style is none: set hybrid subflag, wipe out map @@ -141,8 +412,6 @@ void PairHybridScaled::coeff(int narg, char **arg) if (map[i][j][k] == m) break; if (k == nmap[i][j]) map[i][j][nmap[i][j]++] = m; setflag[i][j] = 1; - scaleval[i][j] = scale; - scaleidx[i][j] = -1; count++; } } @@ -151,7 +420,6 @@ void PairHybridScaled::coeff(int narg, char **arg) if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } - /* ---------------------------------------------------------------------- we need to handle Pair::svector special for hybrid/scaled ------------------------------------------------------------------------- */ diff --git a/src/pair_hybrid_scaled.h b/src/pair_hybrid_scaled.h index d0d2d8838a..a8bd5517ec 100644 --- a/src/pair_hybrid_scaled.h +++ b/src/pair_hybrid_scaled.h @@ -22,27 +22,28 @@ PairStyle(hybrid/scaled,PairHybridScaled) #include "pair_hybrid.h" +#include +#include + namespace LAMMPS_NS { class PairHybridScaled : public PairHybrid { public: PairHybridScaled(class LAMMPS *); virtual ~PairHybridScaled(); - void coeff(int, char **); - //void compute(int, int); + virtual void compute(int, int); + virtual void settings(int, char**); + virtual void coeff(int, char **); void init_svector(); void copy_svector(int,int); protected: double **fsum; - double **scaleval; - int **scaleidx; - char **scalevar; - int nscalevar; + double *scaleval; + int *scaleidx; + std::vector scalevars; int nmaxfsum; - - void allocate(); }; } From ae27d3bf4cdea5a86b8f3bfe1681c2cb3d990fc6 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 9 Apr 2021 00:31:43 -0400 Subject: [PATCH 03/14] add support for single() and read/write_restart() --- src/pair_hybrid.h | 6 +- src/pair_hybrid_scaled.cpp | 126 ++++++++++++++++++++++++++++++++++++- src/pair_hybrid_scaled.h | 4 ++ 3 files changed, 130 insertions(+), 6 deletions(-) diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h index 3bde620514..ca79163fc2 100644 --- a/src/pair_hybrid.h +++ b/src/pair_hybrid.h @@ -42,9 +42,9 @@ class PairHybrid : public Pair { void init_style(); double init_one(int, int); void setup(); - void write_restart(FILE *); - void read_restart(FILE *); - double single(int, int, int, int, double, double, double, double &); + virtual void write_restart(FILE *); + virtual void read_restart(FILE *); + virtual double single(int, int, int, int, double, double, double, double &); void modify_params(int narg, char **arg); double memory_usage(); diff --git a/src/pair_hybrid_scaled.cpp b/src/pair_hybrid_scaled.cpp index 7e1d634627..da8631dcf8 100644 --- a/src/pair_hybrid_scaled.cpp +++ b/src/pair_hybrid_scaled.cpp @@ -14,11 +14,13 @@ #include "pair_hybrid_scaled.h" #include "atom.h" +#include "comm.h" #include "error.h" #include "force.h" #include "input.h" #include "memory.h" #include "respa.h" +#include "suffix.h" #include "update.h" #include "variable.h" @@ -39,7 +41,9 @@ PairHybridScaled::PairHybridScaled(LAMMPS *lmp) PairHybridScaled::~PairHybridScaled() { memory->destroy(fsum); - memory->destroy(scaleval); + memory->destroy(tsum); + delete[] scaleval; + delete[] scaleidx; } /* ---------------------------------------------------------------------- @@ -299,15 +303,20 @@ void PairHybridScaled::settings(int narg, char **arg) ++iarg; if (utils::strmatch(arg[iarg],"^hybrid")) - error->all(FLERR,"Pair style hybrid cannot have hybrid as an argument"); + error->all(FLERR,"Pair style hybrid/scaled cannot have hybrid as an argument"); if (strcmp(arg[iarg],"none") == 0) - error->all(FLERR,"Pair style hybrid cannot have none as an argument"); + error->all(FLERR,"Pair style hybrid/scaled cannot have none as an argument"); styles[nstyles] = force->new_pair(arg[iarg],1,dummy); force->store_style(keywords[nstyles],arg[iarg],0); special_lj[nstyles] = special_coul[nstyles] = nullptr; compute_tally[nstyles] = 1; + if ((((PairHybridScaled *)styles[nstyles])->suffix_flag + & (Suffix::INTEL|Suffix::GPU|Suffix::OMP)) != 0) + error->all(FLERR,"Pair style hybrid/scaled does not support " + "accelerator styles"); + // determine list of arguments for pair style settings // by looking for the next known pair style name. @@ -341,6 +350,60 @@ void PairHybridScaled::settings(int narg, char **arg) flags(); } +/* ---------------------------------------------------------------------- + call sub-style to compute single interaction + error if sub-style does not support single() call + since overlay could have multiple sub-styles, sum results explicitly +------------------------------------------------------------------------- */ + +double PairHybridScaled::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, double &fforce) +{ + if (nmap[itype][jtype] == 0) + error->one(FLERR,"Invoked pair single on pair style none"); + + // update scale values from variables where needed + + const int nvars = scalevars.size(); + if (nvars > 0) { + double *vals = new double[nvars]; + for (i = 0; i < nvars; ++i) { + j = input->variable->find(scalevars[i].c_str()); + vals[i] = input->variable->compute_equal(j); + } + for (i = 0; i < nstyles; ++i) { + if (scaleidx[i] >= 0) + scaleval[i] = vals[scaleidx[i]]; + } + delete[] vals; + } + + double fone; + fforce = 0.0; + double esum = 0.0; + double scale; + + for (int m = 0; m < nmap[itype][jtype]; m++) { + if (rsq < styles[map[itype][jtype][m]]->cutsq[itype][jtype]) { + if (styles[map[itype][jtype][m]]->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 single calls do not support" + " per sub-style special bond values"); + + scale = scaleval[map[itype][jtype][m]]; + esum += scale * styles[map[itype][jtype][m]]->single(i,j,itype,jtype,rsq, + factor_coul,factor_lj,fone); + fforce += scale * fone; + } + } + + if (single_extra) copy_svector(itype,jtype); + return esum; +} + /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ @@ -420,6 +483,63 @@ void PairHybridScaled::coeff(int narg, char **arg) if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairHybridScaled::write_restart(FILE *fp) +{ + PairHybrid::write_restart(fp); + + fwrite(scaleval,sizeof(double),nstyles,fp); + fwrite(scaleidx,sizeof(int),nstyles,fp); + + int n = scalevars.size(); + fwrite(&n,sizeof(int),1,fp); + for (auto var : scalevars) { + n = var.size() + 1; + fwrite(&n,sizeof(int),1,fp); + fwrite(var.c_str(),sizeof(char),n,fp); + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairHybridScaled::read_restart(FILE *fp) +{ + PairHybrid::read_restart(fp); + + delete[] scaleval; + delete[] scaleidx; + scalevars.clear(); + scaleval = new double[nstyles]; + scaleidx = 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); + } + MPI_Bcast(scaleval,nstyles,MPI_DOUBLE,0,world); + MPI_Bcast(scaleidx,nstyles,MPI_INT,0,world); + + char *tmp; + if (me == 0) utils::sfread(FLERR,&n,sizeof(int),1,fp,nullptr,error); + MPI_Bcast(&n,1,MPI_INT,0,world); + scalevars.resize(n); + for (size_t j=0; j < scalevars.size(); ++j) { + if (me == 0) utils::sfread(FLERR,&n,sizeof(int),1,fp,nullptr,error); + MPI_Bcast(&n,1,MPI_INT,0,world); + tmp = new char[n]; + if (me == 0) utils::sfread(FLERR,tmp,sizeof(char),n,fp,nullptr,error); + MPI_Bcast(tmp,n,MPI_CHAR,0,world); + scalevars[j] = tmp; + delete[] tmp; + } +} + /* ---------------------------------------------------------------------- we need to handle Pair::svector special for hybrid/scaled ------------------------------------------------------------------------- */ diff --git a/src/pair_hybrid_scaled.h b/src/pair_hybrid_scaled.h index a8bd5517ec..fbfcdab23a 100644 --- a/src/pair_hybrid_scaled.h +++ b/src/pair_hybrid_scaled.h @@ -35,6 +35,10 @@ class PairHybridScaled : public PairHybrid { virtual void settings(int, char**); virtual void coeff(int, char **); + virtual void write_restart(FILE *); + virtual void read_restart(FILE *); + virtual double single(int, int, int, int, double, double, double, double &); + void init_svector(); void copy_svector(int,int); From 73a33abb44fe8bd2907c044a1652867db3d82e83 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 9 Apr 2021 00:32:06 -0400 Subject: [PATCH 04/14] add unit test for hybrid/scaled --- .../tests/mol-pair-hybrid-scaled.yaml | 104 ++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 unittest/force-styles/tests/mol-pair-hybrid-scaled.yaml diff --git a/unittest/force-styles/tests/mol-pair-hybrid-scaled.yaml b/unittest/force-styles/tests/mol-pair-hybrid-scaled.yaml new file mode 100644 index 0000000000..cc3acb50f0 --- /dev/null +++ b/unittest/force-styles/tests/mol-pair-hybrid-scaled.yaml @@ -0,0 +1,104 @@ +--- +lammps_version: 10 Feb 2021 +date_generated: Fri Feb 26 23:08:45 2021 +epsilon: 5e-14 +skip_tests: gpu intel omp +prerequisites: ! | + atom full + pair lj/cut + pair coul/cut +pre_commands: ! | + variable scale equal 1.0 +post_commands: ! | + pair_modify mix arithmetic +input_file: in.fourmol +pair_style: hybrid/scaled v_scale lj/cut 8.0 0.6 coul/cut 8.0 0.4 coul/cut 8.0 +pair_coeff: ! | + 1 1 lj/cut 0.02 2.5 8 + 1 2 lj/cut 0.01 1.75 8 + 1 3 lj/cut 0.02 2.85 8 + 1 4 lj/cut 0.0173205 2.8 8 + 1 5 lj/cut 0.0173205 2.8 8 + 2 2 lj/cut 0.005 1 8 + 2 3 lj/cut 0.01 2.1 8 + 2 4 lj/cut 0.005 0.5 8 + 2 5 lj/cut 0.00866025 2.05 8 + 3 3 lj/cut 0.02 3.2 8 + 3 4 lj/cut 0.0173205 3.15 8 + 3 5 lj/cut 0.0173205 3.15 8 + 4 4 lj/cut 0.015 3.1 8 + 4 5 lj/cut 0.015 3.1 8 + 5 5 lj/cut 0.015 3.1 8 + * * coul/cut 1 + * * coul/cut 2 +extract: ! "" +natoms: 29 +init_vdwl: 749.237031537357 +init_coul: -127.494586297384 +init_stress: ! |2- + 2.1525607963688685e+03 2.1557327421151899e+03 4.6078904881742919e+03 -7.6038602729615206e+02 1.6844266627640316e+01 6.6957549356541904e+02 +init_forces: ! |2 + 1 -2.1092656751925425e+01 2.6988675971196511e+02 3.3315496490210148e+02 + 2 1.5859534558925552e+02 1.2807631885753918e+02 -1.8817306436807144e+02 + 3 -1.3530454720678361e+02 -3.8712939850050407e+02 -1.4565941679363837e+02 + 4 -7.8195539840070643e+00 2.1451967639963558e+00 -5.9041143405612999e+00 + 5 -2.9163954623584245e+00 -3.3469203159528891e+00 1.2074681739853981e+01 + 6 -8.2989063447195736e+02 9.6019318342576571e+02 1.1479359629470548e+03 + 7 5.7874538635311936e+01 -3.3533985555183068e+02 -1.7140659049826711e+03 + 8 1.4280513303191131e+02 -1.0509295075299345e+02 4.0233495763755388e+02 + 9 8.0984846358566287e+01 7.9600519879262990e+01 3.5197302607961126e+02 + 10 5.3089511229361369e+02 -6.0998478582862322e+02 -1.8376190026890427e+02 + 11 -3.3416993160125812e+00 -4.7792759715873308e+00 -1.0199030124309976e+01 + 12 2.0837574127335213e+01 9.8678992274266921e+00 -6.6547856883058829e+00 + 13 7.7163253261199216e+00 -3.2213746930547997e+00 -1.5767800864580894e-01 + 14 -4.6138299494911639e+00 1.1336312962250332e+00 -8.7660603717255832e+00 + 15 1.6301594996052212e-02 8.3212544078493291e+00 2.0473863128880430e+00 + 16 4.6221152690976908e+02 -3.3124444344467344e+02 -1.1865036959698600e+03 + 17 -4.5568726200724092e+02 3.2159231068141992e+02 1.1980747895060381e+03 + 18 1.2559081069243214e+00 6.6417071126352401e+00 -9.8829024661057083e+00 + 19 1.6184514948299680e+00 -1.6594104323923884e+00 5.6561121961572223e+00 + 20 -3.4526823962510336e+00 -3.1794201827804485e+00 4.2593058942069533e+00 + 21 -6.9075184494915916e+01 -8.0130885501011278e+01 2.1539206802020570e+02 + 22 -1.0659100672969126e+02 -2.5122518903211912e+01 -1.6283765584018167e+02 + 23 1.7515797811309091e+02 1.0400246780074602e+02 -5.2024018223038112e+01 + 24 3.4171625917777114e+01 -2.0194713552213176e+02 1.0982444762500101e+02 + 25 -1.4493448920889654e+02 2.0799041369281703e+01 -1.2091050237305346e+02 + 26 1.0983611557367320e+02 1.8026252731144598e+02 1.2199612526237862e+01 + 27 4.8962849172262665e+01 -2.1594262411895852e+02 8.6423873663236122e+01 + 28 -1.7556665080686602e+02 7.2243004627719102e+01 -1.1798867746650107e+02 + 29 1.2734696054095977e+02 1.4335517724642804e+02 3.2138218235426962e+01 +run_vdwl: 719.583657033589 +run_coul: -127.40544584254 +run_stress: ! |2- + 2.1066855251881925e+03 2.1118463017620702e+03 4.3411898896739367e+03 -7.3939094916433964e+02 3.4004309224046892e+01 6.3091802194682043e+02 +run_forces: ! |2 + 1 -1.8063372896871861e+01 2.6678105157873705e+02 3.2402996659149238e+02 + 2 1.5330358878115447e+02 1.2380492572678898e+02 -1.8151333240574237e+02 + 3 -1.3354888440944052e+02 -3.7931758440809585e+02 -1.4288689214683646e+02 + 4 -7.7881294728555828e+00 2.1395223669670709e+00 -5.8946911982403414e+00 + 5 -2.9015406841040750e+00 -3.3190775902304690e+00 1.2028378254388521e+01 + 6 -8.0488833369818803e+02 9.1802981835006187e+02 1.0244099127408372e+03 + 7 5.5465440662485150e+01 -3.1049131627300432e+02 -1.5711945284966396e+03 + 8 1.3295629283853211e+02 -9.6566834572636509e+01 3.9097872808487460e+02 + 9 7.8594917874857543e+01 7.6787239820699739e+01 3.4114513928465578e+02 + 10 5.2093084326233679e+02 -5.9871672888830824e+02 -1.8144904320802175e+02 + 11 -3.3489474910616370e+00 -4.7299066233626039e+00 -1.0148722292306179e+01 + 12 2.0817110693939330e+01 9.8621648346024777e+00 -6.7801624810903709e+00 + 13 7.6705047254095406e+00 -3.1868508087899996e+00 -1.5820764985177732e-01 + 14 -4.5784791310044675e+00 1.1138053855319887e+00 -8.6502065778611730e+00 + 15 -2.0858645012343142e-03 8.3343285345071436e+00 2.0653788728248101e+00 + 16 4.3381526742384807e+02 -3.1216388880293573e+02 -1.1109931745334770e+03 + 17 -4.2715774864577224e+02 3.0231264864237801e+02 1.1227484174344033e+03 + 18 1.2031503133104606e+00 6.6109154581424221e+00 -9.8172457746610178e+00 + 19 1.6542029696015907e+00 -1.6435312394752812e+00 5.6634735276627497e+00 + 20 -3.4397850729417945e+00 -3.1640002526012512e+00 4.1983600861482540e+00 + 21 -6.8065111490654829e+01 -7.8373161130023504e+01 2.1145341222255522e+02 + 22 -1.0497862711706458e+02 -2.4878742273401844e+01 -1.5988817620288421e+02 + 23 1.7253257365878264e+02 1.0200250230245655e+02 -5.1030905034776815e+01 + 24 3.5759299481226734e+01 -2.0057859782619599e+02 1.1032111627497152e+02 + 25 -1.4570195714964908e+02 2.0679748005808605e+01 -1.2162175868970056e+02 + 26 1.0901403460528100e+02 1.7901644500696690e+02 1.2412674623332103e+01 + 27 4.8035883250870448e+01 -2.1205445789284894e+02 8.4315888267103702e+01 + 28 -1.7229323056476886e+02 7.0823266235363889e+01 -1.1557273097021344e+02 + 29 1.2500312314724302e+02 1.4088629633289813e+02 3.1828931397054006e+01 +... From f2039b56675cce11f90c09b0cd4ad23d0f4b4e00 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 9 Apr 2021 00:36:29 -0400 Subject: [PATCH 05/14] add hybrid/scaled pair style to summary tables --- doc/src/Commands_pair.rst | 2 +- doc/src/pair_style.rst | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index 080f3eff20..e82713f8a4 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -26,6 +26,7 @@ OPT. * :doc:`zero ` * :doc:`hybrid (k) ` * :doc:`hybrid/overlay (k) ` + * :doc:`hybrid/scaled ` * :doc:`kim ` * :doc:`list ` * @@ -33,7 +34,6 @@ OPT. * * * - * * :doc:`adp (o) ` * :doc:`agni (o) ` * :doc:`airebo (io) ` diff --git a/doc/src/pair_style.rst b/doc/src/pair_style.rst index 89530895c4..bc2340d729 100644 --- a/doc/src/pair_style.rst +++ b/doc/src/pair_style.rst @@ -95,6 +95,7 @@ accelerated styles exist. * :doc:`none ` - turn off pairwise interactions * :doc:`hybrid ` - multiple styles of pairwise interactions * :doc:`hybrid/overlay ` - multiple styles of superposed pairwise interactions +* :doc:`hybrid/scaled ` - multiple styles of scaled superposed pairwise interactions * :doc:`zero ` - neighbor list but no interactions * :doc:`adp ` - angular dependent potential (ADP) of Mishin From de158c40ad83a815546a1e07884c1ca37eb2a310 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 9 Apr 2021 00:38:58 -0400 Subject: [PATCH 06/14] add support for torque --- src/pair_hybrid_scaled.cpp | 21 ++++++++++++++++++++- src/pair_hybrid_scaled.h | 2 +- 2 files changed, 21 insertions(+), 2 deletions(-) diff --git a/src/pair_hybrid_scaled.cpp b/src/pair_hybrid_scaled.cpp index da8631dcf8..92cdc767fe 100644 --- a/src/pair_hybrid_scaled.cpp +++ b/src/pair_hybrid_scaled.cpp @@ -31,7 +31,8 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ PairHybridScaled::PairHybridScaled(LAMMPS *lmp) - : PairHybrid(lmp), fsum(nullptr), scaleval(nullptr), scaleidx(nullptr) + : PairHybrid(lmp), fsum(nullptr), tsum(nullptr), + scaleval(nullptr), scaleidx(nullptr) { nmaxfsum = -1; } @@ -92,15 +93,23 @@ void PairHybridScaled::compute(int eflag, int vflag) if (atom->nmax > nmaxfsum) { memory->destroy(fsum); + if (atom->torque_flag) memory->destroy(tsum); nmaxfsum = atom->nmax; memory->create(fsum,nmaxfsum,3,"pair:fsum"); + if (atom->torque_flag) memory->create(tsum,nmaxfsum,3,"pair:tsum"); } const int nall = atom->nlocal + atom->nghost; auto f = atom->f; + auto t = atom->torque; for (i = 0; i < nall; ++i) { fsum[i][0] = f[i][0]; fsum[i][1] = f[i][1]; fsum[i][2] = f[i][2]; + if (atom->torque_flag) { + tsum[i][0] = t[i][0]; + tsum[i][1] = t[i][1]; + tsum[i][2] = t[i][2]; + } } // check if global component of incoming vflag = VIRIAL_FDOTR @@ -147,6 +156,11 @@ void PairHybridScaled::compute(int eflag, int vflag) 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]; + } } restore_special(saved_special); @@ -207,6 +221,11 @@ void PairHybridScaled::compute(int eflag, int vflag) f[i][0] = fsum[i][0]; f[i][1] = fsum[i][1]; f[i][2] = fsum[i][2]; + if (atom->torque_flag) { + t[i][0] = tsum[i][0]; + t[i][1] = tsum[i][1]; + t[i][2] = tsum[i][2]; + } } delete [] saved_special; diff --git a/src/pair_hybrid_scaled.h b/src/pair_hybrid_scaled.h index fbfcdab23a..38a031ad84 100644 --- a/src/pair_hybrid_scaled.h +++ b/src/pair_hybrid_scaled.h @@ -43,7 +43,7 @@ class PairHybridScaled : public PairHybrid { void copy_svector(int,int); protected: - double **fsum; + double **fsum, **tsum; double *scaleval; int *scaleidx; std::vector scalevars; From 4c23ecfd4ff4c856e3b65d48fa4aad3f36409b46 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 9 Apr 2021 00:47:22 -0400 Subject: [PATCH 07/14] error out on special atom styles --- src/pair_hybrid_scaled.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/pair_hybrid_scaled.cpp b/src/pair_hybrid_scaled.cpp index 92cdc767fe..27a02a8d90 100644 --- a/src/pair_hybrid_scaled.cpp +++ b/src/pair_hybrid_scaled.cpp @@ -14,6 +14,7 @@ #include "pair_hybrid_scaled.h" #include "atom.h" +#include "atom_vec.h" #include "comm.h" #include "error.h" #include "force.h" @@ -243,6 +244,9 @@ void PairHybridScaled::settings(int narg, char **arg) error->all(FLERR,fmt::format("Must use pair_style {}/kk with Kokkos", force->pair_style)); + if (atom->avec->forceclearflag) + error->all(FLERR,"Atom style is not compatible with pair_style hybrid/scaled"); + // delete old lists, since cannot just change settings if (nstyles > 0) { From d507c57c8a0b7c22b978e49a03a750433d4073d7 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 9 Apr 2021 09:45:09 -0400 Subject: [PATCH 08/14] add missing zeroing of the torque array, reformat code --- src/pair_hybrid_scaled.cpp | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/pair_hybrid_scaled.cpp b/src/pair_hybrid_scaled.cpp index 27a02a8d90..62b94d9ba9 100644 --- a/src/pair_hybrid_scaled.cpp +++ b/src/pair_hybrid_scaled.cpp @@ -134,9 +134,10 @@ void PairHybridScaled::compute(int eflag, int vflag) for (m = 0; m < nstyles; m++) { - // clear forces + // clear forces and torques memset(&f[0][0],0,nall*3*sizeof(double)); + if (atom->torque_flag) memset(&t[0][0],0,nall*3*sizeof(double)); set_special(m); @@ -171,17 +172,17 @@ void PairHybridScaled::compute(int eflag, int vflag) if (respaflag && !respa->tally_global) continue; if (eflag_global) { - eng_vdwl += scale * styles[m]->eng_vdwl; - eng_coul += scale * styles[m]->eng_coul; + eng_vdwl += scale*styles[m]->eng_vdwl; + eng_coul += scale*styles[m]->eng_coul; } if (vflag_global) { - for (n = 0; n < 6; n++) virial[n] += scale * styles[m]->virial[n]; + for (n = 0; n < 6; n++) virial[n] += scale*styles[m]->virial[n]; } if (eflag_atom) { n = atom->nlocal; if (force->newton_pair) n += atom->nghost; double *eatom_substyle = styles[m]->eatom; - for (i = 0; i < n; i++) eatom[i] += scale * eatom_substyle[i]; + for (i = 0; i < n; i++) eatom[i] += scale*eatom_substyle[i]; } if (vflag_atom) { n = atom->nlocal; @@ -189,7 +190,7 @@ void PairHybridScaled::compute(int eflag, int vflag) double **vatom_substyle = styles[m]->vatom; for (i = 0; i < n; i++) for (j = 0; j < 6; j++) - vatom[i][j] += scale * vatom_substyle[i][j]; + vatom[i][j] += scale*vatom_substyle[i][j]; } // substyles may be CENTROID_SAME or CENTROID_AVAIL @@ -201,22 +202,22 @@ void PairHybridScaled::compute(int eflag, int vflag) double **cvatom_substyle = styles[m]->cvatom; for (i = 0; i < n; i++) for (j = 0; j < 9; j++) - cvatom[i][j] += scale * cvatom_substyle[i][j]; + cvatom[i][j] += scale*cvatom_substyle[i][j]; } else { double **vatom_substyle = styles[m]->vatom; for (i = 0; i < n; i++) { for (j = 0; j < 6; j++) { - cvatom[i][j] += scale * vatom_substyle[i][j]; + cvatom[i][j] += scale*vatom_substyle[i][j]; } for (j = 6; j < 9; j++) { - cvatom[i][j] += scale * vatom_substyle[i][j-3]; + cvatom[i][j] += scale*vatom_substyle[i][j-3]; } } } } } - // copy accumulated forces to original force array + // copy accumulated scaled forces to original force array for (i = 0; i < nall; ++i) { f[i][0] = fsum[i][0]; From 0496fd27dbb8cace53090f3e4fb81abdad489816 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 9 Apr 2021 22:49:25 -0400 Subject: [PATCH 09/14] reorder include files --- src/pair_hybrid.cpp | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index c12ec67351..48ba1ccff7 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -14,20 +14,19 @@ #include "pair_hybrid.h" -#include -#include #include "atom.h" -#include "force.h" -#include "pair.h" -#include "neighbor.h" -#include "neigh_request.h" -#include "update.h" #include "comm.h" -#include "memory.h" #include "error.h" +#include "force.h" +#include "memory.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "pair.h" #include "respa.h" - #include "suffix.h" +#include "update.h" + +#include using namespace LAMMPS_NS; From 0d325f22219239a52d3b2b7ca8344749814b51bc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 9 Apr 2021 22:50:16 -0400 Subject: [PATCH 10/14] error out when scale factor variables do not exist --- src/pair_hybrid_scaled.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/pair_hybrid_scaled.cpp b/src/pair_hybrid_scaled.cpp index 62b94d9ba9..7a30c9a3ee 100644 --- a/src/pair_hybrid_scaled.cpp +++ b/src/pair_hybrid_scaled.cpp @@ -71,6 +71,9 @@ void PairHybridScaled::compute(int eflag, int vflag) double *vals = new double[nvars]; for (i = 0; i < nvars; ++i) { j = input->variable->find(scalevars[i].c_str()); + if (j < 0) + error->all(FLERR,fmt::format("Variable '{}' not found when updating " + "scale factors",scalevars[i])); vals[i] = input->variable->compute_equal(j); } for (i = 0; i < nstyles; ++i) { From ec6e2d35cbd46fe5cdbcb3bacfebc0bd6c06dbf6 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 9 Apr 2021 22:50:51 -0400 Subject: [PATCH 11/14] complete update of the hybrid documentation for hybrid/scaled --- doc/src/pair_hybrid.rst | 157 +++++++++++++++++++++------------------- 1 file changed, 83 insertions(+), 74 deletions(-) diff --git a/doc/src/pair_hybrid.rst b/doc/src/pair_hybrid.rst index c851423029..9fdba318d9 100644 --- a/doc/src/pair_hybrid.rst +++ b/doc/src/pair_hybrid.rst @@ -59,35 +59,40 @@ be assigned to each pair of atom types. The assignment of pair styles to type pairs is made via the :doc:`pair_coeff ` command. The *hybrid/scaled* style differs from the *hybrid/overlay* style by requiring a factor for each pair style that is used to scale all -forces and energies computed by the pair style. +forces, energies and stresses computed by each sub-style. Because of +the additional complexity, the *hybrid/scaled* style will have more +overhead and thus will be a bit slower than *hybrid/overlay*. Here are two examples of hybrid simulations. The *hybrid* style could -be used for a simulation of a metal droplet on a LJ surface. The -metal atoms interact with each other via an *eam* potential, the -surface atoms interact with each other via a *lj/cut* potential, and -the metal/surface interaction is also computed via a *lj/cut* -potential. The *hybrid/overlay* style could be used as in the second -example above, where multiple potentials are superposed in an additive -fashion to compute the interaction between atoms. In this example, -using *lj/cut* and *coul/long* together gives the same result as if -the *lj/cut/coul/long* potential were used by itself. In this case, -it would be more efficient to use the single combined potential, but -in general any combination of pair potentials can be used together in -to produce an interaction that is not encoded in any single pair_style -file, e.g. adding Coulombic forces between granular particles. -The *hybrid/scaled* style enables more complex combinations of pair -styles than a simple sum as *hybrid/overlay* does. Furthermore, since -the scale factors can be variables, they can change during a simulation -which would allow to smoothly switch between two different pair styles -or two different parameter sets. +be used for a simulation of a metal droplet on a LJ surface. The metal +atoms interact with each other via an *eam* potential, the surface atoms +interact with each other via a *lj/cut* potential, and the metal/surface +interaction is also computed via a *lj/cut* potential. The +*hybrid/overlay* style could be used as in the second example above, +where multiple potentials are superposed in an additive fashion to +compute the interaction between atoms. In this example, using *lj/cut* +and *coul/long* together gives the same result as if the +*lj/cut/coul/long* potential were used by itself. In this case, it +would be more efficient to use the single combined potential, but in +general any combination of pair potentials can be used together in to +produce an interaction that is not encoded in any single pair_style +file, e.g. adding Coulombic forces between granular particles. The +*hybrid/scaled* style enables more complex combinations of pair styles +than a simple sum as *hybrid/overlay* does; there may be fractional +contributions from sub-styles or contributions may be subtracted with a +negative scale factor. Furthermore, since the scale factors can be +variables that may change during a simulation, which would allow, for +instance, to smoothly switch between two different pair styles or two +different parameter sets. 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 its scale -factor. The scale factor may be an equal style (or equivalent) -variable. Each sub-style's name is followed by its usual arguments, as -illustrated in the example above. See the doc pages of individual pair -styles for a listing and explanation of the appropriate arguments. +*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 +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 individual pair styles for a listing and explanation of the +appropriate arguments. Note that an individual pair style can be used multiple times as a sub-style. For efficiency this should only be done if your model @@ -164,7 +169,7 @@ one sub-style. Just as with a simulation using a single pair style, if you specify the same atom type pair in a second pair_coeff command, the previous assignment will be overwritten. -For the *hybrid/overlay* and *hybrid/scaled* style, each atom type pair +For the *hybrid/overlay* and *hybrid/scaled* styles, each atom type pair I,J can be assigned to one or more sub-styles. If you specify the same atom type pair in a second pair_coeff command with a new sub-style, then the second sub-style is added to the list of potentials that will be @@ -187,15 +192,15 @@ same: Coefficients must be defined for each pair of atoms types via the :doc:`pair_coeff ` command as described above, or in the -data file or restart files read by the :doc:`read_data ` or -:doc:`read_restart ` commands, or by mixing as described -below. +data file read by the :doc:`read_data ` commands, or by +mixing as described below. For all of the *hybrid*, *hybrid/overlay*, and *hybrid/scaled* styles, every atom type pair I,J (where I <= J) must be assigned to at least one sub-style via the :doc:`pair_coeff ` command as in the examples above, or in the data file read by the :doc:`read_data -`, or by mixing as described below. +`, or by mixing as described below. Also all sub-styles +must be used at least once in a :doc:`pair_coeff ` command. If you want there to be no interactions between a particular pair of atom types, you have 3 choices. You can assign the type pair to some @@ -208,22 +213,22 @@ input script: .. code-block:: LAMMPS - pair_coeff 2 3 none + pair_coeff 2 3 none or this form in the "Pair Coeffs" section of the data file: .. parsed-literal:: - 3 none + 3 none If an assignment to *none* is made in a simulation with the -*hybrid/overlay* pair style, it wipes out all previous assignments of -that atom type pair to sub-styles. +*hybrid/overlay* or *hybrid/scaled* pair style, it wipes out all +previous assignments of that pair of atom types to sub-styles. Note that you may need to use an :doc:`atom_style ` hybrid command in your input script, if atoms in the simulation will need attributes from several atom styles, due to using multiple pair -potentials. +styles with different requirements. ---------- @@ -232,8 +237,9 @@ for applying weightings that change the strength of pairwise interactions between pairs of atoms that are also 1-2, 1-3, and 1-4 neighbors in the molecular bond topology, as normally set by the :doc:`special_bonds ` command. Different weights can be -assigned to different pair hybrid sub-styles via the :doc:`pair_modify special ` command. This allows multiple force fields -to be used in a model of a hybrid system, however, there is no consistent +assigned to different pair hybrid sub-styles via the :doc:`pair_modify +special ` command. This allows multiple force fields to be +used in a model of a hybrid system, however, there is no consistent approach to determine parameters automatically for the interactions between the two force fields, this is only recommended when particles described by the different force fields do not mix. @@ -307,28 +313,27 @@ Pair_style hybrid allows interactions between type pairs 2-2, 1-2, could even add a second interaction for 1-1 to be computed by another 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 interactions, e.g. exclude 1-1 or 1-3 interactions. That -is not conceptually well-defined for many-body interactions, since the +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 +interactions, e.g. exclude 1-1 or 1-3 interactions. That is not +conceptually well-defined for many-body interactions, since the potential will typically calculate energies and foces for small groups of atoms, e.g. 3 or 4 atoms, using the neighbor lists of the atoms to -find the additional atoms in the group. It is typically non-physical -to think of excluding an interaction between a particular pair of -atoms when the potential computes 3-body or 4-body interactions. +find the additional atoms in the group. However, you can still use the pair_coeff none setting or the :doc:`neigh_modify exclude ` command to exclude certain type pairs from the neighbor list that will be passed to a many-body sub-style. This will alter the calculations made by a many-body -potential, since it builds its list of 3-body, 4-body, etc -interactions from the pair list. You will need to think carefully as -to whether it produces a physically meaningful result for your model. +potential beyond the specific pairs, since it builds its list of 3-body, +4-body, etc interactions from the pair lists. You will need to think +**carefully** as to whether excluding such pairs produces a physically +meaningful result for your model. For example, imagine you have two atom types in your model, type 1 for atoms in one surface, and type 2 for atoms in the other, and you wish to use a Tersoff potential to compute interactions within each -surface, but not between surfaces. Then either of these two command +surface, but not between the surfaces. Then either of these two command sequences would implement that model: .. code-block:: LAMMPS @@ -345,9 +350,9 @@ Either way, only neighbor lists with 1-1 or 2-2 interactions would be 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, using hybrid/overlay, to use 2 many-body -potentials together, in an overlapping manner. Imagine you have CNT -(C atoms) on a Si surface. You want to use Tersoff for Si/Si and Si/C +Here is another example to use 2 many-body potentials together in an +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: @@ -358,9 +363,9 @@ atoms are type 2. Something like this will work: pair_coeff * * airebo CH.airebo NULL C Note that to prevent the Tersoff potential from computing C/C -interactions, you would need to modify the SiC.tersoff file to turn -off C/C interaction, i.e. by setting the appropriate coefficients to -0.0. +interactions, you would need to **modify** the SiC.tersoff potential +file to turn off C/C interaction, i.e. by setting the appropriate +coefficients to 0.0. ---------- @@ -368,18 +373,19 @@ Styles with a *gpu*\ , *intel*\ , *kk*\ , *omp*\ , or *opt* suffix are functionally the same as the corresponding style without the suffix. They have been optimized to run faster, depending on your available hardware, as discussed on the :doc:`Speed packages ` doc -page. +page. Pair style *hybrid/scaled* does (currently) not support the +*gpu*, *omp*, *kk*, or *intel* suffix. -Since the *hybrid* and *hybrid/overlay* styles delegate computation to -the individual sub-styles, the suffix versions of the *hybrid* and -*hybrid/overlay* styles are used to propagate the corresponding suffix -to all sub-styles, if those versions exist. Otherwise the -non-accelerated version will be used. +Since the *hybrid*, *hybrid/overlay*, *hybrid/scaled* styles delegate +computation to the individual sub-styles, the suffix versions of the +*hybrid* and *hybrid/overlay* styles are used to propagate the +corresponding suffix to all sub-styles, if those versions +exist. Otherwise the non-accelerated version will be used. -The individual accelerated sub-styles are part of the GPU, USER-OMP -and OPT packages, respectively. They are only enabled if LAMMPS was -built with those packages. See the :doc:`Build package ` -doc page for more info. +The individual accelerated sub-styles are part of the GPU, KOKKOS, +USER-INTEL, USER-OMP, and OPT packages, respectively. They are only +enabled if LAMMPS was built with those packages. See the :doc:`Build +package ` doc page for more info. You can specify the accelerated styles explicitly in your input script by including their suffix, or you can use the :doc:`-suffix command-line switch ` when you invoke LAMMPS, or you can use the @@ -397,17 +403,17 @@ Any pair potential settings made via the :doc:`pair_modify ` command are passed along to all sub-styles of the hybrid potential. -For atom type pairs I,J and I != J, if the sub-style assigned to I,I -and J,J is the same, and if the sub-style allows for mixing, then the +For atom type pairs I,J and I != J, if the sub-style assigned to I,I and +J,J is the same, and if the sub-style allows for mixing, then the coefficients for I,J can be mixed. This means you do not have to specify a pair_coeff command for I,J since the I,J type pair will be -assigned automatically to the sub-style defined for both I,I and J,J -and its coefficients generated by the mixing rule used by that -sub-style. For the *hybrid/overlay* style, there is an additional -requirement that both the I,I and J,J pairs are assigned to a single -sub-style. See the "pair_modify" command for details of mixing rules. -See the See the doc page for the sub-style to see if allows for -mixing. +assigned automatically to the sub-style defined for both I,I and J,J and +its coefficients generated by the mixing rule used by that sub-style. +For the *hybrid/overlay* and *hybrid/scaled* style, there is an +additional requirement that both the I,I and J,J pairs are assigned to a +single sub-style. See the :doc:`pair_modify ` command for +details of mixing rules. See the See the doc page for the sub-style to +see if allows for mixing. The hybrid pair styles supports the :doc:`pair_modify ` shift, table, and tail options for an I,J pair interaction, if the @@ -418,7 +424,10 @@ settings are written to :doc:`binary restart files `, so a :doc:`pair_style ` command does not need to specified in an input script that reads a restart file. However, the coefficient information is not stored in the restart file. Thus, pair_coeff -commands need to be re-specified in the restart input script. +commands need to be re-specified in the restart input script. For pair +style *hybrid/scaled* also the names of any variables used as scale +factors are restored, but not the variables themselves, so those may +need to be redefined when continuing from a restart. These pair styles support the use of the *inner*\ , *middle*\ , and *outer* keywords of the :doc:`run_style respa ` command, if @@ -435,7 +444,7 @@ short-range Coulombic cutoff used by each of these long pair styles is the same or else LAMMPS will generate an error. Pair style *hybrid/scaled* currently only works for non-accelerated -pair styles. +pair styles and pair styles from the OPT package. Related commands """""""""""""""" From d88cf587b2743e1b94436acc293a14ad5def69ee Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 10 Apr 2021 00:25:00 -0400 Subject: [PATCH 12/14] add pair style coul/cut/global and fix restart/data bugs in coul/cut --- doc/src/Commands_pair.rst | 1 + doc/src/pair_coul.rst | 14 ++- doc/src/pair_style.rst | 1 + src/pair_coul_cut.cpp | 61 ++++++++-- src/pair_coul_cut.h | 6 +- src/pair_coul_cut_global.cpp | 44 +++++++ src/pair_coul_cut_global.h | 55 +++++++++ unittest/force-styles/test_pair_style.cpp | 1 - .../force-styles/tests/mol-pair-coul_cut.yaml | 114 +++++++++--------- .../tests/mol-pair-coul_cut_global.yaml | 85 +++++++++++++ 10 files changed, 311 insertions(+), 71 deletions(-) create mode 100644 src/pair_coul_cut_global.cpp create mode 100644 src/pair_coul_cut_global.h create mode 100644 unittest/force-styles/tests/mol-pair-coul_cut_global.yaml diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index e82713f8a4..60996b9c65 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -69,6 +69,7 @@ OPT. * :doc:`comb3 ` * :doc:`cosine/squared ` * :doc:`coul/cut (gko) ` + * :doc:`coul/cut/global (o) ` * :doc:`coul/cut/soft (o) ` * :doc:`coul/debye (gko) ` * :doc:`coul/diel (o) ` diff --git a/doc/src/pair_coul.rst b/doc/src/pair_coul.rst index 4def43647f..b8303622aa 100644 --- a/doc/src/pair_coul.rst +++ b/doc/src/pair_coul.rst @@ -10,6 +10,8 @@ .. index:: pair_style coul/dsf/gpu .. index:: pair_style coul/dsf/kk .. index:: pair_style coul/dsf/omp +.. index:: pair_style coul/cut/global +.. index:: pair_style coul/cut/global/omp .. index:: pair_style coul/long .. index:: pair_style coul/long/omp .. index:: pair_style coul/long/kk @@ -40,6 +42,11 @@ pair_style coul/dsf command Accelerator Variants: *coul/dsf/gpu*, *coul/dsf/kk*, *coul/dsf/omp* +pair_style coul/cut/global command +================================== + +Accelerator Variants: *coul/cut/omp* + pair_style coul/long command ============================ @@ -76,8 +83,8 @@ Syntax pair_style coul/cut cutoff pair_style coul/debye kappa cutoff pair_style coul/dsf alpha cutoff + pair_style coul/cut/global cutoff pair_style coul/long cutoff - pair_style coul/long/gpu cutoff pair_style coul/wolf alpha cutoff pair_style coul/streitz cutoff keyword alpha pair_style tip4p/cut otype htype btype atype qdist cutoff @@ -245,6 +252,11 @@ Streitz-Mintmire parameterization for the material being modeled. ---------- +Pair style *coul/cut/global* computes the same Coulombic interactions +as style *coul/cut* except that it allows only a single global cutoff +and thus makes it compatible for use in combination with long-range +coulomb styles in :doc:`hybrid pair styles `. + Styles *coul/long* and *coul/msm* compute the same Coulombic interactions as style *coul/cut* except that an additional damping factor is applied so it can be used in conjunction with the diff --git a/doc/src/pair_style.rst b/doc/src/pair_style.rst index bc2340d729..49eac18aa8 100644 --- a/doc/src/pair_style.rst +++ b/doc/src/pair_style.rst @@ -133,6 +133,7 @@ accelerated styles exist. * :doc:`comb3 ` - charge-optimized many-body (COMB3) potential * :doc:`cosine/squared ` - Cooke-Kremer-Deserno membrane model potential * :doc:`coul/cut ` - cutoff Coulomb potential +* :doc:`coul/cut/global ` - cutoff Coulomb potential * :doc:`coul/cut/soft ` - Coulomb potential with a soft core * :doc:`coul/debye ` - cutoff Coulomb potential with Debye screening * :doc:`coul/diel ` - Coulomb potential with dielectric permittivity diff --git a/src/pair_coul_cut.cpp b/src/pair_coul_cut.cpp index 44ddd424d9..428c12c2e0 100644 --- a/src/pair_coul_cut.cpp +++ b/src/pair_coul_cut.cpp @@ -13,22 +13,24 @@ #include "pair_coul_cut.h" -#include -#include #include "atom.h" #include "comm.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "memory.h" #include "error.h" +#include "force.h" +#include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" +#include +#include using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairCoulCut::PairCoulCut(LAMMPS *lmp) : Pair(lmp) {} +PairCoulCut::PairCoulCut(LAMMPS *lmp) : Pair(lmp) { + writedata = 1; +} /* ---------------------------------------------------------------------- */ @@ -208,8 +210,10 @@ void PairCoulCut::init_style() double PairCoulCut::init_one(int i, int j) { - if (setflag[i][j] == 0) + if (setflag[i][j] == 0) { cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + scale[i][j] = 1.0; + } scale[j][i] = scale[i][j]; @@ -225,11 +229,15 @@ void PairCoulCut::write_restart(FILE *fp) write_restart_settings(fp); int i,j; - for (i = 1; i <= atom->ntypes; i++) + for (i = 1; i <= atom->ntypes; i++) { for (j = i; j <= atom->ntypes; j++) { + fwrite(&scale[i][j],sizeof(double),1,fp); fwrite(&setflag[i][j],sizeof(int),1,fp); - if (setflag[i][j]) fwrite(&cut[i][j],sizeof(double),1,fp); + if (setflag[i][j]) { + fwrite(&cut[i][j],sizeof(double),1,fp); + } } + } } /* ---------------------------------------------------------------------- @@ -243,15 +251,21 @@ void PairCoulCut::read_restart(FILE *fp) int i,j; int me = comm->me; - for (i = 1; i <= atom->ntypes; i++) + for (i = 1; i <= atom->ntypes; i++) { for (j = i; j <= atom->ntypes; j++) { - if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); + if (me == 0) { + utils::sfread(FLERR,&scale[i][j],sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); + } + MPI_Bcast(&scale[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); if (setflag[i][j]) { - if (me == 0) utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error); + if (me == 0) + utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error); MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); } } + } } /* ---------------------------------------------------------------------- @@ -281,6 +295,27 @@ void PairCoulCut::read_restart_settings(FILE *fp) MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairCoulCut::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp,"%d\n",i); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairCoulCut::write_data_all(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) + fprintf(fp,"%d %d %g\n",i,j,cut[i][j]); +} + /* ---------------------------------------------------------------------- */ double PairCoulCut::single(int i, int j, int /*itype*/, int /*jtype*/, diff --git a/src/pair_coul_cut.h b/src/pair_coul_cut.h index 5fed8344e2..fe44c0d8a4 100644 --- a/src/pair_coul_cut.h +++ b/src/pair_coul_cut.h @@ -30,15 +30,17 @@ class PairCoulCut : public Pair { virtual ~PairCoulCut(); virtual void compute(int, int); virtual void settings(int, char **); - void coeff(int, char **); + virtual void coeff(int, char **); void init_style(); double init_one(int, int); void write_restart(FILE *); void read_restart(FILE *); virtual void write_restart_settings(FILE *); virtual void read_restart_settings(FILE *); + virtual void write_data(FILE *); + virtual void write_data_all(FILE *); virtual double single(int, int, int, int, double, double, double, double &); - void *extract(const char *, int &); + virtual void *extract(const char *, int &); protected: double cut_global; diff --git a/src/pair_coul_cut_global.cpp b/src/pair_coul_cut_global.cpp new file mode 100644 index 0000000000..fa3d394092 --- /dev/null +++ b/src/pair_coul_cut_global.cpp @@ -0,0 +1,44 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://lammps.sandia.gov/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "pair_coul_cut_global.h" + +#include "error.h" + +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairCoulCutGlobal::coeff(int narg, char **arg) +{ + if (narg != 2) + error->all(FLERR,"Incorrect args for pair coefficients"); + + PairCoulCut::coeff(narg,arg); +} + + +/* ---------------------------------------------------------------------- */ + +void *PairCoulCutGlobal::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"cut_coul") == 0) return (void *) &cut_global; + dim = 2; + if (strcmp(str,"scale") == 0) return (void *) scale; + return nullptr; +} diff --git a/src/pair_coul_cut_global.h b/src/pair_coul_cut_global.h new file mode 100644 index 0000000000..4b3a8ddbaa --- /dev/null +++ b/src/pair_coul_cut_global.h @@ -0,0 +1,55 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(coul/cut/global,PairCoulCutGlobal) + +#else + +#ifndef LMP_PAIR_COUL_CUT_GLOBAL_H +#define LMP_PAIR_COUL_CUT_GLOBAL_H + +#include "pair_coul_cut.h" + +namespace LAMMPS_NS { + +class PairCoulCutGlobal : public PairCoulCut { + public: + PairCoulCutGlobal(class LAMMPS *lmp) : PairCoulCut(lmp) {} + void coeff(int, char **); + void *extract(const char *, int &); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair style coul/cut requires atom attribute q + +The atom style defined does not have these attributes. + +*/ diff --git a/unittest/force-styles/test_pair_style.cpp b/unittest/force-styles/test_pair_style.cpp index 5ea7fc371c..057f2b5352 100644 --- a/unittest/force-styles/test_pair_style.cpp +++ b/unittest/force-styles/test_pair_style.cpp @@ -25,7 +25,6 @@ #include "atom.h" #include "compute.h" -#include "fmt/format.h" #include "force.h" #include "info.h" #include "input.h" diff --git a/unittest/force-styles/tests/mol-pair-coul_cut.yaml b/unittest/force-styles/tests/mol-pair-coul_cut.yaml index 921e682a9a..09e97a3a8d 100644 --- a/unittest/force-styles/tests/mol-pair-coul_cut.yaml +++ b/unittest/force-styles/tests/mol-pair-coul_cut.yaml @@ -1,7 +1,8 @@ --- -lammps_version: 10 Feb 2021 -date_generated: Fri Feb 26 23:08:42 2021 +lammps_version: 8 Apr 2021 +date_generated: Sat Apr 10 00:22:24 2021 epsilon: 1e-13 +skip_tests: prerequisites: ! | atom full pair coul/cut @@ -11,17 +12,22 @@ post_commands: ! | input_file: in.fourmol pair_style: coul/cut 8.0 pair_coeff: ! | - * * + 1 1 + 2 2 7.0 + 2 3 7.5 + 3 3 + 4 4 8.0 + 5 5 extract: ! | cut_coul 2 natoms: 29 init_vdwl: 0 -init_coul: -127.494586297384 +init_coul: -93.55664334073657 init_stress: ! |- - -2.6824547034957277e+01 -4.3162775104089022e+01 -5.7507264158338124e+01 -8.2055093564543602e-01 -7.9072701063929465e+00 3.0552092014126533e+00 + -2.0535385669481379e+01 -1.6145387915767863e+01 -5.6875869755487301e+01 -4.1550674863524968e+00 -7.0519267473420957e+00 2.6305371873347489e+00 init_forces: ! |2 1 2.2407335289699457e+00 -5.8916421262140251e-02 4.2668639853362533e-01 - 2 3.0979928510686561e-01 -2.1737695778193569e+00 -1.8762407787145119e+00 + 2 -1.2152150177208687e-01 -3.4650575893033708e+00 -1.3360477574570087e+00 3 -1.5509825092714681e-02 -8.6264917303471161e-02 3.0367577692877928e-02 4 5.1555686582278520e-02 1.0144901459000669e-02 -3.0866112200648410e-01 5 -3.9871973553557233e-01 7.0523075214913145e-01 -7.8022318023030612e-02 @@ -30,57 +36,57 @@ init_forces: ! |2 8 -1.7087898110045159e+00 4.1818078578969917e+00 2.4290227133489601e+00 9 1.8279010754687195e+00 -5.6724899047236770e+00 1.6512690951588267e+00 10 -2.9363989744075952e-01 4.2512276557131090e-01 -2.0317384271111175e-01 - 11 -9.8868358921602750e-01 1.1284880357946445e+00 -5.3995772876854398e-01 + 11 -9.6854530039287257e-01 1.2818484089934798e+00 -4.0964874266313867e-01 12 3.3104189805348057e+00 -7.6522029601079322e-01 1.2706541181424338e+00 - 13 -3.8231563183336364e-01 -1.1565866576649098e-02 -8.7140102078706039e-03 - 14 -1.2285578202646528e+00 4.4726948380599840e-01 -1.5341285475604788e-02 - 15 2.2085158688210516e-01 -1.6336214445566244e-01 -9.6577522905262758e-01 + 13 -3.0889352958622951e-01 2.5103266106999866e-01 1.0669227367631845e-01 + 14 -1.2634615134073266e+00 5.9356329877397740e-01 2.0097321421686021e-01 + 15 3.4164194917662777e-01 -1.3968217225075463e-02 -1.0754490238146823e+00 16 -1.0515762083518698e+00 -3.6728607969167482e-01 2.7987602083946292e+00 - 17 -2.3442527695736954e+00 6.0494781225943264e+00 -7.7669898420813883e+00 - 18 1.2747707302049942e+00 6.6751091235431907e+00 -9.9139029454048035e+00 - 19 1.6181330640335625e+00 -1.6591712461142003e+00 5.6543694708933723e+00 - 20 -3.4516847879389365e+00 -3.1783992642977736e+00 4.2589367844805865e+00 - 21 2.4909407783496178e+00 1.4847928289095484e+00 -1.0503546063193097e+01 - 22 1.4973505700836185e+00 1.0712683327319574e+00 6.7413935914321712e+00 - 23 -4.4865766347042468e+00 -3.8185091520236116e+00 4.2817682561025867e+00 - 24 -2.4197806588078703e+00 9.8687406957238029e+00 -2.3585711007230041e+00 - 25 3.5804022684917332e+00 -3.1080767536674005e+00 3.9458463586096970e+00 - 26 -2.0751789623087848e+00 -7.6352193307876766e+00 -4.5052467808202223e-01 - 27 -2.8475395052833852e+00 1.1111959093179417e+01 -4.4252374188335644e+00 - 28 4.8464204075747359e+00 -5.2910383050538492e+00 4.0808901405648896e+00 - 29 -1.2636120082903082e+00 -6.1719354975176035e+00 9.2219267915908476e-01 + 17 -2.4272925242862877e+00 5.1482976415440387e+00 -6.9745861917436480e+00 + 18 -3.3422585449197340e-01 1.0468316830549977e+01 -9.5556804718861983e+00 + 19 1.8166530000757015e+00 -1.4461616747597046e+00 7.5536140661465714e+00 + 20 -3.2356613999008608e+00 -3.8668819868868143e+00 5.3737386596536343e+00 + 21 4.0222635030104259e+00 1.6132349119387515e+00 -1.0241718963154419e+01 + 22 1.4374544034911272e+00 -9.5656370964487103e-01 6.5980789669241569e+00 + 23 -4.4376821117692176e+00 -3.8947548010087796e+00 4.1208787227921038e+00 + 24 -1.5247815933638265e+00 1.1713843272625617e+01 -6.0248223323521604e+00 + 25 3.1898863386468577e+00 -4.1364515802851436e+00 3.9097585091440594e+00 + 26 -2.2629415636619381e+00 -8.5609820249335087e+00 -1.1581224704665538e+00 + 27 -3.6829593583586151e+00 1.0444039166475267e+01 -3.7821850434646627e+00 + 28 5.4031537104784855e+00 -4.4754845440910707e+00 3.6816742644554843e+00 + 29 -1.2926003313820287e+00 -6.0718114858636643e+00 7.3448520698640607e-02 run_vdwl: 0 -run_coul: -127.551092913469 +run_coul: -93.61491356191763 run_stress: ! |- - -2.6928416018548067e+01 -4.3177725729834854e+01 -5.7444951165085541e+01 -8.5814323779813129e-01 -7.9912081348585087e+00 2.9897259857467389e+00 + -2.0637957613066451e+01 -1.6150296923318230e+01 -5.6826659025532969e+01 -4.1943798944102859e+00 -7.1336719665040942e+00 2.5570315775315215e+00 run_forces: ! |2 - 1 2.2424368291189611e+00 -5.7157845887417236e-02 4.3210281185121069e-01 - 2 3.0183044940246062e-01 -2.1770211242302246e+00 -1.8796316672405435e+00 - 3 -1.5579134024760494e-02 -8.6577075665555045e-02 3.0374580927846943e-02 - 4 5.2605597845996610e-02 1.0081774332824270e-02 -3.0969764881211892e-01 - 5 -3.9802063955902828e-01 7.0541551998225649e-01 -7.9105180426906119e-02 - 6 2.0115016302231399e+00 -3.7433268305651115e+00 -2.9688564502920158e+00 - 7 -3.2650169036425236e-01 7.4540116984953042e-01 3.8864116015671808e+00 - 8 -1.6992053224967520e+00 4.1770597030834642e+00 2.4329396564406349e+00 - 9 1.8250267620013743e+00 -5.6775189439742455e+00 1.6562371734234631e+00 - 10 -2.9411864176775349e-01 4.2484150168808632e-01 -2.0454722162066086e-01 - 11 -9.8864600485826815e-01 1.1292383541031983e+00 -5.3894807546123469e-01 - 12 3.3107287761242974e+00 -7.6414239662660099e-01 1.2747075215620911e+00 - 13 -3.8336813011366144e-01 -1.2486400314578122e-02 -9.4718485253563467e-03 - 14 -1.2290837060322926e+00 4.4862678961891689e-01 -1.5069551071827086e-02 - 15 2.2244012017801007e-01 -1.6546887582342559e-01 -9.6944707732754809e-01 - 16 -1.0531155341397860e+00 -3.6654926990083037e-01 2.7960996028622298e+00 - 17 -2.3452245945824948e+00 6.0586966641613316e+00 -7.7701709031633026e+00 - 18 1.2287925114358347e+00 6.6283596097346500e+00 -9.8694991818018014e+00 - 19 1.6548495610099294e+00 -1.6299328734082061e+00 5.6681416753764280e+00 - 20 -3.4405939298489487e+00 -3.1595057592530829e+00 4.2031829943411507e+00 - 21 2.4999236613139768e+00 1.4581441593082431e+00 -1.0498791509490704e+01 - 22 1.5230738475837411e+00 1.0948611353935194e+00 6.7461577757855986e+00 - 23 -4.5208030817248321e+00 -3.8176488387685179e+00 4.2715025892256220e+00 - 24 -2.4424910840260479e+00 9.8889784537097576e+00 -2.3784455147561552e+00 - 25 3.6199382208005479e+00 -3.1023007862101899e+00 3.9803408068580102e+00 - 26 -2.0901080780170158e+00 -7.6586474495008154e+00 -4.6300658746615819e-01 - 27 -2.8661140489132810e+00 1.1127847846706011e+01 -4.4084385064798797e+00 - 28 4.8657388732889295e+00 -5.2969930881855793e+00 4.0790567237989928e+00 - 29 -1.2659132198580254e+00 -6.1822751233574085e+00 9.0587140991575621e-01 + 1 2.2425199864949272e+00 -5.7185268117056043e-02 4.3203484436437956e-01 + 2 -1.2953282231918894e-01 -3.4675950045872330e+00 -1.3398768542921780e+00 + 3 -1.5576933865004368e-02 -8.6576616933096762e-02 3.0374805119446222e-02 + 4 5.2594937495592908e-02 1.0083837924927244e-02 -3.0970246073870022e-01 + 5 -3.9802884083965534e-01 7.0541446618110337e-01 -7.9104505633399519e-02 + 6 2.0114733380305516e+00 -3.7433607398217483e+00 -2.9688665724203247e+00 + 7 -3.2646062589844127e-01 7.4545113602886504e-01 3.8864515374233779e+00 + 8 -1.6992144552176487e+00 4.1770867949498349e+00 2.4329479646797592e+00 + 9 1.8250302102601557e+00 -5.6775183441065051e+00 1.6562334058459376e+00 + 10 -2.9411499784567297e-01 4.2483825717579743e-01 -2.0454944575981404e-01 + 11 -9.6852892734814910e-01 1.2825160904684365e+00 -4.0863298208468957e-01 + 12 3.3107072705734946e+00 -7.6411402305007570e-01 1.2747076834264663e+00 + 13 -3.0999843695510060e-01 2.5003385487095964e-01 1.0600468871727173e-01 + 14 -1.2638813019667337e+00 5.9487350384637094e-01 2.0105698487957110e-01 + 15 3.4335793714273627e-01 -1.5950553742170387e-02 -1.0791712067698431e+00 + 16 -1.0530808398954454e+00 -3.6657367999349177e-01 2.7960862384014673e+00 + 17 -2.4285042656947571e+00 5.1576639561985473e+00 -6.9782029543394053e+00 + 18 -3.7929628350149169e-01 1.0420786070107907e+01 -9.5123267863543006e+00 + 19 1.8527235581340016e+00 -1.4163845910512083e+00 7.5667871079591533e+00 + 20 -3.2257007148929735e+00 -3.8484754900366887e+00 5.3176481256299590e+00 + 21 4.0318672433929734e+00 1.5871903248026678e+00 -1.0237421892763503e+01 + 22 1.4629760972442929e+00 -9.3386175185967235e-01 6.6032592166314519e+00 + 23 -4.4717803858275476e+00 -3.8939153227171523e+00 4.1111064353488365e+00 + 24 -1.5477331915053858e+00 1.1733740828713346e+01 -6.0432763590713154e+00 + 25 3.2290132418876336e+00 -4.1305242042670889e+00 3.9442459043882159e+00 + 26 -2.2770020178904948e+00 -8.5834251121651697e+00 -1.1704818988226089e+00 + 27 -3.7015228676694329e+00 1.0459984537547431e+01 -3.7654963783953805e+00 + 28 5.4222452431134922e+00 -4.4818388539753489e+00 3.6800169942739531e+00 + 29 -1.2945511546367356e+00 -6.0823641023924875e+00 5.8148360356210294e-02 ... diff --git a/unittest/force-styles/tests/mol-pair-coul_cut_global.yaml b/unittest/force-styles/tests/mol-pair-coul_cut_global.yaml new file mode 100644 index 0000000000..817c7c76ca --- /dev/null +++ b/unittest/force-styles/tests/mol-pair-coul_cut_global.yaml @@ -0,0 +1,85 @@ +--- +lammps_version: 10 Feb 2021 +date_generated: Fri Feb 26 23:08:42 2021 +epsilon: 1e-13 +prerequisites: ! | + atom full + pair coul/cut/global +pre_commands: ! "" +post_commands: ! "" +input_file: in.fourmol +pair_style: coul/cut/global 8.0 +pair_coeff: ! | + * * +extract: ! | + cut_coul 0 +natoms: 29 +init_vdwl: 0 +init_coul: -127.494586297384 +init_stress: ! |- + -2.6824547034957277e+01 -4.3162775104089022e+01 -5.7507264158338124e+01 -8.2055093564543602e-01 -7.9072701063929465e+00 3.0552092014126533e+00 +init_forces: ! |2 + 1 2.2407335289699457e+00 -5.8916421262140251e-02 4.2668639853362533e-01 + 2 3.0979928510686561e-01 -2.1737695778193569e+00 -1.8762407787145119e+00 + 3 -1.5509825092714681e-02 -8.6264917303471161e-02 3.0367577692877928e-02 + 4 5.1555686582278520e-02 1.0144901459000669e-02 -3.0866112200648410e-01 + 5 -3.9871973553557233e-01 7.0523075214913145e-01 -7.8022318023030612e-02 + 6 2.0159901805654243e+00 -3.7483112004912753e+00 -2.9733937038705189e+00 + 7 -3.2885029720170206e-01 7.5012396443748963e-01 3.8958946746344409e+00 + 8 -1.7087898110045159e+00 4.1818078578969917e+00 2.4290227133489601e+00 + 9 1.8279010754687195e+00 -5.6724899047236770e+00 1.6512690951588267e+00 + 10 -2.9363989744075952e-01 4.2512276557131090e-01 -2.0317384271111175e-01 + 11 -9.8868358921602750e-01 1.1284880357946445e+00 -5.3995772876854398e-01 + 12 3.3104189805348057e+00 -7.6522029601079322e-01 1.2706541181424338e+00 + 13 -3.8231563183336364e-01 -1.1565866576649098e-02 -8.7140102078706039e-03 + 14 -1.2285578202646528e+00 4.4726948380599840e-01 -1.5341285475604788e-02 + 15 2.2085158688210516e-01 -1.6336214445566244e-01 -9.6577522905262758e-01 + 16 -1.0515762083518698e+00 -3.6728607969167482e-01 2.7987602083946292e+00 + 17 -2.3442527695736954e+00 6.0494781225943264e+00 -7.7669898420813883e+00 + 18 1.2747707302049942e+00 6.6751091235431907e+00 -9.9139029454048035e+00 + 19 1.6181330640335625e+00 -1.6591712461142003e+00 5.6543694708933723e+00 + 20 -3.4516847879389365e+00 -3.1783992642977736e+00 4.2589367844805865e+00 + 21 2.4909407783496178e+00 1.4847928289095484e+00 -1.0503546063193097e+01 + 22 1.4973505700836185e+00 1.0712683327319574e+00 6.7413935914321712e+00 + 23 -4.4865766347042468e+00 -3.8185091520236116e+00 4.2817682561025867e+00 + 24 -2.4197806588078703e+00 9.8687406957238029e+00 -2.3585711007230041e+00 + 25 3.5804022684917332e+00 -3.1080767536674005e+00 3.9458463586096970e+00 + 26 -2.0751789623087848e+00 -7.6352193307876766e+00 -4.5052467808202223e-01 + 27 -2.8475395052833852e+00 1.1111959093179417e+01 -4.4252374188335644e+00 + 28 4.8464204075747359e+00 -5.2910383050538492e+00 4.0808901405648896e+00 + 29 -1.2636120082903082e+00 -6.1719354975176035e+00 9.2219267915908476e-01 +run_vdwl: 0 +run_coul: -127.551092913469 +run_stress: ! |- + -2.6928416018548067e+01 -4.3177725729834854e+01 -5.7444951165085541e+01 -8.5814323779813129e-01 -7.9912081348585087e+00 2.9897259857467389e+00 +run_forces: ! |2 + 1 2.2424368291189611e+00 -5.7157845887417236e-02 4.3210281185121069e-01 + 2 3.0183044940246062e-01 -2.1770211242302246e+00 -1.8796316672405435e+00 + 3 -1.5579134024760494e-02 -8.6577075665555045e-02 3.0374580927846943e-02 + 4 5.2605597845996610e-02 1.0081774332824270e-02 -3.0969764881211892e-01 + 5 -3.9802063955902828e-01 7.0541551998225649e-01 -7.9105180426906119e-02 + 6 2.0115016302231399e+00 -3.7433268305651115e+00 -2.9688564502920158e+00 + 7 -3.2650169036425236e-01 7.4540116984953042e-01 3.8864116015671808e+00 + 8 -1.6992053224967520e+00 4.1770597030834642e+00 2.4329396564406349e+00 + 9 1.8250267620013743e+00 -5.6775189439742455e+00 1.6562371734234631e+00 + 10 -2.9411864176775349e-01 4.2484150168808632e-01 -2.0454722162066086e-01 + 11 -9.8864600485826815e-01 1.1292383541031983e+00 -5.3894807546123469e-01 + 12 3.3107287761242974e+00 -7.6414239662660099e-01 1.2747075215620911e+00 + 13 -3.8336813011366144e-01 -1.2486400314578122e-02 -9.4718485253563467e-03 + 14 -1.2290837060322926e+00 4.4862678961891689e-01 -1.5069551071827086e-02 + 15 2.2244012017801007e-01 -1.6546887582342559e-01 -9.6944707732754809e-01 + 16 -1.0531155341397860e+00 -3.6654926990083037e-01 2.7960996028622298e+00 + 17 -2.3452245945824948e+00 6.0586966641613316e+00 -7.7701709031633026e+00 + 18 1.2287925114358347e+00 6.6283596097346500e+00 -9.8694991818018014e+00 + 19 1.6548495610099294e+00 -1.6299328734082061e+00 5.6681416753764280e+00 + 20 -3.4405939298489487e+00 -3.1595057592530829e+00 4.2031829943411507e+00 + 21 2.4999236613139768e+00 1.4581441593082431e+00 -1.0498791509490704e+01 + 22 1.5230738475837411e+00 1.0948611353935194e+00 6.7461577757855986e+00 + 23 -4.5208030817248321e+00 -3.8176488387685179e+00 4.2715025892256220e+00 + 24 -2.4424910840260479e+00 9.8889784537097576e+00 -2.3784455147561552e+00 + 25 3.6199382208005479e+00 -3.1023007862101899e+00 3.9803408068580102e+00 + 26 -2.0901080780170158e+00 -7.6586474495008154e+00 -4.6300658746615819e-01 + 27 -2.8661140489132810e+00 1.1127847846706011e+01 -4.4084385064798797e+00 + 28 4.8657388732889295e+00 -5.2969930881855793e+00 4.0790567237989928e+00 + 29 -1.2659132198580254e+00 -6.1822751233574085e+00 9.0587140991575621e-01 +... From 0c2fc07cc5462e28cd228028aea4be3070c981d0 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 12 Apr 2021 10:00:52 -0400 Subject: [PATCH 13/14] add USER-OMP version of pair style coul/cut/global --- src/USER-OMP/pair_coul_cut_global_omp.cpp | 44 ++++++++++++++++++ src/USER-OMP/pair_coul_cut_global_omp.h | 55 +++++++++++++++++++++++ 2 files changed, 99 insertions(+) create mode 100644 src/USER-OMP/pair_coul_cut_global_omp.cpp create mode 100644 src/USER-OMP/pair_coul_cut_global_omp.h diff --git a/src/USER-OMP/pair_coul_cut_global_omp.cpp b/src/USER-OMP/pair_coul_cut_global_omp.cpp new file mode 100644 index 0000000000..36fee39fc8 --- /dev/null +++ b/src/USER-OMP/pair_coul_cut_global_omp.cpp @@ -0,0 +1,44 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://lammps.sandia.gov/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "pair_coul_cut_global_omp.h" + +#include "error.h" + +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairCoulCutGlobalOMP::coeff(int narg, char **arg) +{ + if (narg != 2) + error->all(FLERR,"Incorrect args for pair coefficients"); + + PairCoulCut::coeff(narg,arg); +} + + +/* ---------------------------------------------------------------------- */ + +void *PairCoulCutGlobalOMP::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"cut_coul") == 0) return (void *) &cut_global; + dim = 2; + if (strcmp(str,"scale") == 0) return (void *) scale; + return nullptr; +} diff --git a/src/USER-OMP/pair_coul_cut_global_omp.h b/src/USER-OMP/pair_coul_cut_global_omp.h new file mode 100644 index 0000000000..a2cb4a1dbe --- /dev/null +++ b/src/USER-OMP/pair_coul_cut_global_omp.h @@ -0,0 +1,55 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(coul/cut/global/omp,PairCoulCutGlobalOMP) + +#else + +#ifndef LMP_PAIR_COUL_CUT_GLOBAL_OMP_H +#define LMP_PAIR_COUL_CUT_GLOBAL_OMP_H + +#include "pair_coul_cut_omp.h" + +namespace LAMMPS_NS { + +class PairCoulCutGlobalOMP : public PairCoulCutOMP { + public: + PairCoulCutGlobalOMP(class LAMMPS *lmp) : PairCoulCutOMP(lmp) {} + void coeff(int, char **); + void *extract(const char *, int &); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair style coul/cut requires atom attribute q + +The atom style defined does not have these attributes. + +*/ From 3925bcc1dee5be618f25b293ad234ced6e5e5366 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 12 Apr 2021 13:35:42 -0400 Subject: [PATCH 14/14] apply requested changes do pair style hybrid doc --- doc/src/pair_hybrid.rst | 60 ++++++++++++++++++++++------------------- 1 file changed, 33 insertions(+), 27 deletions(-) diff --git a/doc/src/pair_hybrid.rst b/doc/src/pair_hybrid.rst index 9fdba318d9..e4e253caf9 100644 --- a/doc/src/pair_hybrid.rst +++ b/doc/src/pair_hybrid.rst @@ -14,7 +14,8 @@ pair_style hybrid/overlay command Accelerator Variants: *hybrid/overlay/kk* -pair_style hybrid/scale command +pair_style hybrid/scaled command +================================== Syntax """""" @@ -42,6 +43,10 @@ Examples pair_coeff * * lj/cut 1.0 1.0 pair_coeff * * coul/long + pair_style hybrid/scaled 0.5 tersoff 0.5 sw + pair_coeff * * tersoff Si.tersoff Si + pair_coeff * * sw Si.sw Si + variable one equal ramp(1.0,0.0) variable two equal 1.0-v_one pair_style hybrid/scaled v_one lj/cut 2.5 v_two morse 2.5 @@ -57,11 +62,11 @@ 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. The assignment of pair styles to type pairs is made via the :doc:`pair_coeff ` command. -The *hybrid/scaled* style differs from the *hybrid/overlay* style by -requiring a factor for each pair style that is used to scale all -forces, energies and stresses computed by each sub-style. Because of -the additional complexity, the *hybrid/scaled* style will have more -overhead and thus will be a bit slower than *hybrid/overlay*. +The major difference between the *hybrid/overlay* and *hybrid/scaled* +styles is that the *hybrid/scaled* adds a scale factor for each +sub-style contribution to forces, energies and stresses. Because of the +added complexity, the *hybrid/scaled* style has more overhead and thus +may be slower than *hybrid/overlay*. Here are two examples of hybrid simulations. The *hybrid* style could be used for a simulation of a metal droplet on a LJ surface. The metal @@ -76,34 +81,35 @@ and *coul/long* together gives the same result as if the would be more efficient to use the single combined potential, but in general any combination of pair potentials can be used together in to produce an interaction that is not encoded in any single pair_style -file, e.g. adding Coulombic forces between granular particles. The -*hybrid/scaled* style enables more complex combinations of pair styles -than a simple sum as *hybrid/overlay* does; there may be fractional -contributions from sub-styles or contributions may be subtracted with a -negative scale factor. Furthermore, since the scale factors can be -variables that may change during a simulation, which would allow, for -instance, to smoothly switch between two different pair styles or two -different parameter sets. +file, e.g. adding Coulombic forces between granular particles. + +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 +switching smoothly between two different pair styles or two different +parameter sets during a run. 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 -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 individual pair styles for a listing and explanation of the -appropriate arguments. +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 +appropriate arguments for them. Note that an individual pair style can be used multiple times as a -sub-style. For efficiency this should only be done if your model -requires it. E.g. if you have different regions of Si and C atoms and -wish to use a Tersoff potential for pure Si for one set of atoms, and -a Tersoff potential for pure C for the other set (presumably with some -third potential for Si-C interactions), then the sub-style *tersoff* -could be listed twice. But if you just want to use a Lennard-Jones or -other pairwise potential for several different atom type pairs in your -model, then you should just list the sub-style once and use the -pair_coeff command to assign parameters for the different type pairs. +sub-style. For efficiency reasons this should only be done if your +model requires it. E.g. if you have different regions of Si and C atoms +and wish to use a Tersoff potential for pure Si for one set of atoms, +and a Tersoff potential for pure C for the other set (presumably with +some third potential for Si-C interactions), then the sub-style +*tersoff* could be listed twice. But if you just want to use a +Lennard-Jones or other pairwise potential for several different atom +type pairs in your model, then you should just list the sub-style once +and use the pair_coeff command to assign parameters for the different +type pairs. .. note::