From 3f3c481554ee190fe352a4566bb228445b55ed7f Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Thu, 21 Apr 2022 17:00:11 -0600 Subject: [PATCH 1/6] add support to fix adapt for angle coeffs --- doc/src/fix_adapt.rst | 81 +++++++++++++++------ src/EXTRA-MOLECULE/bond_fene_nm.cpp | 2 +- src/EXTRA-MOLECULE/bond_gaussian.cpp | 4 +- src/MOLECULE/angle_harmonic.cpp | 12 ++++ src/MOLECULE/angle_harmonic.h | 1 + src/MOLECULE/bond_fene.cpp | 2 +- src/MOLECULE/bond_gromos.cpp | 7 +- src/MOLECULE/bond_harmonic.cpp | 10 ++- src/angle.cpp | 12 ++++ src/angle.h | 5 ++ src/bond.cpp | 7 +- src/bond.h | 6 +- src/fix_adapt.cpp | 102 ++++++++++++++++++++++++--- src/fix_adapt.h | 8 ++- 14 files changed, 202 insertions(+), 57 deletions(-) diff --git a/doc/src/fix_adapt.rst b/doc/src/fix_adapt.rst index b7151334db..0f68219358 100644 --- a/doc/src/fix_adapt.rst +++ b/doc/src/fix_adapt.rst @@ -14,7 +14,7 @@ Syntax * adapt = style name of this fix command * N = adapt simulation settings every this many timesteps * one or more attribute/arg pairs may be appended -* attribute = *pair* or *bond* or *kspace* or *atom* +* attribute = *pair* or *bond* or *angle* or *kspace* or *atom* .. parsed-literal:: @@ -28,11 +28,16 @@ Syntax bparam = parameter to adapt over time I = type bond to set parameter for v_name = variable with name that calculates value of bparam + *angle* args = astyle aparam I v_name + astyle = angle style name, e.g. harmonic + aparam = parameter to adapt over time + I = type angle to set parameter for + v_name = variable with name that calculates value of aparam *kspace* arg = v_name v_name = variable with name that calculates scale factor on K-space terms - *atom* args = aparam v_name - aparam = parameter to adapt over time - v_name = variable with name that calculates value of aparam + *atom* args = atomparam v_name + atomparam = parameter to adapt over time + v_name = variable with name that calculates value of atomparam * zero or more keyword/value pairs may be appended * keyword = *scale* or *reset* or *mass* @@ -283,30 +288,60 @@ operates. The only difference is that now a bond coefficient for a given bond type is adapted. A wild-card asterisk can be used in place of or in conjunction with -the bond type argument to set the coefficients for multiple bond types. -This takes the form "\*" or "\*n" or "n\*" or "m\*n". If N = the number of -atom types, then an asterisk with no numeric values means all types -from 1 to N. A leading asterisk means all types from 1 to n (inclusive). -A trailing asterisk means all types from n to N (inclusive). A middle -asterisk means all types from m to n (inclusive). +the bond type argument to set the coefficients for multiple bond +types. This takes the form "\*" or "\*n" or "n\*" or "m\*n". If N = +the number of bond types, then an asterisk with no numeric values +means all types from 1 to N. A leading asterisk means all types from +1 to n (inclusive). A trailing asterisk means all types from n to N +(inclusive). A middle asterisk means all types from m to n +(inclusive). Currently *bond* does not support bond_style hybrid nor bond_style hybrid/overlay as bond styles. The only bonds that currently are working with fix_adapt are -+------------------------------------+-------+------------+ -| :doc:`class2 ` | r0 | type bonds | -+------------------------------------+-------+------------+ -| :doc:`fene ` | k, r0 | type bonds | -+------------------------------------+-------+------------+ -| :doc:`gromos ` | k, r0 | type bonds | -+------------------------------------+-------+------------+ -| :doc:`harmonic ` | k,r0 | type bonds | -+------------------------------------+-------+------------+ -| :doc:`morse ` | r0 | type bonds | -+------------------------------------+-------+------------+ -| :doc:`nonlinear ` | r0 | type bonds | -+------------------------------------+-------+------------+ ++------------------------------------+-------+-----------------+ +| :doc:`class2 ` | r0 | type bonds | ++------------------------------------+-------+-----------------+ +| :doc:`fene ` | k, r0 | type bonds | ++------------------------------------+-------+-----------------+ +| :doc:`fene/nm ` | k, r0 | type bonds | ++------------------------------------+-------+-----------------+ +| :doc:`gromos ` | k, r0 | type bonds | ++------------------------------------+-------+-----------------+ +| :doc:`harmonic ` | k,r0 | type bonds | ++------------------------------------+-------+-----------------+ +| :doc:`morse ` | r0 | type bonds | ++------------------------------------+-------+-----------------+ +| :doc:`nonlinear ` | epsilon,r0 | type bonds | ++------------------------------------+-------+-----------------+ + +---------- + +The *angle* keyword uses the specified variable to change the value of +an angle coefficient over time, very similar to how the *pair* keyword +operates. The only difference is that now an angle coefficient for a +given angle type is adapted. + +A wild-card asterisk can be used in place of or in conjunction with +the angle type argument to set the coefficients for multiple angle +types. This takes the form "\*" or "\*n" or "n\*" or "m\*n". If N = +the number of angle types, then an asterisk with no numeric values +means all types from 1 to N. A leading asterisk means all types from +1 to n (inclusive). A trailing asterisk means all types from n to N +(inclusive). A middle asterisk means all types from m to n +(inclusive). + +Currently *angle* does not support angle_style hybrid nor angle_style +hybrid/overlay as angle styles. The only angles that currently are +working with fix_adapt are + ++------------------------------------+-------+-----------------+ +| :doc:`harmonic ` | k,theta0 | type angles | ++------------------------------------+-------+-----------------+ + +Note that internally, theta0 is stored in radians, so the variable +this fix uses to reset theta0 needs to generate values in radians. ---------- diff --git a/src/EXTRA-MOLECULE/bond_fene_nm.cpp b/src/EXTRA-MOLECULE/bond_fene_nm.cpp index 839217df78..bd9d16b95d 100644 --- a/src/EXTRA-MOLECULE/bond_fene_nm.cpp +++ b/src/EXTRA-MOLECULE/bond_fene_nm.cpp @@ -273,7 +273,7 @@ double BondFENENM::single(int type, double rsq, int /*i*/, int /*j*/, double &ff void *BondFENENM::extract(const char *str, int &dim) { dim = 1; - if (strcmp(str, "kappa") == 0) return (void *) k; + if (strcmp(str, "k") == 0) return (void *) k; if (strcmp(str, "r0") == 0) return (void *) r0; return nullptr; } diff --git a/src/EXTRA-MOLECULE/bond_gaussian.cpp b/src/EXTRA-MOLECULE/bond_gaussian.cpp index c2ab00dfde..0ad125aeca 100644 --- a/src/EXTRA-MOLECULE/bond_gaussian.cpp +++ b/src/EXTRA-MOLECULE/bond_gaussian.cpp @@ -34,9 +34,7 @@ using namespace MathConst; BondGaussian::BondGaussian(LAMMPS *lmp) : Bond(lmp), nterms(nullptr), bond_temperature(nullptr), alpha(nullptr), width(nullptr), r0(nullptr) -{ - reinitflag = 1; -} +{} /* ---------------------------------------------------------------------- */ diff --git a/src/MOLECULE/angle_harmonic.cpp b/src/MOLECULE/angle_harmonic.cpp index 261770d9a3..aa24fa27b4 100644 --- a/src/MOLECULE/angle_harmonic.cpp +++ b/src/MOLECULE/angle_harmonic.cpp @@ -264,3 +264,15 @@ double AngleHarmonic::single(int type, int i1, int i2, int i3) double tk = k[type] * dtheta; return tk * dtheta; } + +/* ---------------------------------------------------------------------- + return ptr to internal members upon request +------------------------------------------------------------------------ */ + +void *AngleHarmonic::extract(const char *str, int &dim) +{ + dim = 1; + if (strcmp(str, "k") == 0) return (void *) k; + if (strcmp(str, "theta0") == 0) return (void *) theta0; + return nullptr; +} diff --git a/src/MOLECULE/angle_harmonic.h b/src/MOLECULE/angle_harmonic.h index 718ac4bd0a..2643ea3a17 100644 --- a/src/MOLECULE/angle_harmonic.h +++ b/src/MOLECULE/angle_harmonic.h @@ -35,6 +35,7 @@ class AngleHarmonic : public Angle { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, int, int, int) override; + void *extract(const char *, int &) override; protected: double *k, *theta0; diff --git a/src/MOLECULE/bond_fene.cpp b/src/MOLECULE/bond_fene.cpp index 8daf6e862c..2950d48ca6 100644 --- a/src/MOLECULE/bond_fene.cpp +++ b/src/MOLECULE/bond_fene.cpp @@ -265,7 +265,7 @@ double BondFENE::single(int type, double rsq, int /*i*/, int /*j*/, double &ffor void *BondFENE::extract(const char *str, int &dim) { dim = 1; - if (strcmp(str, "kappa") == 0) return (void *) k; + if (strcmp(str, "k") == 0) return (void *) k; if (strcmp(str, "r0") == 0) return (void *) r0; return nullptr; } diff --git a/src/MOLECULE/bond_gromos.cpp b/src/MOLECULE/bond_gromos.cpp index adf3f91044..4572d2c900 100644 --- a/src/MOLECULE/bond_gromos.cpp +++ b/src/MOLECULE/bond_gromos.cpp @@ -30,10 +30,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -BondGromos::BondGromos(LAMMPS *_lmp) : Bond(_lmp) -{ - reinitflag = 1; -} +BondGromos::BondGromos(LAMMPS *_lmp) : Bond(_lmp) {} /* ---------------------------------------------------------------------- */ @@ -200,7 +197,7 @@ double BondGromos::single(int type, double rsq, int /*i*/, int /*j*/, double &ff void *BondGromos::extract(const char *str, int &dim) { dim = 1; - if (strcmp(str, "kappa") == 0) return (void *) k; + if (strcmp(str, "k") == 0) return (void *) k; if (strcmp(str, "r0") == 0) return (void *) r0; return nullptr; } diff --git a/src/MOLECULE/bond_harmonic.cpp b/src/MOLECULE/bond_harmonic.cpp index cff766aa3b..687a871f17 100644 --- a/src/MOLECULE/bond_harmonic.cpp +++ b/src/MOLECULE/bond_harmonic.cpp @@ -27,10 +27,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -BondHarmonic::BondHarmonic(LAMMPS *_lmp) : Bond(_lmp) -{ - reinitflag = 1; -} +BondHarmonic::BondHarmonic(LAMMPS *_lmp) : Bond(_lmp) {} /* ---------------------------------------------------------------------- */ @@ -201,12 +198,13 @@ double BondHarmonic::single(int type, double rsq, int /*i*/, int /*j*/, double & } /* ---------------------------------------------------------------------- - Return ptr to internal members upon request. + return ptr to internal members upon request ------------------------------------------------------------------------ */ + void *BondHarmonic::extract(const char *str, int &dim) { dim = 1; - if (strcmp(str, "kappa") == 0) return (void *) k; + if (strcmp(str, "k") == 0) return (void *) k; if (strcmp(str, "r0") == 0) return (void *) r0; return nullptr; } diff --git a/src/angle.cpp b/src/angle.cpp index 52d92b72b2..caa86d691e 100644 --- a/src/angle.cpp +++ b/src/angle.cpp @@ -353,3 +353,15 @@ double Angle::memory_usage() bytes += (double) comm->nthreads * maxcvatom * 9 * sizeof(double); return bytes; } + +/* ----------------------------------------------------------------------- + reset all type-based angle params via init() +-------------------------------------------------------------------------- */ + +void Angle::reinit() +{ + if (!reinitflag) + error->all(FLERR, "Fix adapt interface to this angle style not supported"); + + init(); +} diff --git a/src/angle.h b/src/angle.h index 12443fa4f3..8e200ce37b 100644 --- a/src/angle.h +++ b/src/angle.h @@ -36,6 +36,9 @@ class Angle : protected Pointers { // CENTROID_AVAIL = different and implemented // CENTROID_NOTAVAIL = different, not yet implemented + int reinitflag; // 0 if not compatible with fix adapt + // extract() method may still need to be added + // KOKKOS host/device flag and data masks ExecutionSpace execution_space; @@ -57,6 +60,8 @@ class Angle : protected Pointers { virtual void write_data(FILE *) {} virtual double single(int, int, int, int) = 0; virtual double memory_usage(); + virtual void *extract(const char *, int &) { return nullptr; } + void reinit(); protected: int suffix_flag; // suffix compatibility flag diff --git a/src/bond.cpp b/src/bond.cpp index 6ae78bc429..7f7b6aa3a3 100644 --- a/src/bond.cpp +++ b/src/bond.cpp @@ -43,6 +43,7 @@ Bond::Bond(LAMMPS *_lmp) : Pointers(_lmp) energy = 0.0; virial[0] = virial[1] = virial[2] = virial[3] = virial[4] = virial[5] = 0.0; writedata = 1; + reinitflag = 1; comm_forward = comm_reverse = comm_reverse_off = 0; @@ -336,11 +337,13 @@ double Bond::memory_usage() } /* ----------------------------------------------------------------------- - Reset all type-based bond params via init. + reset all type-based bond params via init() -------------------------------------------------------------------------- */ + void Bond::reinit() { - if (!reinitflag) error->all(FLERR, "Fix adapt interface to this bond style not supported"); + if (!reinitflag) + error->all(FLERR, "Fix adapt interface to this bond style not supported"); init(); } diff --git a/src/bond.h b/src/bond.h index 47c8687f50..7b6a70dc7d 100644 --- a/src/bond.h +++ b/src/bond.h @@ -37,7 +37,8 @@ class Bond : protected Pointers { int comm_reverse; // size of reverse communication (0 if none) int comm_reverse_off; // size of reverse comm even if newton off - int reinitflag; // 1 if compatible with fix adapt and alike + int reinitflag; // 0 if not compatible with fix adapt + // extract() method may still need to be added // KOKKOS host/device flag and data masks @@ -61,7 +62,8 @@ class Bond : protected Pointers { virtual double single(int, double, int, int, double &) = 0; virtual double memory_usage(); virtual void *extract(const char *, int &) { return nullptr; } - virtual void reinit(); + void reinit(); + virtual int pack_forward_comm(int, int *, double *, int, int *) {return 0;} virtual void unpack_forward_comm(int, int, double *) {} virtual int pack_reverse_comm(int, int, double *) {return 0;} diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index 2632ccf597..0d9aa00b7a 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -14,6 +14,7 @@ #include "fix_adapt.h" +#include "angle.h" #include "atom.h" #include "bond.h" #include "domain.h" @@ -38,7 +39,7 @@ using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; -enum{PAIR,KSPACE,ATOM,BOND}; +enum{PAIR,KSPACE,ATOM,BOND,ANGLE}; enum{DIAMETER,CHARGE}; /* ---------------------------------------------------------------------- */ @@ -75,6 +76,10 @@ nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) if (iarg+5 > narg) error->all(FLERR,"Illegal fix adapt command"); nadapt++; iarg += 5; + } else if (strcmp(arg[iarg],"angle") == 0) { + if (iarg+5 > narg) error->all(FLERR,"Illegal fix adapt command"); + nadapt++; + iarg += 5; } else break; } @@ -119,6 +124,20 @@ nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) nadapt++; iarg += 5; + } else if (strcmp(arg[iarg],"angle") == 0) { + if (iarg+5 > narg) error->all(FLERR, "Illegal fix adapt command"); + adapt[nadapt].which = ANGLE; + adapt[nadapt].angle = nullptr; + adapt[nadapt].astyle = utils::strdup(arg[iarg+1]); + adapt[nadapt].aparam = utils::strdup(arg[iarg+2]); + utils::bounds(FLERR,arg[iarg+3],1,atom->nangletypes, + adapt[nadapt].ilo,adapt[nadapt].ihi,error); + if (utils::strmatch(arg[iarg+4],"^v_")) { + adapt[nadapt].var = utils::strdup(arg[iarg+4]+2); + } else error->all(FLERR,"Illegal fix adapt command"); + nadapt++; + iarg += 5; + } else if (strcmp(arg[iarg],"kspace") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); adapt[nadapt].which = KSPACE; @@ -133,12 +152,12 @@ nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) adapt[nadapt].which = ATOM; if (strcmp(arg[iarg+1],"diameter") == 0 || strcmp(arg[iarg+1],"diameter/disc") == 0) { - adapt[nadapt].aparam = DIAMETER; + adapt[nadapt].atomparam = DIAMETER; diamflag = 1; discflag = 0; if (strcmp(arg[iarg+1],"diameter/disc") == 0) discflag = 1; } else if (strcmp(arg[iarg+1],"charge") == 0) { - adapt[nadapt].aparam = CHARGE; + adapt[nadapt].atomparam = CHARGE; chgflag = 1; } else error->all(FLERR,"Illegal fix adapt command"); if (utils::strmatch(arg[iarg+2],"^v_")) { @@ -191,6 +210,13 @@ nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) for (int m = 0; m < nadapt; ++m) if (adapt[m].which == BOND) memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig"); + + // allocate angle style arrays: + + n = atom->nbondtypes; + for (int m = 0; m < nadapt; ++m) + if (adapt[m].which == ANGLE) + memory->create(adapt[m].vector_orig,n+1,"adapt:vector_orig"); } /* ---------------------------------------------------------------------- */ @@ -207,6 +233,10 @@ FixAdapt::~FixAdapt() delete [] adapt[m].bstyle; delete [] adapt[m].bparam; memory->destroy(adapt[m].vector_orig); + } else if (adapt[m].which == ANGLE) { + delete [] adapt[m].astyle; + delete [] adapt[m].aparam; + memory->destroy(adapt[m].vector_orig); } } delete [] adapt; @@ -357,6 +387,7 @@ void FixAdapt::init() } delete [] pstyle; + } else if (ad->which == BOND) { ad->bond = nullptr; anybond = 1; @@ -383,13 +414,39 @@ void FixAdapt::init() delete [] bstyle; + } else if (ad->which == ANGLE) { + ad->angle = nullptr; + anyangle = 1; + + char *astyle = utils::strdup(ad->astyle); + if (lmp->suffix_enable) + ad->angle = force->angle_match(fmt::format("{}/{}",astyle,lmp->suffix)); + + if (ad->angle == nullptr) ad->angle = force->angle_match(astyle); + if (ad->angle == nullptr ) + error->all(FLERR,"Fix adapt angle style does not exist"); + + void *ptr = ad->angle->extract(ad->aparam,ad->bdim); + + if (ptr == nullptr) + error->all(FLERR,"Fix adapt angle style param not supported"); + + // for angle styles, use a vector + + if (ad->adim == 1) ad->vector = (double *) ptr; + + if (utils::strmatch(force->angle_style,"^hybrid")) + error->all(FLERR,"Fix adapt does not support angle_style hybrid"); + + delete [] astyle; + } else if (ad->which == KSPACE) { if (force->kspace == nullptr) error->all(FLERR,"Fix adapt kspace style does not exist"); kspace_scale = (double *) force->kspace->extract("scale"); } else if (ad->which == ATOM) { - if (ad->aparam == DIAMETER) { + if (ad->atomparam == DIAMETER) { if (!atom->radius_flag) error->all(FLERR,"Fix adapt requires atom attribute diameter"); if (!atom->rmass_flag) @@ -398,7 +455,7 @@ void FixAdapt::init() error->all(FLERR,"Fix adapt requires 2d simulation"); if (!restart_reset) previous_diam_scale = 1.0; } - if (ad->aparam == CHARGE) { + if (ad->atomparam == CHARGE) { if (!atom->q_flag) error->all(FLERR,"Fix adapt requires atom attribute charge"); if (!restart_reset) previous_chg_scale = 1.0; @@ -408,7 +465,7 @@ void FixAdapt::init() if (restart_reset) restart_reset = 0; - // make copy of original pair/bond array values + // make copy of original pair/bond/angle array values for (int m = 0; m < nadapt; m++) { Adapt *ad = &adapt[m]; @@ -422,6 +479,10 @@ void FixAdapt::init() } else if (ad->which == BOND && ad->bdim == 1) { for (i = ad->ilo; i <= ad->ihi; ++i ) ad->vector_orig[i] = ad->vector[i]; + + } else if (ad->which == ANGLE && ad->adim == 1) { + for (i = ad->ilo; i <= ad->ihi; ++i ) + ad->vector_orig[i] = ad->vector[i]; } } @@ -483,7 +544,7 @@ void FixAdapt::post_run() } /* ---------------------------------------------------------------------- - change pair,kspace,atom parameters based on variable evaluation + change pair,bond,angle,kspace,atom parameters based on variable evaluation ------------------------------------------------------------------------- */ void FixAdapt::change_settings() @@ -527,6 +588,18 @@ void FixAdapt::change_settings() ad->vector[i] = value; } + // set angle type array values: + + } else if (ad->which == ANGLE) { + if (ad->adim == 1) { + if (scaleflag) + for (i = ad->ilo; i <= ad->ihi; ++i ) + ad->vector[i] = value*ad->vector_orig[i]; + else + for (i = ad->ilo; i <= ad->ihi; ++i ) + ad->vector[i] = value; + } + // set kspace scale factor } else if (ad->which == KSPACE) { @@ -540,7 +613,7 @@ void FixAdapt::change_settings() // also reset rmass to new value assuming density remains constant // for scaleflag, previous_diam_scale is the scale factor on previous step - if (ad->aparam == DIAMETER) { + if (ad->atomparam == DIAMETER) { double scale; double *radius = atom->radius; double *rmass = atom->rmass; @@ -567,7 +640,7 @@ void FixAdapt::change_settings() // reset charge to new value, for both owned and ghost atoms // for scaleflag, previous_chg_scale is the scale factor on previous step - } else if (ad->aparam == CHARGE) { + } else if (ad->atomparam == CHARGE) { double scale; double *q = atom->q; int *mask = atom->mask; @@ -591,7 +664,7 @@ void FixAdapt::change_settings() modify->addstep_compute(update->ntimestep + nevery); // re-initialize pair styles if any PAIR settings were changed - // ditto for bond styles if any BOND settings were changed + // ditto for bond/angle styles if any BOND/ANGLE settings were changed // this resets other coeffs that may depend on changed values, // and also offset and tail corrections // we must call force->pair->reinit() instead of the individual @@ -601,6 +674,7 @@ void FixAdapt::change_settings() if (anypair) force->pair->reinit(); if (anybond) force->bond->reinit(); + if (anyangle) force->angle->reinit(); // reset KSpace charges if charges have changed @@ -624,7 +698,13 @@ void FixAdapt::restore_settings() } } else if (ad->which == BOND) { - if (ad->pdim == 1) { + if (ad->bdim == 1) { + for (int i = ad->ilo; i <= ad->ihi; i++) + ad->vector[i] = ad->vector_orig[i]; + } + + } else if (ad->which == ANGLE) { + if (ad->adim == 1) { for (int i = ad->ilo; i <= ad->ihi; i++) ad->vector[i] = ad->vector_orig[i]; } diff --git a/src/fix_adapt.h b/src/fix_adapt.h index 939f46f8a2..121ef06ece 100644 --- a/src/fix_adapt.h +++ b/src/fix_adapt.h @@ -45,7 +45,7 @@ class FixAdapt : public Fix { private: int nadapt, resetflag, scaleflag, massflag; - int anypair, anybond; + int anypair, anybond, anyangle; int nlevels_respa; char *id_fix_diam, *id_fix_chg; class FixStore *fix_diam, *fix_chg; @@ -57,14 +57,16 @@ class FixAdapt : public Fix { char *var; char *pstyle, *pparam; char *bstyle, *bparam; + char *astyle, *aparam; int ilo, ihi, jlo, jhi; - int pdim, bdim; + int pdim, bdim, adim; double *scalar, scalar_orig; double *vector, *vector_orig; double **array, **array_orig; - int aparam; + int atomparam; class Pair *pair; class Bond *bond; + class Angle *angle; }; Adapt *adapt; From b4e2e2ec34a359079cae7bfe2d3d429cce7ad679 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 22 Apr 2022 09:12:12 -0600 Subject: [PATCH 2/6] add support for angle cosine --- doc/src/fix_adapt.rst | 16 +++++++++------- src/MOLECULE/angle_cosine.cpp | 11 +++++++++++ src/MOLECULE/angle_cosine.h | 1 + 3 files changed, 21 insertions(+), 7 deletions(-) diff --git a/doc/src/fix_adapt.rst b/doc/src/fix_adapt.rst index 0f68219358..1276adf444 100644 --- a/doc/src/fix_adapt.rst +++ b/doc/src/fix_adapt.rst @@ -297,17 +297,17 @@ means all types from 1 to N. A leading asterisk means all types from (inclusive). Currently *bond* does not support bond_style hybrid nor bond_style -hybrid/overlay as bond styles. The only bonds that currently are -working with fix_adapt are +hybrid/overlay as bond styles. The bond styles that currently work +with fix_adapt are +------------------------------------+-------+-----------------+ | :doc:`class2 ` | r0 | type bonds | +------------------------------------+-------+-----------------+ -| :doc:`fene ` | k, r0 | type bonds | +| :doc:`fene ` | k,r0 | type bonds | +------------------------------------+-------+-----------------+ -| :doc:`fene/nm ` | k, r0 | type bonds | +| :doc:`fene/nm ` | k,r0 | type bonds | +------------------------------------+-------+-----------------+ -| :doc:`gromos ` | k, r0 | type bonds | +| :doc:`gromos ` | k,r0 | type bonds | +------------------------------------+-------+-----------------+ | :doc:`harmonic ` | k,r0 | type bonds | +------------------------------------+-------+-----------------+ @@ -333,12 +333,14 @@ means all types from 1 to N. A leading asterisk means all types from (inclusive). Currently *angle* does not support angle_style hybrid nor angle_style -hybrid/overlay as angle styles. The only angles that currently are -working with fix_adapt are +hybrid/overlay as angle styles. The angle styles that currently work +with fix_adapt are +------------------------------------+-------+-----------------+ | :doc:`harmonic ` | k,theta0 | type angles | +------------------------------------+-------+-----------------+ +| :doc:`cosine ` | k | type angles | ++------------------------------------+-------+-----------------+ Note that internally, theta0 is stored in radians, so the variable this fix uses to reset theta0 needs to generate values in radians. diff --git a/src/MOLECULE/angle_cosine.cpp b/src/MOLECULE/angle_cosine.cpp index d32dc4559d..260ebc3948 100644 --- a/src/MOLECULE/angle_cosine.cpp +++ b/src/MOLECULE/angle_cosine.cpp @@ -234,3 +234,14 @@ double AngleCosine::single(int type, int i1, int i2, int i3) return k[type] * (1.0 + c); } + +/* ---------------------------------------------------------------------- + return ptr to internal members upon request +------------------------------------------------------------------------ */ + +void *AngleCosine::extract(const char *str, int &dim) +{ + dim = 1; + if (strcmp(str, "k") == 0) return (void *) k; + return nullptr; +} diff --git a/src/MOLECULE/angle_cosine.h b/src/MOLECULE/angle_cosine.h index 19b6222c87..e249e7a44c 100644 --- a/src/MOLECULE/angle_cosine.h +++ b/src/MOLECULE/angle_cosine.h @@ -35,6 +35,7 @@ class AngleCosine : public Angle { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, int, int, int) override; + void *extract(const char *, int &) override; protected: double *k; From c0d0c84f7d6735263a4560327f96e0e9a61e627b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 22 Apr 2022 12:08:23 -0400 Subject: [PATCH 3/6] update unit test files to implementation changes --- unittest/force-styles/tests/angle-cosine.yaml | 3 ++- unittest/force-styles/tests/angle-harmonic.yaml | 4 +++- unittest/force-styles/tests/bond-fene.yaml | 2 +- unittest/force-styles/tests/bond-fene_nm.yaml | 2 +- unittest/force-styles/tests/bond-gromos.yaml | 2 +- unittest/force-styles/tests/bond-harmonic.yaml | 2 +- 6 files changed, 9 insertions(+), 6 deletions(-) diff --git a/unittest/force-styles/tests/angle-cosine.yaml b/unittest/force-styles/tests/angle-cosine.yaml index 5a59fcce86..43629712d4 100644 --- a/unittest/force-styles/tests/angle-cosine.yaml +++ b/unittest/force-styles/tests/angle-cosine.yaml @@ -16,7 +16,8 @@ angle_coeff: ! | 3 50.0 4 100.0 equilibrium: 4 3.141592653589793 3.141592653589793 3.141592653589793 3.141592653589793 -extract: ! "" +extract: ! | + k 1 natoms: 29 init_energy: 1347.8670856939623 init_stress: ! |2- diff --git a/unittest/force-styles/tests/angle-harmonic.yaml b/unittest/force-styles/tests/angle-harmonic.yaml index dee700aa2c..d2164af8c5 100644 --- a/unittest/force-styles/tests/angle-harmonic.yaml +++ b/unittest/force-styles/tests/angle-harmonic.yaml @@ -16,7 +16,9 @@ angle_coeff: ! | 3 50.0 120.0 4 100.0 108.5 equilibrium: 4 1.9216075064457567 1.9373154697137058 2.0943951023931953 1.8936822384138476 -extract: ! "" +extract: ! | + k 1 + theta0 1 natoms: 29 init_energy: 41.53081789649104 init_stress: ! |2- diff --git a/unittest/force-styles/tests/bond-fene.yaml b/unittest/force-styles/tests/bond-fene.yaml index 88c02574cb..e5077eda0e 100644 --- a/unittest/force-styles/tests/bond-fene.yaml +++ b/unittest/force-styles/tests/bond-fene.yaml @@ -18,7 +18,7 @@ bond_coeff: ! | 5 450 2 0.018 1 equilibrium: 5 1.455 1.067 1.261 1.164 0.97 extract: ! | - kappa 1 + k 1 r0 1 natoms: 29 init_energy: 7104.900486467235 diff --git a/unittest/force-styles/tests/bond-fene_nm.yaml b/unittest/force-styles/tests/bond-fene_nm.yaml index c6be31a1c3..892d26c7aa 100644 --- a/unittest/force-styles/tests/bond-fene_nm.yaml +++ b/unittest/force-styles/tests/bond-fene_nm.yaml @@ -18,7 +18,7 @@ bond_coeff: ! | 5 450 2 0.018 1 12 6 equilibrium: 5 1.455 1.067 1.261 1.164 0.97 extract: ! | - kappa 1 + k 1 r0 1 natoms: 29 init_energy: 7104.538647187164 diff --git a/unittest/force-styles/tests/bond-gromos.yaml b/unittest/force-styles/tests/bond-gromos.yaml index a2f7e7ef3e..18abc99a3c 100644 --- a/unittest/force-styles/tests/bond-gromos.yaml +++ b/unittest/force-styles/tests/bond-gromos.yaml @@ -18,7 +18,7 @@ bond_coeff: ! | 5 450.0 1.0 equilibrium: 5 1.5 1.1 1.3 1.2 1 extract: ! | - kappa 1 + k 1 r0 1 natoms: 29 init_energy: 33.70930417641326 diff --git a/unittest/force-styles/tests/bond-harmonic.yaml b/unittest/force-styles/tests/bond-harmonic.yaml index 9ee14c07b9..bf686558d7 100644 --- a/unittest/force-styles/tests/bond-harmonic.yaml +++ b/unittest/force-styles/tests/bond-harmonic.yaml @@ -18,7 +18,7 @@ bond_coeff: ! | 5 450.0 1.0 equilibrium: 5 1.5 1.1 1.3 1.2 1 extract: ! | - kappa 1 + k 1 r0 1 natoms: 29 init_energy: 4.789374024601252 From c054edda6b2b4b317180ddb96d1485498095ac3e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 22 Apr 2022 13:22:01 -0400 Subject: [PATCH 4/6] allow larger error margin for pressure computes --- unittest/commands/test_compute_global.cpp | 28 +++++++++++------------ 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/unittest/commands/test_compute_global.cpp b/unittest/commands/test_compute_global.cpp index 1c9be99ba4..e21acdbca0 100644 --- a/unittest/commands/test_compute_global.cpp +++ b/unittest/commands/test_compute_global.cpp @@ -102,24 +102,24 @@ TEST_F(ComputeGlobalTest, Energy) EXPECT_DOUBLE_EQ(get_scalar("pe1"), 24155.155261642241); EXPECT_DOUBLE_EQ(get_scalar("pe2"), 361.37528652881286); EXPECT_DOUBLE_EQ(get_scalar("pe3"), 0.0); - EXPECT_DOUBLE_EQ(get_scalar("pr1"), 1956948.4735454607); - EXPECT_DOUBLE_EQ(get_scalar("pr2"), 1956916.7725807722); + EXPECT_NEAR(get_scalar("pr1"), 1956948.4735454607, 0.0000000005); + EXPECT_NEAR(get_scalar("pr2"), 1956916.7725807722, 0.0000000005); EXPECT_DOUBLE_EQ(get_scalar("pr3"), 0.0); auto pr1 = get_vector("pr1"); auto pr2 = get_vector("pr2"); auto pr3 = get_vector("pr3"); - EXPECT_DOUBLE_EQ(pr1[0], 2150600.9207200543); - EXPECT_DOUBLE_EQ(pr1[1], 1466949.7512112649); - EXPECT_DOUBLE_EQ(pr1[2], 2253294.7487050635); - EXPECT_DOUBLE_EQ(pr1[3], 856643.16926486336); - EXPECT_DOUBLE_EQ(pr1[4], 692710.86929464422); - EXPECT_DOUBLE_EQ(pr1[5], -44403.909298603547); - EXPECT_DOUBLE_EQ(pr2[0], 2150575.6989334146); - EXPECT_DOUBLE_EQ(pr2[1], 1466911.3911461537); - EXPECT_DOUBLE_EQ(pr2[2], 2253263.2276627473); - EXPECT_DOUBLE_EQ(pr2[3], 856632.34707690508); - EXPECT_DOUBLE_EQ(pr2[4], 692712.89222328411); - EXPECT_DOUBLE_EQ(pr2[5], -44399.277068014424); + EXPECT_NEAR(pr1[0], 2150600.9207200543, 0.0000000005); + EXPECT_NEAR(pr1[1], 1466949.7512112649, 0.0000000005); + EXPECT_NEAR(pr1[2], 2253294.7487050635, 0.0000000005); + EXPECT_NEAR(pr1[3], 856643.16926486336, 0.0000000005); + EXPECT_NEAR(pr1[4], 692710.86929464422, 0.0000000005); + EXPECT_NEAR(pr1[5], -44403.909298603547, 0.0000000005); + EXPECT_NEAR(pr2[0], 2150575.6989334146, 0.0000000005); + EXPECT_NEAR(pr2[1], 1466911.3911461537, 0.0000000005); + EXPECT_NEAR(pr2[2], 2253263.2276627473, 0.0000000005); + EXPECT_NEAR(pr2[3], 856632.34707690508, 0.0000000005); + EXPECT_NEAR(pr2[4], 692712.89222328411, 0.0000000005); + EXPECT_NEAR(pr2[5], -44399.277068014424, 0.0000000005); EXPECT_DOUBLE_EQ(pr3[0], 0.0); EXPECT_DOUBLE_EQ(pr3[1], 0.0); EXPECT_DOUBLE_EQ(pr3[2], 0.0); From 1568974e8e3216d3dcef85e2c22cd5f4ddc5e828 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 22 Apr 2022 13:39:05 -0400 Subject: [PATCH 5/6] whitespace --- src/angle.cpp | 2 +- src/bond.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/angle.cpp b/src/angle.cpp index caa86d691e..8132146ba4 100644 --- a/src/angle.cpp +++ b/src/angle.cpp @@ -360,7 +360,7 @@ double Angle::memory_usage() void Angle::reinit() { - if (!reinitflag) + if (!reinitflag) error->all(FLERR, "Fix adapt interface to this angle style not supported"); init(); diff --git a/src/bond.cpp b/src/bond.cpp index 7f7b6aa3a3..a208b5a7e4 100644 --- a/src/bond.cpp +++ b/src/bond.cpp @@ -342,7 +342,7 @@ double Bond::memory_usage() void Bond::reinit() { - if (!reinitflag) + if (!reinitflag) error->all(FLERR, "Fix adapt interface to this bond style not supported"); init(); From fec5538d3cbbb1c00d987389a8565c8fb39ffd97 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 22 Apr 2022 13:52:15 -0400 Subject: [PATCH 6/6] fix initialization bugs --- src/fix_adapt.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index 0d9aa00b7a..edaa5c9866 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -44,8 +44,9 @@ enum{DIAMETER,CHARGE}; /* ---------------------------------------------------------------------- */ -FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), -nadapt(0), id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) +FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), nadapt(0), anypair(0), anybond(0), anyangle(0), + id_fix_diam(nullptr), id_fix_chg(nullptr), adapt(nullptr) { if (narg < 5) error->all(FLERR,"Illegal fix adapt command"); nevery = utils::inumeric(FLERR,arg[3],false,lmp); @@ -329,6 +330,7 @@ void FixAdapt::init() anypair = 0; anybond = 0; + anyangle = 0; for (int m = 0; m < nadapt; m++) { Adapt *ad = &adapt[m]; @@ -748,6 +750,7 @@ void FixAdapt::restore_settings() if (anypair) force->pair->reinit(); if (anybond) force->bond->reinit(); + if (anyangle) force->angle->reinit(); if (chgflag && force->kspace) force->kspace->qsum_qsq(); }