ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

upgrades to global and local hyper, including more output

Browse Source
This commit is contained in:
Steve Plimpton
2019-06-18 12:29:46 -06:00
parent a39a35af20
commit 08471684f3
10 changed files with 447 additions and 190 deletions
Download Patch File Download Diff File Expand all files Collapse all files

153
doc/src/fix_hyper_local.txt
View File

@ -22,14 +22,19 @@ Dcut = minimum distance between boosted bonds (distance units) :l
alpha = boostostat relaxation time (time units) :l
Btarget = desired time boost factor (unitless) :l
zero or more keyword/value pairs may be appended :l
keyword = {check/ghost} or {check/bias} :l
keyword = {bound} or {reset} or {check/ghost} or {check/bias} :l
{bound} value = Bfrac
Bfrac = -1 or a value >= 0.0
{reset} value = Rfreq
Rfreq = -1 or 0 or timestep value > 0
{check/ghost} values = none
{check/bias} values = Nevery error/warn/ignore :pre
:ule
[Examples:]
fix 1 all hyper/local 1.0 0.3 0.8 300.0 :pre
fix 1 all hyper/local 1.0 0.3 0.8 300.0
fix 1 all hyper/local 1.0 0.3 0.8 300.0 bound 0.1 reset 0 :pre
[Description:]
@ -191,8 +196,8 @@ compensate to reduce the overall prefactor if necessary. However the
Cij are initialized to 1.0 and the boostostatting procedure typically
operates slowly enough that there can be a time period of bad dynamics
if {Vmax} is set too large. A better strategy is to set {Vmax} to the
smallest barrier height for an event (the same as for GHD), so that
the Cij remain near unity.
slightly smaller than the lowest barrier height for an event (the same
as for GHD), so that the Cij remain near unity.
The {Tequil} argument is the temperature at which the system is
simulated; see the comment above about the "fix
@ -274,6 +279,42 @@ rate does not change (as a function of hyper time).
Here is additional information on the optional keywords for this fix.
The {bound} keyword turns on min/max bounds for bias coefficients Cij
for all bonds. Cij is a prefactor for each bond on the bias potential
of maximum strength Vmax. Depending on the choice of {alpha} and
{Btarget} and {Vmax}, the booststatting can cause individual Cij
values to fluctuate. If the fluctuations are too large Cij*Vmax can
exceed low barrier heights and induce bad event dynamics. Bounding
the Cij values is a way to prevent this. If {Bfrac} is set to -1 or
any negative value (the default) then no bounds are enforced on Cij
values (except they must always be >= 0.0). A {Bfrac} setting >= 0.0
sets a lower bound of 1.0 - Bfrac and upper bound of 1.0 + Bfrac on
each Cij value. Note that all Cij values are initialized to 1.0 when
a bond is created for the first time. Thus {Bfrac} limits the bias
potential height to {Vmax} +/- {Bfrac}*{Vmax}.
The {reset} keyword allow {Vmax} to be adjusted dynamically depending
on the average value of all Cij prefactors. This can be useful if you
are unsure what value of {Vmax} will match the {Btarget} boost for the
system. The Cij values will then adjust in aggregate (up or down) so
that Cij*Vmax produces a boost of {Btarget}, but this may conflict
with the {bound} keyword settings. By using {bound} and {reset}
together, {Vmax} itself can be reset, and desired bounds still applied
to the Cij values.
A setting for {Rfreq} of -1 (the default) means {Vmax} never changes.
A setting of 0 means {Vmax} is adjusted every time an event occurs and
bond pairs are recalculated. A setting of N > 0 timesteps means
{Vmax} is adjusted on the first time an event occurs on a timestep >=
N steps after the previous adjustment. The adjustement to {Vmax} is
computed as follows. The current average of all Cij*Vmax values is
computed and the {Vmax} is reset to that value. All Cij values are
changed to new prefactors such the new Cij*Vmax is the same as it was
previously. If the {bound} keyword was used, those bounds are
enforced on the new Cij values. Henceforth, new bonds are assigned a
Cij = 1.0, which means their bias potential magnitude is the new
{Vmax}.
The {check/ghost} keyword turns on extra computation each timestep to
compute statistics about ghost atoms used to determine which bonds to
bias. The output of these stats are the vector values 14 and 15,
@ -301,83 +342,93 @@ No information about this fix is written to "binary restart
files"_restart.html.
The "fix_modify"_fix_modify.html {energy} option is supported by this
fix to add the energy of the bias potential to the system's
potential energy as part of "thermodynamic output"_thermo_style.html.
fix to add the energy of the bias potential to the system's potential
energy as part of "thermodynamic output"_thermo_style.html.
This fix computes a global scalar and global vector of length 21,
This fix computes a global scalar and global vector of length 26,
which can be accessed by various "output commands"_Howto_output.html.
The scalar is the magnitude of the bias potential (energy units)
applied on the current timestep, summed over all biased bonds. The
vector stores the following quantities:
1 = # of biased bonds on this step
2 = max strain Eij of any bond on this step (absolute value, unitless)
3 = average bias coeff for all bonds on this step (unitless)
4 = average # of bonds/atom on this step
5 = average neighbor bonds/bond on this step within {Dcut} :ul
1 = average boost for all bonds on this step (unitless)
2 = # of biased bonds on this step
3 = max strain Eij of any bond on this step (absolute value, unitless)
4 = value of Vmax on this step (energy units)
5 = average bias coeff for all bonds on this step (unitless)
6 = min bias coeff for all bonds on this step (unitless)
7 = max bias coeff for all bonds on this step (unitless)
8 = average # of bonds/atom on this step
9 = average neighbor bonds/bond on this step within {Dcut} :ul
6 = max bond length during this run (distance units)
7 = average # of biased bonds/step during this run
8 = fraction of biased bonds with no bias during this run
9 = fraction of biased bonds with negative strain during this run
10 = average bias coeff for all bonds during this run (unitless)
11 = min bias coeff for any bond during this run (unitless)
12 = max bias coeff for any bond during this run (unitless) :ul
10 = average boost for all bonds during this run (unitless)
11 = average # of biased bonds/step during this run
12 = fraction of biased bonds with no bias during this run
13 = fraction of biased bonds with negative strain during this run
14 = max bond length during this run (distance units)
15 = average bias coeff for all bonds during this run (unitless)
16 = min bias coeff for any bond during this run (unitless)
17 = max bias coeff for any bond during this run (unitless)
13 = max drift distance of any bond atom during this run (distance units)
14 = max distance from proc subbox of any ghost atom with maxstrain < qfactor during this run (distance units)
15 = max distance outside my box of any ghost atom with any maxstrain during this run (distance units)
16 = count of ghost atoms that could not be found on reneighbor steps during this run
17 = count of bias overlaps (< Dcut) found during this run :ul
18 = max drift distance of any bond atom during this run (distance units)
19 = max distance from proc subbox of any ghost atom with maxstrain < qfactor during this run (distance units)
20 = max distance outside my box of any ghost atom with any maxstrain during this run (distance units)
21 = count of ghost atoms that could not be found on reneighbor steps during this run
22 = count of bias overlaps (< Dcut) found during this run
18 = cumulative hyper time since fix created (time units)
19 = cumulative count of event timesteps since fix created
20 = cumulative count of atoms in events since fix created
21 = cumulative # of new bonds formed since fix created :ul
23 = cumulative hyper time since fix created (time units)
24 = cumulative count of event timesteps since fix created
25 = cumulative count of atoms in events since fix created
26 = cumulative # of new bonds formed since fix created :ul
The first quantities (1-5) are for the current timestep. Quantities
6-17 are for the current hyper run. They are reset each time a new
hyper run is performed. Quantities 18-21 are cumulative across
The first quantities 1-9 are for the current timestep. Quantities
10-22 are for the current hyper run. They are reset each time a new
hyper run is performed. Quantities 23-26 are cumulative across
multiple runs (since the point in the input script the fix was
defined).
For value 8, the numerator is a count of all biased bonds on each
For value 10, each bond instantaneous boost factor is given by the
equation for Bij above. The total system boost (average across all
bonds) fluctuates, but should average to a value close to the
speficied Btarget.
For value 12, the numerator is a count of all biased bonds on each
timestep whose bias energy = 0.0 due to Eij >= {qfactor}. The
denominator is the count of all biased bonds on all timesteps.
For value 9, the numerator is a count of all biased bonds on each
For value 13, the numerator is a count of all biased bonds on each
timestep with negative strain. The denominator is the count of all
biased bonds on all timesteps.
Values 13-17 are mostly useful for debugging and diagnostic purposes.
Values 18-22 are mostly useful for debugging and diagnostic purposes.
For value 13, drift is the distance an atom moves between two quenched
For value 18, drift is the distance an atom moves between two quenched
states when the second quench determines an event has occurred. Atoms
involved in an event will typically move the greatest distance since
others typically remain near their original quenched position.
For values 14-16, neighbor atoms in the full neighbor list with cutoff
For values 19-21, neighbor atoms in the full neighbor list with cutoff
{Dcut} may be ghost atoms outside a processor's sub-box. Before the
next event occurs they may move further than {Dcut} away from the
sub-box boundary. Value 14 is the furthest (from the sub-box) any
sub-box boundary. Value 19 is the furthest (from the sub-box) any
ghost atom in the neighbor list with maxstrain < {qfactor} was
accessed during the run. Value 15 is the same except that the ghost
accessed during the run. Value 20 is the same except that the ghost
atom's maxstrain may be >= {qfactor}, which may mean it is about to
participate in an event. Value 16 is a count of how many ghost atoms
participate in an event. Value 21 is a count of how many ghost atoms
could not be found on reneighbor steps, presumably because they moved
too far away due to their participation in an event (which will likely
be detected at the next quench).
Typical values for 14 and 15 should be slightly larger than {Dcut},
Typical values for 19 and 20 should be slightly larger than {Dcut},
which accounts for ghost atoms initially at a {Dcut} distance moving
thermally before the next event takes place.
Note that for values 14 and 15 to be computed, the optional keyword
Note that for values 19 and 20 to be computed, the optional keyword
{check/ghost} must be specified. Otherwise these values will be zero.
This is because computing them incurs overhead, so the values are only
computed if requested.
Value 16 should be zero or small. As explained above a small count
Value 21 should be zero or small. As explained above a small count
likely means some ghost atoms were participating in their own events
and moved a longer distance. If the value is large, it likely means
the communication cutoff for ghosts is too close to {Dcut} leading to
@ -387,11 +438,11 @@ assumes those atoms are part of highly strained bonds. As explained
above, the "comm_modify cutoff"_comm_modify.html command can be used
to set a longer cutoff.
For value 17, no two bonds should be biased if they are within a
For value 22, no two bonds should be biased if they are within a
{Dcut} distance of each other. This value should be zero, indicating
that no pair of biased bonds are closer than {Dcut} from each other.
Note that for values 17 to be computed, the optional keyword
Note that for value 22 to be computed, the optional keyword
{check/bias} must be specified and it determines how often this check
is performed. This is because performing the check incurs overhead,
so if only computed as often as requested.
@ -401,23 +452,23 @@ timestep the check was made. Note that the value is a count of atoms
in bonds which found other atoms in bonds too close, so it is almost
always an over-count of the number of too-close bonds.
Value 18 is simply the specified {boost} factor times the number of
Value 23 is simply the specified {boost} factor times the number of
timesteps times the timestep size.
For value 19, events are checked for by the "hyper"_hyper.html command
For value 24, events are checked for by the "hyper"_hyper.html command
once every {Nevent} timesteps. This value is the count of those
timesteps on which one (or more) events was detected. It is NOT the
number of distinct events, since more than one event may occur in the
same {Nevent} time window.
For value 20, each time the "hyper"_hyper.html command checks for an
For value 25, each time the "hyper"_hyper.html command checks for an
event, it invokes a compute to flag zero or more atoms as
participating in one or more events. E.g. atoms that have displaced
more than some distance from the previous quench state. Value 20 is
more than some distance from the previous quench state. Value 25 is
the cumulative count of the number of atoms participating in any of
the events that were found.
Value 21 tallies the number of new bonds created by the bond reset
Value 26 tallies the number of new bonds created by the bond reset
operation. Bonds between a specific I,J pair of atoms may persist for
the entire hyperdynamics simulation if neither I or J are involved in
an event.
@ -451,7 +502,9 @@ doc page for more info.
[Default:]
The check/ghost and check/bias keywords are not enabled by default.
The default settings for optinal keywords are bounds = -1 and reset =
-1. The check/ghost and check/bias keywords are not enabled by
default.
:line

24
src/REPLICA/fix_event_hyper.cpp
View File

@ -63,12 +63,12 @@ void FixEventHyper::write_restart(FILE *fp)
{
int n = 0;
double list[6];
list[n++] = event_number;
list[n++] = event_timestep;
list[n++] = clock;
list[n++] = replica_number;
list[n++] = correlated_event;
list[n++] = ncoincident;
list[n++] = ubuf(event_number).d;
list[n++] = ubuf(event_timestep).d;
list[n++] = ubuf(clock).d;
list[n++] = ubuf(replica_number).d;
list[n++] = ubuf(correlated_event).d;
list[n++] = ubuf(ncoincident).d;
if (comm->me == 0) {
int size = n * sizeof(double);
@ -86,10 +86,10 @@ void FixEventHyper::restart(char *buf)
int n = 0;
double *list = (double *) buf;
event_number = static_cast<int> (list[n++]);
event_timestep = static_cast<bigint> (list[n++]);
clock = static_cast<bigint> (list[n++]);
replica_number = static_cast<int> (list[n++]);
correlated_event = static_cast<int> (list[n++]);
ncoincident = static_cast<int> (list[n++]);
event_number = (int) ubuf(list[n++]).i;
event_timestep = (bigint) ubuf(list[n++]).i;
clock = (bigint) ubuf(list[n++]).i;
replica_number = (int) ubuf(list[n++]).i;
correlated_event = (int) ubuf(list[n++]).i;
ncoincident = (int) ubuf(list[n++]).i;
}

8
src/REPLICA/fix_hyper.cpp
View File

@ -13,13 +13,17 @@
#include <cstring>
#include "fix_hyper.h"
#include "update.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixHyper::FixHyper(LAMMPS *lmp, int narg, char **arg)
: Fix(lmp, narg, arg), hyperflag(0) {}
FixHyper::FixHyper(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
hyperflag = 0;
ntimestep_initial = update->ntimestep;
}
/* ----------------------------------------------------------------------
extract hyper flag setting for all Fixes that perform hyperdynamics

2
src/REPLICA/fix_hyper.h
View File

@ -20,6 +20,8 @@ namespace LAMMPS_NS {
class FixHyper : public Fix {
public:
bigint ntimestep_initial;
FixHyper(class LAMMPS *, int, char **);
virtual ~FixHyper() {}
void *extract(const char *, int &);

7
src/REPLICA/fix_hyper_global.cpp
View File

@ -18,6 +18,7 @@
#include "fix_hyper_global.h"
#include "atom.h"
#include "update.h"
#include "group.h"
#include "force.h"
#include "domain.h"
#include "comm.h"
@ -132,6 +133,10 @@ void FixHyperGlobal::init()
dt = update->dt;
// count of atoms in fix group
groupatoms = group->count(igroup);
// need an occasional half neighbor list
int irequest = neighbor->request(this,instance_me);
@ -493,7 +498,7 @@ double FixHyperGlobal::compute_vector(int i)
bigint mybonds = nblocal;
bigint allbonds;
MPI_Allreduce(&mybonds,&allbonds,1,MPI_LMP_BIGINT,MPI_SUM,world);
return 2.0*allbonds/atom->natoms;
return 1.0*allbonds/groupatoms;
}
if (i == 5) {

1
src/REPLICA/fix_hyper_global.h
View File

@ -57,6 +57,7 @@ class FixHyperGlobal : public FixHyper {
double maxdriftsq; // max distance any atom drifts from original pos
int nobias; // # of steps when bias = 0, b/c bond too long
int negstrain; // # of steps when biased bond has negative strain
bigint groupatoms; // # of atoms in fix group
class NeighList *list;

371
src/REPLICA/fix_hyper_local.cpp
View File

@ -18,6 +18,7 @@
#include "fix_hyper_local.h"
#include "atom.h"
#include "update.h"
#include "group.h"
#include "force.h"
#include "pair.h"
#include "domain.h"
@ -62,7 +63,8 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) :
hyperflag = 2;
scalar_flag = 1;
vector_flag = 1;
size_vector = 21;
// DEBUG - changed 26 to 28
size_vector = 28;
local_flag = 1;
size_local_rows = 0;
size_local_cols = 0;
@ -84,6 +86,7 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) :
tequil <= 0.0 || dcut <= 0.0 || alpha_user <= 0.0 || boost_target < 1.0)
error->all(FLERR,"Illegal fix hyper/local command");
invvmax = 1.0 / vmax;
invqfactorsq = 1.0 / (qfactor*qfactor);
cutbondsq = cutbond*cutbond;
dcutsq = dcut*dcut;
@ -91,12 +94,26 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) :
// optional args
boundflag = 0;
boundfrac = 0.0;
resetfreq = -1;
checkghost = 0;
checkbias = 0;
int iarg = 10;
while (iarg < narg) {
if (strcmp(arg[iarg],"check/ghost") == 0) {
if (strcmp(arg[iarg],"bound") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix hyper/local command");
boundfrac = force->numeric(FLERR,arg[iarg+1]);
if (boundfrac < 0.0) boundflag = 0;
else boundflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"reset") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix hyper/local command");
resetfreq = force->inumeric(FLERR,arg[iarg+1]);
if (resetfreq < -1) error->all(FLERR,"Illegal fix hyper/local command");
iarg += 2;
} else if (strcmp(arg[iarg],"check/ghost") == 0) {
checkghost = 1;
iarg++;
} else if (strcmp(arg[iarg],"check/bias") == 0) {
@ -159,8 +176,12 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) :
me = comm->me;
firstflag = 1;
sumboost = 0.0;
aveboost_running = 0.0;
sumbiascoeff = 0.0;
avebiascoeff = 0.0;
avebiascoeff_running = 0.0;
minbiascoeff = 0.0;
maxbiascoeff = 0.0;
starttime = update->ntimestep;
nostrainyet = 1;
@ -168,6 +189,15 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) :
nevent = 0;
nevent_atom = 0;
mybias = 0.0;
// bias bounding and reset params
bound_lower = 1.0 - boundfrac;
bound_upper = 1.0 + boundfrac;
lastreset = update->ntimestep;
// DEBUG - one line
overcount = 0;
}
/* ---------------------------------------------------------------------- */
@ -217,9 +247,9 @@ void FixHyperLocal::init_hyper()
checkbias_count = 0;
maxdriftsq = 0.0;
maxbondlen = 0.0;
avebiascoeff = 0.0;
minbiascoeff = BIG;
maxbiascoeff = 0.0;
avebiascoeff_running = 0.0;
minbiascoeff_running = BIG;
maxbiascoeff_running = 0.0;
nbias_running = 0;
nobias_running = 0;
negstrain_running = 0;
@ -267,6 +297,10 @@ void FixHyperLocal::init()
alpha = update->dt / alpha_user;
// count of atoms in fix group
groupatoms = group->count(igroup);
// need occasional full neighbor list with cutoff = Dcut
// used for finding maxstrain of neighbor bonds out to Dcut
// do not need to include neigh skin in cutoff,
@ -451,6 +485,9 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */)
double **x = atom->x;
// DEBUG - one line
overcount = 0;
m = 0;
for (iold = 0; iold < nlocal_old; iold++) {
nbond = numbond[iold];
@ -469,6 +506,8 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */)
maxbondlen = MAX(r,maxbondlen);
r0 = blist[m].r0;
estrain = fabs(r-r0) / r0;
// DEBUG - one line
if (estrain >= qfactor) overcount++;
maxstrain[i] = MAX(maxstrain[i],estrain);
maxstrain[j] = MAX(maxstrain[j],estrain);
if (estrain > halfstrain) {
@ -613,7 +652,7 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */)
// maxstrain[i] = maxstrain_domain[i] (checked in stage 2)
// maxstrain[j] = maxstrain_domain[j] (checked here)
// I is not part of an I,J bond with > strain owned by some J (checked in 2)
// no ties with other maxstrain bonds in atom I's domain (chedcked in 2)
// no ties with other maxstrain bonds in atom I's domain (checked in 2)
nbias = 0;
for (iold = 0; iold < nlocal_old; iold++) {
@ -640,12 +679,17 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */)
int negstrain = 0;
mybias = 0.0;
// DEBUG - one line
myboost = 0;
for (int ibias = 0; ibias < nbias; ibias++) {
m = bias[ibias];
i = blist[m].i;
j = blist[m].j;
if (maxstrain[i] >= qfactor) {
// DEBUG - one line
myboost += 1.0;
nobias++;
continue;
}
@ -657,8 +701,7 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */)
r0 = blist[m].r0;
ebias = (r-r0) / r0;
vbias = biascoeff[m] * vmax * (1.0 - ebias*ebias*invqfactorsq);
fbias = biascoeff[m] * 2.0 * vmax * ebias * invqfactorsq;
fbias = 2.0 * biascoeff[m] * vmax * ebias * invqfactorsq;
fbiasr = fbias / r0 / r;
f[i][0] += delx*fbiasr;
f[i][1] += dely*fbiasr;
@ -670,11 +713,15 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */)
if (ebias < 0.0) negstrain++;
mybias += vbias;
// DEBUG - one line
myboost += exp(beta * biascoeff[m]*vbias);
}
//time7 = MPI_Wtime();
// -------------------------------------------------------------
// stage 5:
// apply boostostat to bias coeffs of all bonds I own
// -------------------------------------------------------------
@ -683,16 +730,22 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */)
// on run steps
if (setupflag) return;
nbias_running += nbias;
nobias_running += nobias;
negstrain_running += negstrain;
// loop over bonds I own to adjust bias coeff
// delta in boost coeff is function of maxboost_domain vs target boost
// maxboost_domain is function of two maxstrain_domains for I,J
// delta in boost coeff is function of boost_domain vs target boost
// boost_domain is function of two maxstrain_domains for I,J
// NOTE: biascoeff update is now scaled by 1/Vmax
double emaxi,emaxj,maxboost_domain,bc;
double mybiascoeff = 0.0;
minbiascoeff = BIG;
maxbiascoeff = 0.0;
double emaxi,emaxj,boost_domain,bc;
double sumcoeff_me = 0.0;
double sumboost_me = 0.0;
for (m = 0; m < nblocal; m++) {
i = blist[m].i;
@ -703,23 +756,37 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */)
if (emax < qfactor) vbias = vmax * (1.0 - emax*emax*invqfactorsq);
else vbias = 0.0;
maxboost_domain = exp(beta * biascoeff[m]*vbias);
biascoeff[m] -= alpha * (maxboost_domain-boost_target) / boost_target;
boost_domain = exp(beta * biascoeff[m]*vbias);
biascoeff[m] -= alpha*invvmax * (boost_domain-boost_target) / boost_target;
// enforce biascoeff bounds
// min value must always be >= 0.0
biascoeff[m] = MAX(biascoeff[m],0.0);
if (boundflag) {
biascoeff[m] = MAX(biascoeff[m],bound_lower);
biascoeff[m] = MIN(biascoeff[m],bound_upper);
}
// stats
bc = biascoeff[m];
mybiascoeff += bc;
sumcoeff_me += bc;
minbiascoeff = MIN(minbiascoeff,bc);
maxbiascoeff = MAX(maxbiascoeff,bc);
sumboost_me += boost_domain;
}
// -------------------------------------------------------------
// diagnostics, some optional
// -------------------------------------------------------------
MPI_Allreduce(&mybiascoeff,&sumbiascoeff,1,MPI_DOUBLE,MPI_SUM,world);
if (allbonds) avebiascoeff += sumbiascoeff/allbonds;
MPI_Allreduce(&sumcoeff_me,&sumbiascoeff,1,MPI_DOUBLE,MPI_SUM,world);
if (allbonds) avebiascoeff_running += sumbiascoeff/allbonds;
minbiascoeff_running = MIN(minbiascoeff_running,minbiascoeff);
maxbiascoeff_running = MAX(maxbiascoeff_running,maxbiascoeff);
MPI_Allreduce(&sumboost_me,&sumboost,1,MPI_DOUBLE,MPI_SUM,world);
if (allbonds) aveboost_running += sumboost/allbonds;
// if requested, monitor ghost distance from processor sub-boxes
@ -821,6 +888,42 @@ void FixHyperLocal::build_bond_list(int natom)
nevent_atom += natom;
}
// reset Vmax to current bias coeff average
// only if requested and elapsed time >= resetfreq
// ave = curent ave of all bias coeffs
// if reset, adjust all Cij to keep Cij*Vmax unchanged
if (resetfreq >= 0) {
int flag = 0;
bigint elapsed = update->ntimestep - lastreset;
if (resetfreq == 0) {
if (elapsed) flag = 1;
} else {
if (elapsed >= resetfreq) flag = 1;
}
if (flag) {
lastreset = update->ntimestep;
double vmaxold = vmax;
double ave = sumbiascoeff / allbonds;
vmax *= ave;
// adjust all Cij to keep Cij * Vmax = Cijold * Vmaxold
for (m = 0; m < nblocal; m++) biascoeff[m] *= vmaxold/vmax;
// enforce bounds for new Cij
if (boundflag) {
for (m = 0; m < nblocal; m++) {
biascoeff[m] = MAX(biascoeff[m],bound_lower);
biascoeff[m] = MIN(biascoeff[m],bound_upper);
}
}
}
}
// compute max distance any bond atom has moved between 2 quenched states
// xold[iold] = last quenched coord for iold
// x[ilocal] = current quenched coord for same atom
@ -1073,6 +1176,9 @@ void FixHyperLocal::build_bond_list(int natom)
bigint bondcount = nblocal;
MPI_Allreduce(&bondcount,&allbonds,1,MPI_LMP_BIGINT,MPI_SUM,world);
// DEBUG
//if (me == 0) printf("TOTAL BOND COUNT = %ld\n",allbonds);
time2 = MPI_Wtime();
if (firstflag) nnewbond = 0;
@ -1360,43 +1466,57 @@ double FixHyperLocal::compute_scalar()
double FixHyperLocal::compute_vector(int i)
{
// 21 vector outputs returned for i = 0-20
// 26 vector outputs returned for i = 0-25
// DEBUG = 28
// i = 0 = # of biased bonds on this step
// i = 1 = max strain of any bond on this step
// i = 2 = average bias coeff for all bonds on this step
// i = 3 = ave bonds/atom on this step
// i = 4 = ave neighbor bonds/bond on this step
// i = 0 = average boost for all bonds on this step
// i = 1 = # of biased bonds on this step
// i = 2 = max strain of any bond on this step
// i = 3 = value of Vmax on this step
// i = 4 = average bias coeff for all bonds on this step
// i = 5 = min bias coeff for any bond on this step
// i = 6 = max bias coeff for any bond on this step
// i = 7 = ave bonds/atom on this step
// i = 8 = ave neighbor bonds/bond on this step
// i = 5 = max bond length during this run
// i = 6 = average # of biased bonds/step during this run
// i = 7 = fraction of biased bonds with no bias during this run
// i = 8 = fraction of biased bonds with negative strain during this run
// i = 9 = average bias coeff for all bonds during this run
// i = 10 = min bias coeff for any bond during this run
// i = 11 = max bias coeff for any bond during this run
// i = 9 = average boost for all bonds during this run
// i = 10 = average # of biased bonds/step during this run
// i = 11 = fraction of biased bonds with no bias during this run
// i = 12 = fraction of biased bonds with negative strain during this run
// i = 13 = max bond length during this run
// i = 14 = average bias coeff for all bonds during this run
// i = 15 = min bias coeff for any bond during this run
// i = 16 = max bias coeff for any bond during this run
// i = 12 = max drift distance of any atom during this run
// i = 13 = max distance from proc subbox of any ghost atom with
// i = 17 = max drift distance of any atom during this run
// i = 18 = max distance from proc subbox of any ghost atom with
// maxstrain < qfactor during this run
// i = 14 = max distance from proc subbox of any ghost atom with
// i = 19 = max distance from proc subbox of any ghost atom with
// any maxstrain during this run
// i = 15 = count of ghost atoms that could not be found
// i = 20 = count of ghost atoms that could not be found
// on reneighbor steps during this run
// i = 16 = count of bias overlaps (< Dcut) found during this run
// i = 21 = count of bias overlaps (< Dcut) found during this run
// i = 17 = cumulative hyper time since fix created
// i = 18 = cumulative # of event timesteps since fix created
// i = 19 = cumulative # of atoms in events since fix created
// i = 20 = cumulative # of new bonds formed since fix created
// i = 22 = cumulative hyper time since fix created
// i = 23 = cumulative # of event timesteps since fix created
// i = 24 = cumulative # of atoms in events since fix created
// i = 25 = cumulative # of new bonds formed since fix created
// i = 26 = average boost for biased bonds on this step (unitless)
// i = 27 = current count of bonds with strain >= q
if (i == 0) {
if (allbonds) return sumboost/allbonds;
return 1.0;
}
if (i == 1) {
int nbiasall;
MPI_Allreduce(&nbias,&nbiasall,1,MPI_INT,MPI_SUM,world);
return (double) nbiasall;
}
if (i == 1) {
if (i == 2) {
if (nostrainyet) return 0.0;
int nlocal = atom->nlocal;
double emax = 0.0;
@ -1407,39 +1527,58 @@ double FixHyperLocal::compute_vector(int i)
return eall;
}
if (i == 2) {
if (i == 3) {
return vmax;
}
if (i == 4) {
if (allbonds) return sumbiascoeff/allbonds;
return 1.0;
}
if (i == 3) {
return 2.0*allbonds/atom->natoms;
}
if (i == 4) {
bigint allneigh,thisneigh;
thisneigh = listfull->ipage->ndatum;
MPI_Allreduce(&thisneigh,&allneigh,1,MPI_LMP_BIGINT,MPI_SUM,world);
const double natoms = atom->natoms;
const double neighsperatom = static_cast<double>(allneigh)/natoms;
const double bondsperatom = static_cast<double>(allbonds)/natoms;
return neighsperatom * bondsperatom;
}
if (i == 5) {
double allbondlen;
MPI_Allreduce(&maxbondlen,&allbondlen,1,MPI_DOUBLE,MPI_MAX,world);
return allbondlen;
double coeff;
MPI_Allreduce(&minbiascoeff,&coeff,1,MPI_DOUBLE,MPI_MIN,world);
return coeff;
}
if (i == 6) {
double coeff;
MPI_Allreduce(&maxbiascoeff,&coeff,1,MPI_DOUBLE,MPI_MAX,world);
return coeff;
}
if (i == 7) return 1.0*allbonds/groupatoms;
if (i == 8) {
bigint allneigh,thisneigh;
thisneigh = listfull->ipage->ndatum;
MPI_Allreduce(&thisneigh,&allneigh,1,MPI_LMP_BIGINT,MPI_SUM,world);
double natoms = atom->natoms;
double neighsperatom = 1.0*allneigh/natoms;
double bondsperatom = 1.0*allbonds/groupatoms;
return neighsperatom * bondsperatom;
}
// during minimization, just output previous value
if (i == 9) {
if (update->ntimestep == update->firststep)
aveboost_running_output = 0.0;
else if (update->whichflag == 1)
aveboost_running_output =
aveboost_running / (update->ntimestep - update->firststep);
return aveboost_running_output;
}
if (i == 10) {
if (update->ntimestep == update->firststep) return 0.0;
int allbias_running;
MPI_Allreduce(&nbias_running,&allbias_running,1,MPI_INT,MPI_SUM,world);
return 1.0*allbias_running / (update->ntimestep - update->firststep);
}
if (i == 7) {
if (i == 11) {
int allbias_running,allnobias_running;
MPI_Allreduce(&nbias_running,&allbias_running,1,MPI_INT,MPI_SUM,world);
MPI_Allreduce(&nobias_running,&allnobias_running,1,MPI_INT,MPI_SUM,world);
@ -1447,7 +1586,7 @@ double FixHyperLocal::compute_vector(int i)
return 0.0;
}
if (i == 8) {
if (i == 12) {
int allbias_running,allnegstrain_running;
MPI_Allreduce(&nbias_running,&allbias_running,1,MPI_INT,MPI_SUM,world);
MPI_Allreduce(&negstrain_running,&allnegstrain_running,1,MPI_INT,
@ -1456,57 +1595,86 @@ double FixHyperLocal::compute_vector(int i)
return 0.0;
}
if (i == 9) {
if (update->ntimestep == update->firststep) return 0.0;
return avebiascoeff / (update->ntimestep - update->firststep);
if (i == 13) {
double allbondlen;
MPI_Allreduce(&maxbondlen,&allbondlen,1,MPI_DOUBLE,MPI_MAX,world);
return allbondlen;
}
if (i == 10) {
double biascoeff;
MPI_Allreduce(&minbiascoeff,&biascoeff,1,MPI_DOUBLE,MPI_MIN,world);
return biascoeff;
// during minimization, just output previous value
if (i == 14) {
if (update->ntimestep == update->firststep)
avebiascoeff_running_output = 0.0;
else if (update->whichflag == 1)
avebiascoeff_running_output =
avebiascoeff_running / (update->ntimestep - update->firststep);
return avebiascoeff_running_output;
}
if (i == 11) {
double biascoeff;
MPI_Allreduce(&maxbiascoeff,&biascoeff,1,MPI_DOUBLE,MPI_MAX,world);
return biascoeff;
if (i == 15) {
double coeff;
MPI_Allreduce(&minbiascoeff_running,&coeff,1,MPI_DOUBLE,MPI_MIN,world);
return coeff;
}
if (i == 12) {
if (i == 16) {
double coeff;
MPI_Allreduce(&maxbiascoeff_running,&coeff,1,MPI_DOUBLE,MPI_MAX,world);
return coeff;
}
if (i == 17) {
double alldriftsq;
MPI_Allreduce(&maxdriftsq,&alldriftsq,1,MPI_DOUBLE,MPI_MAX,world);
return (double) sqrt(alldriftsq);
}
if (i == 13) return rmaxever;
if (i == 14) return rmaxeverbig;
if (i == 18) return rmaxever;
if (i == 19) return rmaxeverbig;
if (i == 15) {
if (i == 20) {
int allghost_toofar;
MPI_Allreduce(&ghost_toofar,&allghost_toofar,1,MPI_INT,MPI_SUM,world);
return 1.0*allghost_toofar;
}
if (i == 16) {
if (i == 21) {
int allclose;
MPI_Allreduce(&checkbias_count,&allclose,1,MPI_INT,MPI_SUM,world);
return 1.0*allclose;
}
if (i == 17) {
if (i == 22) {
return boost_target * update->dt * (update->ntimestep - starttime);
}
if (i == 18) return (double) nevent;
if (i == 19) return (double) nevent_atom;
if (i == 23) return (double) nevent;
if (i == 24) return (double) nevent_atom;
if (i == 20) {
if (i == 25) {
int allnewbond;
MPI_Allreduce(&nnewbond,&allnewbond,1,MPI_INT,MPI_SUM,world);
return (double) allnewbond;
}
// DEBUG - added two options
if (i == 26) {
double allboost;
MPI_Allreduce(&myboost,&allboost,1,MPI_DOUBLE,MPI_SUM,world);
int nbiasall;
MPI_Allreduce(&nbias,&nbiasall,1,MPI_INT,MPI_SUM,world);
if (nbiasall) return (double) allboost/nbiasall;
return 1.0;
}
if (i == 27) {
int allovercount;
MPI_Allreduce(&overcount,&allovercount,1,MPI_INT,MPI_SUM,world);
return (double) allovercount;
}
return 0.0;
}
@ -1517,30 +1685,31 @@ double FixHyperLocal::compute_vector(int i)
double FixHyperLocal::query(int i)
{
if (i == 1) return compute_vector(17); // cummulative hyper time
if (i == 2) return compute_vector(18); // nevent
if (i == 3) return compute_vector(19); // nevent_atom
if (i == 4) return compute_vector(3); // ave bonds/atom
if (i == 5) return compute_vector(12); // maxdrift
if (i == 6) return compute_vector(5); // maxbondlen
if (i == 7) return compute_vector(7); // fraction with zero bias
if (i == 8) return compute_vector(8); // fraction with negative strain
if (i == 1) return compute_vector(22); // cummulative hyper time
if (i == 2) return compute_vector(23); // nevent
if (i == 3) return compute_vector(24); // nevent_atom
if (i == 4) return compute_vector(7); // ave bonds/atom
if (i == 5) return compute_vector(17); // maxdrift
if (i == 6) return compute_vector(13); // maxbondlen
if (i == 7) return compute_vector(11); // fraction with zero bias
if (i == 8) return compute_vector(12); // fraction with negative strain
// unique to local hyper
if (i == 9) return compute_vector(20); // number of new bonds
if (i == 9) return compute_vector(25); // number of new bonds
if (i == 10) return 1.0*maxbondperatom; // max bonds/atom
if (i == 11) return compute_vector(6); // ave # of biased bonds/step
if (i == 12) return compute_vector(9); // ave bias coeff over all bonds
if (i == 13) return compute_vector(10); // min bias cooef for any bond
if (i == 14) return compute_vector(11); // max bias cooef for any bond
if (i == 15) return compute_vector(4); // neighbor bonds/bond
if (i == 16) return compute_vector(2); // ave bias coeff now
if (i == 17) return time_bondbuild; // CPU time for build_bond calls
if (i == 18) return rmaxever; // ghost atom distance for < maxstrain
if (i == 19) return rmaxeverbig; // ghost atom distance for any strain
if (i == 20) return compute_vector(15); // count of ghost atoms not found
if (i == 21) return compute_vector(16); // count of bias overlaps
if (i == 11) return compute_vector(9); // ave boost/step
if (i == 12) return compute_vector(10); // ave # of biased bonds/step
if (i == 13) return compute_vector(14); // ave bias coeff over all bonds
if (i == 14) return compute_vector(15); // min bias cooef for any bond
if (i == 15) return compute_vector(16); // max bias cooef for any bond
if (i == 16) return compute_vector(8); // neighbor bonds/bond
if (i == 17) return compute_vector(4); // ave bias coeff now
if (i == 18) return time_bondbuild; // CPU time for build_bond calls
if (i == 19) return rmaxever; // ghost atom distance for < maxstrain
if (i == 20) return rmaxeverbig; // ghost atom distance for any strain
if (i == 21) return compute_vector(20); // count of ghost atoms not found
if (i == 22) return compute_vector(21); // count of bias overlaps
error->all(FLERR,"Invalid query to fix hyper/local");

31
src/REPLICA/fix_hyper_local.h
View File

@ -64,9 +64,18 @@ class FixHyperLocal : public FixHyper {
double alpha; // unitless dt/alpha_user
double boost_target; // target value of boost
int checkghost,checkbias; // flags for optional stats
int checkbias_every;
int checkbias_flag;
int boundflag,resetfreq; // bias coeff bounding and reset settings
double boundfrac;
bigint groupatoms; // # of atoms in fix group
double cutbondsq,dcutsq;
double beta,invqfactorsq;
double beta,invvmax,invqfactorsq;
// DEBUG - 2 lines
int overcount;
double myboost;
// flags
@ -76,6 +85,11 @@ class FixHyperLocal : public FixHyper {
bigint starttime; // timestep when this fix was invoked
int commflag; // flag for communication mode
// bias coeff bounds and reset
double bound_lower,bound_upper;
bigint lastreset;
// stats
int nbondbuild; // # of rebuilds of bond list
@ -95,15 +109,20 @@ class FixHyperLocal : public FixHyper {
double maxbondlen; // cummulative max length of any bond
double maxdriftsq; // max distance any bond atom drifts from quenched x
double sumboost; // sum of all bond boosts at each timestep
double aveboost_running; // cummulative sumboost/allbonds across steps
double aveboost_running_output; // most recent output of ab_running
double sumbiascoeff; // sum of all bond bias coeffs at each timestep
double avebiascoeff; // cummulative sumbiascoeff/allbonds across steps
double minbiascoeff; // cummulative min bias coeff for any bond
double maxbiascoeff; // cummulative max bias coeff for any bond
double avebiascoeff_running; // cummulative sumbiascoeff/allbonds across steps
double avebiascoeff_running_output; // most recent output of abc_running
double minbiascoeff; // min bias coeff on this step for my bonds
double maxbiascoeff; // max bias coeff on this step for my bonds
double minbiascoeff_running; // cummulative min bias coeff for any bond
double maxbiascoeff_running; // cummulative max bias coeff for any bond
double rmaxever,rmaxeverbig;
int ghost_toofar; // # of ghost atoms not found in Dcut neigh list
int checkbias_every,checkbias_flag,checkbias_count;
int checkbias_count; // count of too-close biased bonds
// 2 neighbor lists

34
src/REPLICA/hyper.cpp
View File

@ -142,7 +142,7 @@ void Hyper::command(int narg, char **arg)
update->whichflag = 1;
update->nsteps = nsteps;
update->beginstep = update->firststep = update->ntimestep;
update->endstep = update->laststep = update->firststep + nsteps;
update->endstep = update->laststep = update->beginstep + nsteps;
if (update->laststep < 0)
error->all(FLERR,"Too many timesteps");
@ -265,7 +265,7 @@ void Hyper::command(int narg, char **arg)
double fraczero = 1.0;
double fracneg = 1.0;
double nnewbond,avenbias,avebiascoeff,minbiascoeff,maxbiascoeff;
double nnewbond,aveboost,avenbias,avebiascoeff,minbiascoeff,maxbiascoeff;
double maxbondperatom,neighbondperbond,avebiasnow;
double tbondbuild,rmaxever,rmaxeverbig,allghost_toofar;
double biasoverlap;
@ -283,17 +283,18 @@ void Hyper::command(int narg, char **arg)
if (hyperstyle == LOCAL) {
nnewbond = fix_hyper->query(9);
maxbondperatom = fix_hyper->query(10);
avenbias = fix_hyper->query(11);
avebiascoeff = fix_hyper->query(12);
minbiascoeff = fix_hyper->query(13);
maxbiascoeff = fix_hyper->query(14);
neighbondperbond = fix_hyper->query(15);
avebiasnow = fix_hyper->query(16);
tbondbuild = fix_hyper->query(17);
rmaxever = fix_hyper->query(18);
rmaxeverbig = fix_hyper->query(19);
allghost_toofar = fix_hyper->query(20);
biasoverlap = fix_hyper->query(21);
aveboost = fix_hyper->query(11);
avenbias = fix_hyper->query(12);
avebiascoeff = fix_hyper->query(13);
minbiascoeff = fix_hyper->query(14);
maxbiascoeff = fix_hyper->query(15);
neighbondperbond = fix_hyper->query(16);
avebiasnow = fix_hyper->query(17);
tbondbuild = fix_hyper->query(18);
rmaxever = fix_hyper->query(19);
rmaxeverbig = fix_hyper->query(20);
allghost_toofar = fix_hyper->query(21);
biasoverlap = fix_hyper->query(22);
}
}
@ -305,7 +306,8 @@ void Hyper::command(int narg, char **arg)
if (!out) continue;
fprintf(out,"Cummulative quantities for fix hyper:\n");
fprintf(out," hyper time = %g\n",t_hyper);
fprintf(out," time boost factor = %g\n",t_hyper/(nsteps*update->dt));
fprintf(out," time boost factor = %g\n", t_hyper /
((update->ntimestep-fix_hyper->ntimestep_initial)*update->dt));
fprintf(out," event timesteps = %d\n",nevent_running);
fprintf(out," # of atoms in events = %d\n",nevent_atoms_running);
fprintf(out,"Quantities for this hyper run:\n");
@ -325,6 +327,7 @@ void Hyper::command(int narg, char **arg)
fprintf(out," max bonds/atom = %g\n",maxbondperatom);
fprintf(out,"Quantities for this hyper run specific to "
"fix hyper/local:\n");
fprintf(out," ave boost for all bonds/step = %g\n",aveboost);
fprintf(out," ave biased bonds/step = %g\n",avenbias);
fprintf(out," ave bias coeff of all bonds = %g\n",avebiascoeff);
fprintf(out," min bias coeff of any bond = %g\n",minbiascoeff);
@ -382,6 +385,7 @@ void Hyper::dynamics(int nsteps, double & /* time_category */)
lmp->init();
update->integrate->setup(0);
// this may be needed if don't do full init
//modify->addstep_compute_all(update->ntimestep);
bigint ncalls = neighbor->ncalls;
@ -413,7 +417,7 @@ void Hyper::quench(int flag)
update->whichflag = 2;
update->nsteps = maxiter;
update->endstep = update->laststep = update->firststep + maxiter;
update->endstep = update->laststep = update->ntimestep + maxiter;
if (update->laststep < 0)
error->all(FLERR,"Too many iterations");
update->restrict_output = 1;

4
src/fix_store.cpp
View File

@ -280,7 +280,7 @@ int FixStore::pack_restart(int i, double *buf)
return 1;
}
buf[0] = nvalues+1;
buf[0] = ubuf(nvalues+1).d;
if (vecflag) buf[1] = vstore[i];
else
for (int m = 0; m < nvalues; m++)
@ -301,7 +301,7 @@ void FixStore::unpack_restart(int nlocal, int nth)
// skip to Nth set of extra values
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
for (int i = 0; i < nth; i++) m += (int) ubuf(extra[nlocal][m]).i;
m++;
if (vecflag) vstore[nlocal] = extra[nlocal][m];