From 8ed6d80b851a4dbe02aa61bcd50ca12e5d3c15b1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 8 Apr 2021 21:02:28 -0400 Subject: [PATCH] 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(); }; }