tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'develop' into region-lookup-refactor

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-04-22 22:53:02 -04:00
parent 47564ab6bc 9be8ef09a2
commit 9f9a5ca76f
29 changed files with 249 additions and 87 deletions
Download Patch File Download Diff File Expand all files Collapse all files

87
doc/src/fix_adapt.rst
View File

@ -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,62 @@ 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
hybrid/overlay as bond styles. The bond styles that currently work
with fix_adapt are
+------------------------------------+-------+------------+
| :doc:`class2 <bond_class2>` | r0 | type bonds |
+------------------------------------+-------+------------+
| :doc:`fene <bond_fene>` | k, r0 | type bonds |
+------------------------------------+-------+------------+
| :doc:`gromos <bond_gromos>` | k, r0 | type bonds |
+------------------------------------+-------+------------+
| :doc:`harmonic <bond_harmonic>` | k,r0 | type bonds |
+------------------------------------+-------+------------+
| :doc:`morse <bond_morse>` | r0 | type bonds |
+------------------------------------+-------+------------+
| :doc:`nonlinear <bond_nonlinear>` | r0 | type bonds |
+------------------------------------+-------+------------+
+------------------------------------+-------+-----------------+
| :doc:`class2 <bond_class2>` | r0 | type bonds |
+------------------------------------+-------+-----------------+
| :doc:`fene <bond_fene>` | k,r0 | type bonds |
+------------------------------------+-------+-----------------+
| :doc:`fene/nm <bond_fene_nm>` | k,r0 | type bonds |
+------------------------------------+-------+-----------------+
| :doc:`gromos <bond_gromos>` | k,r0 | type bonds |
+------------------------------------+-------+-----------------+
| :doc:`harmonic <bond_harmonic>` | k,r0 | type bonds |
+------------------------------------+-------+-----------------+
| :doc:`morse <bond_morse>` | r0 | type bonds |
+------------------------------------+-------+-----------------+
| :doc:`nonlinear <bond_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 angle styles that currently work
with fix_adapt are
+------------------------------------+-------+-----------------+
| :doc:`harmonic <angle_harmonic>` | k,theta0 | type angles |
+------------------------------------+-------+-----------------+
| :doc:`cosine <angle_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.
----------

2
src/EXTRA-COMPUTE/compute_ave_sphere_atom.h
View File

@ -48,7 +48,7 @@ class ComputeAveSphereAtom : public Compute {
#endif
#endif
/* ERROR/WARNING messages:
/* ERROR/WARNING messages:
E: Illegal ... command

3
src/EXTRA-FIX/fix_oneway.cpp
View File

@ -80,8 +80,7 @@ int FixOneWay::setmask()
void FixOneWay::init()
{
region = domain->get_region_by_id(idregion);
if (!region)
error->all(FLERR, "Region {} for fix oneway does not exist", idregion);
if (!region) error->all(FLERR, "Region {} for fix oneway does not exist", idregion);
}
/* ---------------------------------------------------------------------- */

2
src/EXTRA-MOLECULE/bond_fene_nm.cpp
View File

@ -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;
}

1
src/EXTRA-MOLECULE/bond_gaussian.cpp
View File

@ -35,7 +35,6 @@ BondGaussian::BondGaussian(LAMMPS *lmp) :
Bond(lmp), nterms(nullptr), bond_temperature(nullptr), alpha(nullptr), width(nullptr),
r0(nullptr)
{
reinitflag = 1;
}
/* ---------------------------------------------------------------------- */

11
src/MOLECULE/angle_cosine.cpp
View File

@ -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;
}

1
src/MOLECULE/angle_cosine.h
View File

@ -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;

12
src/MOLECULE/angle_harmonic.cpp
View File

@ -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;
}

1
src/MOLECULE/angle_harmonic.h
View File

@ -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;

2
src/MOLECULE/bond_fene.cpp
View File

@ -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;
}

7
src/MOLECULE/bond_gromos.cpp
View File

@ -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;
}

10
src/MOLECULE/bond_harmonic.cpp
View File

@ -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;
}

3
src/PLUGIN/plugin.cpp
View File

@ -388,8 +388,7 @@ void plugin_unload(const char *style, const char *name, LAMMPS *lmp)
auto found = region_map->find(name);
if (found != region_map->end()) region_map->erase(name);
for (auto iregion : lmp->domain->get_region_by_style(name))
lmp->domain->delete_region(iregion);
for (auto iregion : lmp->domain->get_region_by_style(name)) lmp->domain->delete_region(iregion);
} else if (pstyle == "command") {

2
src/RIGID/fix_ehex.cpp
View File

@ -162,7 +162,7 @@ void FixEHEX::init()
if (idregion) {
region = domain->get_region_by_id(idregion);
if (!region) error->all(FLERR, "Region {} for fix ehex does not exist",idregion);
if (!region) error->all(FLERR, "Region {} for fix ehex does not exist", idregion);
}
// cannot have 0 atoms in group

11
src/angle.cpp
View File

@ -353,3 +353,14 @@ 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();
}

5
src/angle.h
View File

@ -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

4
src/atom_vec_hybrid.h
View File

@ -74,8 +74,8 @@ class AtomVecHybrid : public AtomVec {
int nstyles_bonus;
class AtomVec **styles_bonus;
void merge_fields(std::vector<std::string> &, const std::vector<std::string> &,
int, std::vector<std::string> &);
void merge_fields(std::vector<std::string> &, const std::vector<std::string> &, int,
std::vector<std::string> &);
void build_styles();
int known_style(char *);
};

4
src/bond.cpp
View File

@ -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,8 +337,9 @@ 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");

4
src/bond.h
View File

@ -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
@ -62,6 +63,7 @@ class Bond : protected Pointers {
virtual double memory_usage();
virtual void *extract(const char *, int &) { return nullptr; }
virtual 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; }

109
src/fix_adapt.cpp
View File

@ -14,6 +14,7 @@
#include "fix_adapt.h"
#include "angle.h"
#include "atom.h"
#include "bond.h"
#include "domain.h"
@ -38,13 +39,14 @@ using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
enum{PAIR,KSPACE,ATOM,BOND};
enum{PAIR,KSPACE,ATOM,BOND,ANGLE};
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);
@ -75,6 +77,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 +125,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 +153,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 +211,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 +234,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;
@ -299,6 +330,7 @@ void FixAdapt::init()
anypair = 0;
anybond = 0;
anyangle = 0;
for (int m = 0; m < nadapt; m++) {
Adapt *ad = &adapt[m];
@ -357,6 +389,7 @@ void FixAdapt::init()
}
delete [] pstyle;
} else if (ad->which == BOND) {
ad->bond = nullptr;
anybond = 1;
@ -383,13 +416,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 +457,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 +467,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 +481,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 +546,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 +590,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 +615,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 +642,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 +666,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 +676,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 +700,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];
}
@ -668,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();
}

8
src/fix_adapt.h
View File

@ -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;

4
src/variable.h
View File

@ -118,8 +118,8 @@ class Variable : protected Pointers {
Tree **extra; // ptrs further down tree for nextra args
Tree() :
array(nullptr), iarray(nullptr), barray(nullptr), selfalloc(0), ivalue(0),
nextra(0), region(nullptr), first(nullptr), second(nullptr), extra(nullptr)
array(nullptr), iarray(nullptr), barray(nullptr), selfalloc(0), ivalue(0), nextra(0),
region(nullptr), first(nullptr), second(nullptr), extra(nullptr)
{
}
};

28
unittest/commands/test_compute_global.cpp
View File

@ -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);

3
unittest/force-styles/tests/angle-cosine.yaml
View File

@ -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-

4
unittest/force-styles/tests/angle-harmonic.yaml
View File

@ -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-

2
unittest/force-styles/tests/bond-fene.yaml
View File

@ -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

2
unittest/force-styles/tests/bond-fene_nm.yaml
View File

@ -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

2
unittest/force-styles/tests/bond-gromos.yaml
View File

@ -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

2
unittest/force-styles/tests/bond-harmonic.yaml
View File

@ -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