From fbd610b8a9573bc6ed360ec599a4196c84a109eb Mon Sep 17 00:00:00 2001
From: Steve Plimpton <sjplimp@sandia.gov>
Date: Wed, 10 Oct 2018 13:39:44 -0600
Subject: [PATCH] global/local hyperdynamics src and doc files
---
doc/fix_hyper_global.txt | 146 +++
doc/fix_hyper_local.txt | 200 ++++
doc/hyper.txt | 117 +++
src/REPLICA/fix_event_hyper.cpp | 95 ++
src/REPLICA/fix_event_hyper.h | 60 ++
src/REPLICA/fix_hyper.cpp | 35 +
src/REPLICA/fix_hyper.h | 43 +
src/REPLICA/fix_hyper_global.cpp | 497 +++++++++
src/REPLICA/fix_hyper_global.h | 110 ++
src/REPLICA/fix_hyper_local.cpp | 1610 ++++++++++++++++++++++++++++++
src/REPLICA/fix_hyper_local.h | 175 ++++
src/REPLICA/hyper.cpp | 548 ++++++++++
src/REPLICA/hyper.h | 65 ++
13 files changed, 3701 insertions(+)
create mode 100644 doc/fix_hyper_global.txt
create mode 100644 doc/fix_hyper_local.txt
create mode 100644 doc/hyper.txt
create mode 100644 src/REPLICA/fix_event_hyper.cpp
create mode 100644 src/REPLICA/fix_event_hyper.h
create mode 100644 src/REPLICA/fix_hyper.cpp
create mode 100644 src/REPLICA/fix_hyper.h
create mode 100644 src/REPLICA/fix_hyper_global.cpp
create mode 100644 src/REPLICA/fix_hyper_global.h
create mode 100644 src/REPLICA/fix_hyper_local.cpp
create mode 100644 src/REPLICA/fix_hyper_local.h
create mode 100644 src/REPLICA/hyper.cpp
create mode 100644 src/REPLICA/hyper.h
diff --git a/doc/fix_hyper_global.txt b/doc/fix_hyper_global.txt
new file mode 100644
index 0000000000..c841f0b891
--- /dev/null
+++ b/doc/fix_hyper_global.txt
@@ -0,0 +1,146 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+fix hyper/global command :h3
+
+[Syntax:]
+
+fix ID group-ID hyper/global cutbond qfactor Vmax Tequil :pre
+
+ID, group-ID are documented in "fix"_fix.html command
+hyper/global = style name of this fix command
+cutbond = max distance at which a pair of atoms is considered bonded (distance units)
+qfactor = max strain at which bias potential goes to 0.0 (unitless)
+Vmax = height of bias potential (energy units)
+Tequil = equilibration temperature (temperature units) :ul
+
+[Examples:]
+
+fix 1 all hyper/global 1.0 0.3 0.8 300.0 :pre
+
+[Description:]
+
+This fix is meant to be used with the "hyper"_hyper.html command to
+perform a bond-boost hyperdynamics simulation. The role of this fix
+is a select a single pair of atoms within the system at each timestep
+to add a global bias potential to, which will alter their dynamics.
+This is in contrast to the "fix hyper/local"_fix_hyper_local.html
+command, which can add a local bias potential to multiple pairs of
+atoms at each timestep.
+
+For a system that undergoes rare transition events, where one or more
+atoms move across an energy barrier to a new potential energy basin,
+the effect of the bias potential is to induce more rapid transitions.
+This can lead to a dramatic speed-up in the rate at which events
+occurs, without altering their relative frequencies, thus leading to
+an overall dramatic speed-up in the effective elapsed time of the
+simulation.
+
+Cite various papers.
+
+The current strain of a bond IJ is defined as
+
+Eij = (Rij - R0ij) / R0ij :pre
+
+Emax = is the max of the absolute value of Eij for all IJ bonds.
+
+dVij = Vmax * (1 - (Eij/q)^2) for abs(Eij) < q
+ = 0 otherwise :pre
+
+Delta Vbias = minimum of dVij for all IJ bonds :pre
+
+The boost factor B = exp(beta * Delta Vbias)
+for a single timestep.
+
+Thus the accumulated hypertime is simply
+
+t_hyper = Sum (i = 1 to Nsteps) Bi * dt :pre
+
+The {cutbond} argument is the cutoff distance for defining bonds
+between pairs of nearby atoms. A pair of atoms in their equilibrium,
+minimum-energy config, which are a distance Rij < cutbond, are
+defined as a bonded pair.
+
+The {qfactor} argument is the limiting strain at which
+the Vbias (the bias potential) goes to 0.0.
+
+If {qfactor} is too big, then transitions are affected b/c
+the bias energy is non-zero at the transitions. If it is
+too small than not must boost is achieved for a large system
+with many bonds (some bonds Eij always exceeds qfactor).
+
+The {Vmax} argument is the prefactor on the bias potential.
+
+The {Tequil} argument is part of the beta term in the
+exponential factor that determines how much boost is
+achieved as a function of the bias potential.
+
+[Restart, fix_modify, output, run start/stop, minimize info:]
+
+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 the system's
+potential energy as part of "thermodynamic output"_thermo_style.html.
+
+This fix computes a global scalar and global vector of length 10,
+which can be accessed by various "output
+commands"_Section_howto.html#howto_15. The scalar is the magnitude of
+the bias potential (energy units) applied on the current timestep.
+The vector stores the following quantities:
+
+1 = boost factor on this step (unitless)
+2 = max strain of any bond on this step (unitless)
+3 = ID of first atom in the max-strain bond
+4 = ID of second atom in the max-strain bond
+5 = average # of bonds/atom on this step :ul
+
+6 = fraction of step with no bias during this run
+7 = max drift distance of any atom during this run (distance units)
+8 = max bond length during this run (distance units) :ul
+
+9 = cummulative hyper time since fix created (time units)
+10 = cummulative count of event timesteps since fix created
+11 = cummulative count of atoms in events since fix created :ul
+
+The first 5 quantities are for the current timestep. The quantities
+6-8 are for the current run. The quantities 9-11 are cummulative
+across multiple runs (since the fix was defined in the input script).
+
+For value 10, events are checked for by the "hyper"_hyper.html command
+once every {Nevent} timesteps. This value is the count of the number
+of timesteps on which one (or more) events was detected. It is NOT
+the number of distinct events, since more than one physical event may
+occur in the same {Nevent} time window.
+
+For value 11, each time the "hyper"_hyper.html command checks for an
+event, the event compute it uses will flag zero or more atoms as
+participating in an event. E.g. atoms that have displaced more than
+some distance from the previous quench state. This value is the count
+of the number of atoms participating in any of the events that were
+found.
+
+The scalar and vector values calculated by this fix are all
+"intensive".
+
+No parameter of this fix can be used with the {start/stop} keywords of
+the "run"_run.html command. This fix is not invoked during "energy
+minimization"_minimize.html.
+
+[Restrictions:]
+
+This fix is part of the REPLICA package. It is only enabled if LAMMPS
+was built with that package. See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info.
+
+[Related commands:]
+
+"fix hyper/local"_fix_hyper_local.html
+
+[Default:] None
diff --git a/doc/fix_hyper_local.txt b/doc/fix_hyper_local.txt
new file mode 100644
index 0000000000..ef98637fa5
--- /dev/null
+++ b/doc/fix_hyper_local.txt
@@ -0,0 +1,200 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+fix hyper/local command :h3
+
+[Syntax:]
+
+fix ID group-ID hyper/local cutbond qfactor Vmax Tequil Dcut alpha boost :pre
+
+ID, group-ID are documented in "fix"_fix.html command :ulb,l
+hyper/local = style name of this fix command :l
+cutbond = max distance at which a pair of atoms is considered bonded (distance units) :l
+qfactor = max strain at which bias potential goes to 0.0 (unitless) :l
+Vmax = estimated height of bias potential (energy units) :l
+Tequil = equilibration temperature (temperature units) :l
+Dcut = min distance between boosted bonds (distance units) :l
+alpha = boostostat relaxation time (time units) :l
+boost = desired time boost factor (unitless) :l
+zero or more keyword/value pairs may be appended :l
+keyword = {histo} or {lost} or {check/bias} or {check/coeff}
+ {histo} values = Nevery Nbin delta Nout
+ Nevery = histogram bond bias coefficients every this many timesteps
+ Nbin = # of histogram bins
+ delta = width of each histogram bin
+ Nout = output histogram counts every this many timesteps
+ {lostbond} value = error/warn/ignore
+ {check/bias} values = Nevery error/warn/ignore
+ {check/coeff} values = Nevery error/warn/ignore :pre
+:ule
+
+[Examples:]
+
+fix 1 all hyper/local 1.0 0.3 0.8 300.0 :pre
+
+[Description:]
+
+This fix is meant to be used with the "hyper"_hyper.html command to
+perform a bond-boost hyperdynamics simulation. The role of this fix
+is a select a multiple pairs of atoms within the system at each
+timestep to add a local bias potential to, which will alter their
+dynamics. This is in contrast to the "fix
+hyper/global"_fix_hyper_global.html command, which adds a global bias
+potential to a single pair of atoms at each timestep.
+
+For a system that undergoes rare transition events, where one or more
+atoms move across an energy barrier to a new potential energy basin,
+the effect of the bias potential is to induce more rapid transitions.
+This can lead to a dramatic speed-up in the rate at which events
+occurs, without altering their relative frequencies, thus leading to
+an overall dramatic speed-up in the effective wall-clock time of the
+simulation.
+
+Cite various papers.
+
+The current strain of a bond IJ is defined as
+
+Eij = (Rij - R0ij) / R0ij
+
+Emax = is the max of the absolute value of Eij for all IJ bonds.
+
+dVij = Vmax * (1 - (Eij/q)^2) for abs(Eij) < q
+ = 0 otherwise
+
+Delta Vbias = minimum of dVij for all IJ bonds
+
+The boost factor B = exp(beta * Delta Vbias)
+for a single timestep.
+
+Thus the accumulated hypertime is simply
+
+t_hyper = Sum (i = 1 to Nsteps) Bi * dt
+
+NOTE: Add eqs for boostostat and boost coeff.
+Explain how many bonds are boosted simultaneously
+and how to choose boost factor and initial Vmax.
+
+The {cutbond} argument is the cutoff distance for defining bonds
+between pairs of nearby atoms. A pair of atoms in their equilibrium,
+minimum-energy config, which are a distance Rij < cutbond, are
+defined as a bonded pair.
+
+The {qfactor} argument is the limiting strain at which
+the Vbias (the bias potential) goes to 0.0.
+
+If qfactor is too big, then transitions are affected b/c
+the bias energy is non-zero at the transitions. If it is
+too small than not must boost is achieved for a large system
+with many bonds (some bonds Eij always exceeds qfactor).
+
+The {Vmax} argument is the initial prefactor on the bias potential.
+Should be chosen as estimate of final. Will be adjusted by
+boost cooeficient
+
+The {Tequil} argument is part of the beta term in the
+exponential factor that determines how much boost is
+achieved as a function of the bias potential.
+
+The {Dcut} argument is the distance required between two bonds for
+them to be selected as both being boosted. The distance is between
+the center points of each bond. Actually between any pair of atoms in
+either bond.
+
+The {alpha} argument is a pre-factor on the update equation
+for each atom's boostostat:
+
+NOTE: give equation above
+
+Note that the units for {alpha} are in time units, similar
+to other thermostat or barostat damping parameters
+
+The {boost} argument is the desired boost factor (e.g. 100x)
+that all the atoms in the system will experience.
+
+NOTE: explain how to choose good value for this. If this parameter is
+chosen to be too large, then the bias potential applied to pairs of
+bonded atoms may become unphysically large, leading to bad dynamics.
+If chosen too small, the hyperdynamics run may be inefficient in the
+sense that events take a long time to occur.
+
+[Restart, fix_modify, output, run start/stop, minimize info:]
+
+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 the system's
+potential energy as part of "thermodynamic output"_thermo_style.html.
+
+This fix computes a global scalar and global vector of length 23,
+which can be accessed by various "output
+commands"_Section_howto.html#howto_15. 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 boosted bonds on this step
+2 = max strain of any bond on this step (unitless)
+3 = average bias potential for all bonds on this step (energy units)
+4 = average bonds/atom on this step
+5 = average neighbor bonds/bond on this step :ul
+
+6 = fraction of steps and bonds with no bias during this run
+7 = max drift distance of any atom during this run (distance units)
+8 = max bond length during this run (distance units)
+9 = average # of boosted bonds/step during this run
+10 = average bias potential for all bonds during this run (energy units)
+11 = max bias potential for any bond during this run (energy units)
+12 = min bias potential for any bond during this run (energy units)
+13 = max distance from my box of any ghost atom with maxstrain < qfactor during th
+is run (distance units)
+14 = max distance outside my box of any ghost atom with any maxstrain during this run (distance units)
+15 = count of ghost neighbor atoms not found on reneighbor steps during this run
+16 = count of lost bond partners during this run
+17 = average bias coeff for lost bond partners during this run
+18 = count of bias overlaps found during this run
+19 = count of non-matching bias coefficients found during this run :ul
+
+20 = cummulative hyper time since fix created (time units)
+21 = cummulative count of event timesteps since fix created
+22 = cummulative count of atoms in events since fix created
+23 = cummulative # of new bonds since fix created :ul
+
+The first quantities (1-5) are for the current timestep. The
+quantities 6-19 are for the current hyper run. They are reset each
+time a new hyper run is performed. The quantities 20-23 are
+cummulative across multiple hyper runs, They are only set to initial
+values once, when this fix is defined in the input script.
+
+For value 6, the numerator is a count of all biased bonds on every
+timestep whose bias value = 0.0. The denominator is the count of all
+biased bonds on all timesteps.
+
+For values 13 and 14, the maxstrain of a ghost atom is the maxstrain
+of any bond it is part of, and it is checked for ghost atoms within
+the bond neighbor cutoff.
+
+The scalar and vector values calculated by this fix are all
+"intensive".
+
+No parameter of this fix can be used with the {start/stop} keywords of
+the "run"_run.html command. This fix is not invoked during "energy
+minimization"_minimize.html.
+
+[Restrictions:]
+
+This fix is part of the REPLICA package. It is only enabled if LAMMPS
+was built with that package. See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info.
+
+[Related commands:]
+
+"hyper"_hyper.html, "fix hyper/global"_fix_hyper_global.html
+
+[Default:] None
+
diff --git a/doc/hyper.txt b/doc/hyper.txt
new file mode 100644
index 0000000000..966d33a30f
--- /dev/null
+++ b/doc/hyper.txt
@@ -0,0 +1,117 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+hyper command :h3
+
+[Syntax:]
+
+hyper N Nevent fixID computeID keyword values ... :pre
+
+N = # of timesteps to run :ulb,l
+Nevent = check for events every this many steps :l
+fixID = ID of a fix that applies a global or local bias potential :l
+computeID = ID of a compute that identifies when an event has occurred :l
+zero or more keyword/value pairs may be appended :l
+keyword = {min} or {time} :l
+ {min} values = etol ftol maxiter maxeval
+ etol = stopping tolerance for energy, used in quenching
+ ftol = stopping tolerance for force, used in quenching
+ maxiter = max iterations of minimize, used in quenching
+ maxeval = max number of force/energy evaluations, used in quenching
+ {time} value = {steps} or {clock}
+ {steps} = simulation runs for N timesteps (default)
+ {clock} = simulation runs until hyper time exceeds N timesteps :pre
+:ule
+
+[Examples:]
+
+hyper 5000 100 global event :pre
+
+[Description:]
+
+Run a bond-boost hyperdynamics (HD) simulation where time is
+accelerated by application of a bias potential to a one or more pairs
+of atoms in the system. This command can be used to run both global
+and local hyperdyamics. In global HD a single bond (nearby pair of
+atoms) within the system is biased on a given timestep. In local HD
+multiple bonds (separated by a sufficient distance) can be biased
+simultaneously at each timestep.
+
+Make next paragraph more HD-specific.
+
+Both global and local HD are described in "this paper"_#Voter2013 by
+Art Voter and collaborators. They are methods for performing
+accelerated dynamics that is suitable for infrequent-event systems
+that obey first-order kinetics. A good overview of accelerated
+dynamics methods for such systems in given in "this review
+paper"_#Voter2002prd from the same group. To quote from the paper:
+"The dynamical evolution is characterized by vibrational excursions
+within a potential basin, punctuated by occasional transitions between
+basins." The transition probability is characterized by p(t) =
+k*exp(-kt) where k is the rate constant. Running multiple replicas
+gives an effective enhancement in the timescale spanned by the
+multiple simulations, while waiting for an event to occur.
+
+How different than PRD.
+
+An HD run has several stages, which are repeated each time an "event"
+occurs, as defined below. The logic for a HD run is as follows:
+
+quench
+reset list of bonds :pre
+
+while (time remains):
+ run dynamics for Nevery steps
+ quench
+ check for an event
+ if event occurred: reset list of bonds
+ restore pre-quench state :pre
+
+Explain each of the steps above. Explain what list of bonds is and
+how event is detected.
+
+Explain the specified fix and compute.
+
+Statistics about the number of events, the number of bonds that a bias
+potential is applied to, the accumulated hyper time, etc are stored by
+the fix and can be output with thermodynamic output.
+
+:line
+
+[Restrictions:]
+
+This command can only be used if LAMMPS was built with the REPLICA
+package. See the "Making LAMMPS"_Section_start.html#start_3 section
+for more info on packages.
+
+NOTE: is this true?
+This command cannot be used when any fixes are defined that keep track
+of elapsed time to perform time-dependent operations. Examples
+include the "ave" fixes such as "fix ave/chunk"_fix_ave_chunk.html.
+Also "fix dt/reset"_fix_dt_reset.html and "fix
+deposit"_fix_deposit.html.
+
+[Related commands:]
+
+"fix hyper/global"_fix_hyper_global.html, "fix
+hyper/local"_fix_hyper_local.html, "compute
+event/displace"_compute_event_displace.html
+
+[Default:]
+
+The option defaults are min = 0.1 0.1 40 50 and time = steps.
+
+:line
+
+:link(Voter2013)
+[(Voter2013)] S. Y. Kim, D. Perez, A. F. Voter, J Chem Phys, 139,
+144110 (2013).
+
+:link(Voter2002prd)
+[(Voter2002)] Voter, Montalenti, Germann, Annual Review of Materials
+Research 32, 321 (2002).
diff --git a/src/REPLICA/fix_event_hyper.cpp b/src/REPLICA/fix_event_hyper.cpp
new file mode 100644
index 0000000000..7dd306c735
--- /dev/null
+++ b/src/REPLICA/fix_event_hyper.cpp
@@ -0,0 +1,95 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include <stdlib.h>
+#include <string.h>
+#include "fix_event_hyper.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "neighbor.h"
+#include "comm.h"
+#include "universe.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixEventHyper::FixEventHyper(LAMMPS *lmp, int narg, char **arg) :
+ FixEvent(lmp, narg, arg)
+{
+ if (narg != 3) error->all(FLERR,"Illegal fix event command");
+
+ restart_global = 1;
+
+ event_number = 0;
+ event_timestep = update->ntimestep;
+ clock = 0;
+}
+
+/* ----------------------------------------------------------------------
+ save current atom coords as an event (via call to base class)
+ called when an event occurs in some replica
+ set event_timestep = when event occurred in a particular replica
+ update clock = elapsed time since last event, across all replicas
+------------------------------------------------------------------------- */
+
+void FixEventHyper::store_event_hyper(bigint ntimestep, int delta_clock)
+{
+ store_event();
+ event_timestep = ntimestep;
+ clock += delta_clock;
+ event_number++;
+}
+
+/* ----------------------------------------------------------------------
+ pack entire state of Fix into one write
+------------------------------------------------------------------------- */
+
+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;
+
+ if (comm->me == 0) {
+ int size = n * sizeof(double);
+ fwrite(&size,sizeof(int),1,fp);
+ fwrite(list,sizeof(double),n,fp);
+ }
+}
+
+/* ----------------------------------------------------------------------
+ use state info from restart file to restart the Fix
+------------------------------------------------------------------------- */
+
+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++]);
+}
diff --git a/src/REPLICA/fix_event_hyper.h b/src/REPLICA/fix_event_hyper.h
new file mode 100644
index 0000000000..4c5d4a93ee
--- /dev/null
+++ b/src/REPLICA/fix_event_hyper.h
@@ -0,0 +1,60 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(EVENT/HYPER,FixEventHyper)
+
+#else
+
+#ifndef LMP_FIX_EVENT_HYPER_H
+#define LMP_FIX_EVENT_HYPER_H
+
+#include "fix_event.h"
+
+namespace LAMMPS_NS {
+
+class FixEventHyper : public FixEvent {
+ public:
+ int event_number; // event counter
+ bigint event_timestep; // timestep of last event on any replica
+ bigint clock; // total elapsed timesteps across all replicas
+ int replica_number; // replica where last event occured
+ int correlated_event; // 1 if last event was correlated, 0 otherwise
+ int ncoincident; // # of simultaneous events on different replicas
+
+ FixEventHyper(class LAMMPS *, int, char **);
+ ~FixEventHyper() {}
+
+ void write_restart(FILE *);
+ void restart(char *);
+
+ // methods specific to FixEventHyper, invoked by hyper
+
+ void store_event_hyper(bigint, int);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+*/
diff --git a/src/REPLICA/fix_hyper.cpp b/src/REPLICA/fix_hyper.cpp
new file mode 100644
index 0000000000..9f01d9542c
--- /dev/null
+++ b/src/REPLICA/fix_hyper.cpp
@@ -0,0 +1,35 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include <string.h>
+#include "fix_hyper.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+FixHyper::FixHyper(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) {}
+
+/* ----------------------------------------------------------------------
+ extract hyper flag setting for all Fixes that perform hyperdynamics
+------------------------------------------------------------------------- */
+
+void *FixHyper::extract(const char *str, int &dim)
+{
+ dim = 0;
+ if (strcmp(str,"hyperflag") == 0) {
+ return &hyperflag;
+ }
+ return NULL;
+}
+
diff --git a/src/REPLICA/fix_hyper.h b/src/REPLICA/fix_hyper.h
new file mode 100644
index 0000000000..7cb6cadf8f
--- /dev/null
+++ b/src/REPLICA/fix_hyper.h
@@ -0,0 +1,43 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef LMP_FIX_HYPER_H
+#define LMP_FIX_HYPER_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixHyper : public Fix {
+ public:
+ FixHyper(class LAMMPS *, int, char **);
+ virtual ~FixHyper() {}
+ void *extract(const char *, int &);
+
+ // must be provided by child class
+
+ virtual void init_hyper() = 0;
+ virtual void build_bond_list(int) = 0;
+ virtual double query(int) = 0;
+
+ protected:
+ int hyperflag;
+};
+
+}
+
+#endif
+
+/* ERROR/WARNING messages:
+
+*/
diff --git a/src/REPLICA/fix_hyper_global.cpp b/src/REPLICA/fix_hyper_global.cpp
new file mode 100644
index 0000000000..ba45e47407
--- /dev/null
+++ b/src/REPLICA/fix_hyper_global.cpp
@@ -0,0 +1,497 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include <mpi.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#include "fix_hyper_global.h"
+#include "atom.h"
+#include "update.h"
+#include "force.h"
+#include "domain.h"
+#include "comm.h"
+#include "neighbor.h"
+#include "neigh_request.h"
+#include "neigh_list.h"
+#include "modify.h"
+#include "math_extra.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+#define DELTA 16384
+#define VECLEN 5
+
+/* ---------------------------------------------------------------------- */
+
+FixHyperGlobal::FixHyperGlobal(LAMMPS *lmp, int narg, char **arg) :
+ FixHyper(lmp, narg, arg)
+{
+ // NOTE: add NULL declarations
+ // NOTE: count/output # of timesteps on which bias is non-zero
+ // NOTE: require newton on?
+
+ if (atom->map_style == 0)
+ error->all(FLERR,"Fix hyper/global command requires atom map");
+
+ if (narg != 7) error->all(FLERR,"Illegal fix hyper/global command");
+
+ hyperflag = 1;
+ scalar_flag = 1;
+ vector_flag = 1;
+ size_vector = 11;
+ global_freq = 1;
+ extscalar = 0;
+ extvector = 0;
+
+ cutbond = force->numeric(FLERR,arg[3]);
+ qfactor = force->numeric(FLERR,arg[4]);
+ vmax = force->numeric(FLERR,arg[5]);
+ tequil = force->numeric(FLERR,arg[6]);
+
+ if (cutbond < 0.0 || qfactor <= 0.0 || vmax < 0.0 || tequil <= 0.0)
+ error->all(FLERR,"Illegal fix hyper/global command");
+
+ invqfactorsq = 1.0 / (qfactor*qfactor);
+ cutbondsq = cutbond*cutbond;
+ beta = 1.0 / (force->boltz * tequil);
+
+ maxbond = 0;
+ nblocal = 0;
+ blist = NULL;
+
+ maxold = 0;
+ xold = NULL;
+ tagold = NULL;
+
+ me = comm->me;
+ firstflag = 1;
+ bcastflag = 0;
+ for (int i = 0; i < VECLEN; i++) outvec[i] = 0.0;
+
+ nevent = 0;
+ nevent_atom = 0;
+ t_hyper = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixHyperGlobal::~FixHyperGlobal()
+{
+ memory->sfree(blist);
+ memory->destroy(xold);
+ memory->destroy(tagold);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixHyperGlobal::setmask()
+{
+ int mask = 0;
+ mask |= PRE_NEIGHBOR;
+ mask |= PRE_FORCE;
+ mask |= PRE_REVERSE;
+ mask |= THERMO_ENERGY;
+ return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperGlobal::init_hyper()
+{
+ maxdriftsq = 0.0;
+ maxbondlen = 0.0;
+ nobias = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperGlobal::init()
+{
+ dt = update->dt;
+
+ // need an occasional half neighbor list
+
+ int irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->pair = 0;
+ neighbor->requests[irequest]->fix = 1;
+ neighbor->requests[irequest]->occasional = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperGlobal::init_list(int id, NeighList *ptr)
+{
+ list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperGlobal::setup_pre_neighbor()
+{
+ pre_neighbor();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperGlobal::setup_pre_reverse(int eflag, int vflag)
+{
+ // no increment in nobias or hyper time when pre-run forces are calculated
+
+ int nobias_hold = nobias;
+ double t_hyper_hold = t_hyper;
+
+ pre_reverse(eflag,vflag);
+
+ nobias = nobias_hold;
+ t_hyper = t_hyper_hold;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperGlobal::pre_neighbor()
+{
+ int i,m,jnum,iold,jold,ilocal,jlocal;
+ double distsq;
+
+ // reset local IDs for owned bond atoms, since atoms have migrated
+ // uses xold and tagold from when bonds were created
+ // must be done after ghost atoms are setup via comm->borders()
+
+ double **x = atom->x;
+
+ int flag = 0;
+
+ for (m = 0; m < nblocal; m++) {
+ iold = blist[m].iold;
+ jold = blist[m].jold;
+ ilocal = atom->map(tagold[iold]);
+ jlocal = atom->map(tagold[jold]);
+ ilocal = domain->closest_image(xold[iold],ilocal);
+ jlocal = domain->closest_image(xold[iold],jlocal);
+ blist[m].i = ilocal;
+ blist[m].j = jlocal;
+
+ if (ilocal < 0 || jlocal < 0) flag++;
+ else {
+ distsq = MathExtra::distsq3(x[ilocal],xold[iold]);
+ maxdriftsq = MAX(distsq,maxdriftsq);
+ }
+ }
+
+ if (flag) error->one(FLERR,"Fix hyper/global bond atom not found");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperGlobal::pre_reverse(int eflag, int vflag)
+{
+ int i,j,m,imax,jmax;
+ double xtmp,ytmp,ztmp,delx,dely,delz;
+ double r,r0,estrain,rmax,r0max,emax,dt_boost;
+ double vbias,fbias,fbiasr;
+
+ // compute current strain of each owned bond
+ // emax = maximum strain of any bond I own
+ // imax,jmax = local indices of my 2 atoms in that bond
+
+ double **x = atom->x;
+ emax = 0.0;
+
+ for (m = 0; m < nblocal; m++) {
+ i = blist[m].i;
+ j = blist[m].j;
+ delx = x[i][0] - x[j][0];
+ dely = x[i][1] - x[j][1];
+ delz = x[i][2] - x[j][2];
+ r = sqrt(delx*delx + dely*dely + delz*delz);
+ maxbondlen = MAX(r,maxbondlen);
+ r0 = blist[m].r0;
+ estrain = fabs(r-r0) / r0;
+
+ if (estrain > emax) {
+ emax = estrain;
+ rmax = r;
+ r0max = r0;
+ imax = i;
+ jmax = j;
+ }
+ }
+
+ // perform Allreduce() on pairme = double/int pair
+ // finds max strain and what proc owns it
+ // owner = proc that owns that bond
+
+ pairme.value = emax;
+ pairme.proc = me;
+ MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MAXLOC,world);
+ owner = pairall.proc;
+
+ bcastflag = 1;
+
+ // all procs acquire t_hyper from owner proc
+ // MPI_Bcast call by owner proc is below
+
+ for (int i = 0; i < VECLEN; i++) outvec[i] = 0.0;
+
+ if (owner != me) {
+ MPI_Bcast(&t_hyper,1,MPI_DOUBLE,owner,world);
+ return;
+ }
+
+ // I own the bond with max strain
+ // compute Vbias and apply force to atoms imax,jmax
+ // NOTE: logic would need to be different for newton off
+
+ double **f = atom->f;
+
+ vbias = fbias = 0.0;
+ dt_boost = 1.0;
+
+ if (emax < qfactor) {
+ vbias = vmax * (1.0 - emax*emax*invqfactorsq);
+ fbias = 2.0 * vmax * emax / (qfactor*qfactor * r0max);
+ dt_boost = exp(beta*vbias);
+
+ delx = x[imax][0] - x[jmax][0];
+ dely = x[imax][1] - x[jmax][1];
+ delz = x[imax][2] - x[jmax][2];
+
+ fbiasr = fbias / rmax;
+ f[imax][0] += delx*fbiasr;
+ f[imax][1] += dely*fbiasr;
+ f[imax][2] += delz*fbiasr;
+
+ f[jmax][0] -= delx*fbiasr;
+ f[jmax][1] -= dely*fbiasr;
+ f[jmax][2] -= delz*fbiasr;
+ } else nobias++;
+
+ // NOTE: should there be a virial contribution from boosted bond?
+
+ // output quantities
+
+ outvec[0] = vbias;
+ outvec[1] = dt_boost;
+ outvec[2] = emax;
+ outvec[3] = atom->tag[imax];
+ outvec[4] = atom->tag[jmax];
+
+ t_hyper += dt_boost*dt;
+ MPI_Bcast(&t_hyper,1,MPI_DOUBLE,owner,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperGlobal::build_bond_list(int natom)
+{
+ int i,j,n,ii,jj,inum,jnum;
+ double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ tagint *bptr;
+ double *rptr;
+
+ if (natom) {
+ nevent++;
+ nevent_atom += natom;
+ }
+
+ // trigger neighbor list build
+
+ neighbor->build_one(list);
+
+ // identify bonds assigned to each owned atom
+ // do not create a bond between two non-group atoms
+
+ double **x = atom->x;
+ int *mask = atom->mask;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ nblocal = 0;
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+
+ // skip if neither atom I or J are in fix group
+
+ if (!(mask[i] & groupbit) && !(mask[j] & groupbit)) continue;
+
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cutbondsq) {
+ if (nblocal == maxbond) grow_bond();
+ blist[nblocal].i = i;
+ blist[nblocal].j = j;
+ blist[nblocal].iold = i;
+ blist[nblocal].jold = j;
+ blist[nblocal].r0 = sqrt(rsq);
+ nblocal++;
+ }
+ }
+ }
+
+ // store IDs and coords for owned+ghost atoms at time of bond creation
+ // realloc xold and tagold as needed
+
+ if (atom->nmax > maxold) {
+ memory->destroy(xold);
+ memory->destroy(tagold);
+ maxold = atom->nmax;
+ memory->create(xold,maxold,3,"hyper/global:xold");
+ memory->create(tagold,maxold,"hyper/global:tagold");
+ }
+
+ tagint *tag = atom->tag;
+ int nall = atom->nlocal + atom->nghost;
+
+ for (i = 0; i < nall; i++) {
+ xold[i][0] = x[i][0];
+ xold[i][1] = x[i][1];
+ xold[i][2] = x[i][2];
+ tagold[i] = tag[i];
+ }
+}
+
+/* ----------------------------------------------------------------------
+ grow bond list by a chunk
+------------------------------------------------------------------------- */
+
+void FixHyperGlobal::grow_bond()
+{
+ // NOTE: could add int arg to do initial large alloc:
+ // maxbond = maxbond/DELTA * DELTA; maxbond += DELTA;
+
+ maxbond += DELTA;
+ if (maxbond < 0 || maxbond > MAXSMALLINT)
+ error->one(FLERR,"Fix hyper/local per-processor bond count is too big");
+ blist = (OneBond *)
+ memory->srealloc(blist,maxbond*sizeof(OneBond),"hyper/local:blist");
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixHyperGlobal::compute_scalar()
+{
+ if (bcastflag) {
+ MPI_Bcast(outvec,VECLEN,MPI_DOUBLE,owner,world);
+ bcastflag = 0;
+ }
+ return outvec[0];
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixHyperGlobal::compute_vector(int i)
+{
+ if (bcastflag) {
+ MPI_Bcast(outvec,VECLEN,MPI_DOUBLE,owner,world);
+ bcastflag = 0;
+ }
+
+ // 11 vector outputs returned for i = 0-10
+
+ // i = 0 = boost factor on this step
+ // i = 1 = max strain of any bond on this step
+ // i = 2 = ID of atom I in max-strain bond on this step
+ // i = 3 = ID of atom J in max-strain bond on this step
+ // i = 4 = ave bonds/atom on this step
+
+ // i = 5 = fraction of steps with no bias during this run
+ // i = 6 = max drift of any atom during this run
+ // i = 7 = max bond length during this run
+
+ // i = 8 = cummulative hyper time since fix created
+ // i = 9 = cummulative # of event timesteps since fix created
+ // i = 10 = cummulative # of atoms in events since fix created
+
+ if (i == 0) return outvec[1];
+ if (i == 1) return outvec[2];
+ if (i == 2) return outvec[3];
+ if (i == 3) return outvec[4];
+
+ if (i == 4) {
+ int allbonds; // NOTE: bigint?
+ MPI_Allreduce(&nblocal,&allbonds,1,MPI_INT,MPI_SUM,world);
+ return 2.0*allbonds/atom->natoms;
+ }
+
+ if (i == 5) {
+ if (update->ntimestep == update->firststep) return 0.0;
+ int allnobias;
+ MPI_Allreduce(&nobias,&allnobias,1,MPI_INT,MPI_SUM,world);
+ return 1.0*allnobias / (update->ntimestep - update->firststep);
+ }
+
+ if (i == 6) {
+ double alldriftsq;
+ MPI_Allreduce(&maxdriftsq,&alldriftsq,1,MPI_DOUBLE,MPI_MAX,world);
+ return sqrt(alldriftsq);
+ }
+
+ if (i == 7) {
+ double allbondlen;
+ MPI_Allreduce(&maxbondlen,&allbondlen,1,MPI_DOUBLE,MPI_MAX,world);
+ return allbondlen;
+ }
+
+ if (i == 8) return t_hyper;
+ if (i == 9) return (double) nevent;
+ if (i == 10) return (double) nevent_atom;
+}
+
+/* ----------------------------------------------------------------------
+ wrapper on compute_vector()
+ used by hyper.cpp to call FixHyper
+------------------------------------------------------------------------- */
+
+double FixHyperGlobal::query(int i)
+{
+ if (i == 1) return compute_vector(8); // cummulative hyper time
+ if (i == 2) return compute_vector(9); // nevent
+ if (i == 3) return compute_vector(10); // nevent_atom
+ if (i == 4) return compute_vector(4); // ave bonds/atom
+ if (i == 5) return compute_vector(6); // maxdrift
+ if (i == 6) return compute_vector(7); // maxbondlen
+ if (i == 7) return compute_vector(5); // fraction with zero bias
+
+ error->all(FLERR,"Invalid query to fix hyper/global");
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+
+double FixHyperGlobal::memory_usage()
+{
+ double bytes = maxbond * sizeof(OneBond); // blist
+ return bytes;
+}
diff --git a/src/REPLICA/fix_hyper_global.h b/src/REPLICA/fix_hyper_global.h
new file mode 100644
index 0000000000..d8b1cef425
--- /dev/null
+++ b/src/REPLICA/fix_hyper_global.h
@@ -0,0 +1,110 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(hyper/global,FixHyperGlobal)
+
+#else
+
+#ifndef LMP_FIX_HYPER_GLOBAL_H
+#define LMP_FIX_HYPER_GLOBAL_H
+
+#include "fix_hyper.h"
+
+namespace LAMMPS_NS {
+
+class FixHyperGlobal : public FixHyper {
+ public:
+ FixHyperGlobal(class LAMMPS *, int, char **);
+ ~FixHyperGlobal();
+ int setmask();
+ void init();
+ void init_list(int, class NeighList *);
+ void setup_pre_neighbor();
+ void setup_pre_reverse(int, int);
+ void pre_neighbor();
+ void pre_reverse(int, int);
+ double compute_scalar();
+ double compute_vector(int);
+ double query(int);
+
+ double memory_usage();
+
+ // extra methods visible to callers
+
+ void init_hyper();
+ void build_bond_list(int);
+
+ private:
+ int me;
+ double cutbond,qfactor,vmax,tequil;
+
+ int firstflag,bcastflag,owner,nevent,nevent_atom;
+ double cutbondsq,beta,dt,t_hyper,invqfactorsq;
+ double outvec[5]; // same as VECLEN in *.cpp
+ double maxbondlen; // max length of any bond
+ double maxdriftsq; // max distance any atom drifts from original pos
+ int nobias; // # of steps when bias = 0, b/c bond too long
+
+ class NeighList *list;
+
+ // list of my owned bonds
+ // persists on a proc from one event until the next
+
+ int maxbond; // allocated size of blist
+
+ struct OneBond { // single IJ bond, atom I is owner
+ int i,j; // current local indices of 2 bond atoms
+ int iold,jold; // local indices when bonds were formed
+ double r0; // relaxed bond length
+ };
+
+ struct OneBond *blist; // list of owned bonds
+ int nblocal; // # of owned bonds
+
+ // coords and IDs of owned+ghost atoms when bonds were formed
+ // persists on a proc from one event until the next
+
+ int maxold; // allocated size of old atoms
+
+ double **xold; // coords of atoms when bonds were formed
+ tagint *tagold; // IDs of atoms when bonds were formed
+
+ // MPI data struct for finding bond with max strain via Allreduce
+
+ struct Two {
+ double value;
+ int proc;
+ };
+ Two pairme,pairall;
+
+ // internal methods
+
+ void grow_bond();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+*/
diff --git a/src/REPLICA/fix_hyper_local.cpp b/src/REPLICA/fix_hyper_local.cpp
new file mode 100644
index 0000000000..b60ba7339f
--- /dev/null
+++ b/src/REPLICA/fix_hyper_local.cpp
@@ -0,0 +1,1610 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU Gene<ral Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include <mpi.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#include "fix_hyper_local.h"
+#include "atom.h"
+#include "update.h"
+#include "force.h"
+#include "pair.h"
+#include "domain.h"
+#include "comm.h"
+#include "neighbor.h"
+#include "neigh_request.h"
+#include "neigh_list.h"
+#include "modify.h"
+#include "math_extra.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+#define DELTABOOST 16
+#define BOOSTINIT 1.0
+#define COEFFMAX 1.2
+#define BIG 1.0e20
+
+enum{STRAIN,STRAINREGION,BIASFLAG};
+enum{IGNORE,WARN,ERROR};
+
+/* ---------------------------------------------------------------------- */
+
+FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) :
+ FixHyper(lmp, narg, arg)
+{
+ // NOTE: add NULL declarations
+
+ // error checks
+ // solution for tagint != int is to worry about storing
+ // local index vs global ID in same variable
+ // maybe need to declare them all tagint, not int
+
+ if (atom->map_style == 0)
+ error->all(FLERR,"Fix hyper/local command requires atom map");
+
+ if (sizeof(tagint) != sizeof(int))
+ error->all(FLERR,"Fix hyper/local requires tagint = int");
+
+ // parse args
+
+ if (narg < 10) error->all(FLERR,"Illegal fix hyper/local command");
+
+ hyperflag = 2;
+ scalar_flag = 1;
+ vector_flag = 1;
+ size_vector = 23;
+
+ global_freq = 1;
+ extscalar = 0;
+ extvector = 0;
+
+ cutbond = force->numeric(FLERR,arg[3]);
+ qfactor = force->numeric(FLERR,arg[4]);
+ vmax = force->numeric(FLERR,arg[5]);
+ tequil = force->numeric(FLERR,arg[6]);
+ dcut = force->numeric(FLERR,arg[7]);
+ alpha_user = force->numeric(FLERR,arg[8]);
+ boosttarget = force->numeric(FLERR,arg[9]);
+
+ if (cutbond < 0.0 || qfactor < 0.0 || vmax < 0.0 ||
+ tequil <= 0.0 || dcut <= 0.0 || alpha_user <= 0.0 || boosttarget < 1.0)
+ error->all(FLERR,"Illegal fix hyper/local command");
+
+ invqfactorsq = 1.0 / (qfactor*qfactor);
+ cutbondsq = cutbond*cutbond;
+ dcutsq = dcut*dcut;
+ beta = 1.0 / (force->boltz * tequil);
+
+ // optional args
+
+ histoflag = 0;
+ lostbond = IGNORE;
+ checkbias = 0;
+ checkcoeff = 0;
+
+ int iarg = 10;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"histo") == 0) {
+ if (iarg+5 > narg) error->all(FLERR,"Illegal fix hyper/local command");
+ histoflag = 1;
+ histo_every = force->inumeric(FLERR,arg[iarg+1]);
+ histo_count = force->inumeric(FLERR,arg[iarg+2]);
+ histo_delta = force->numeric(FLERR,arg[iarg+3]);
+ histo_print = force->inumeric(FLERR,arg[iarg+4]);
+ if (histo_every <= 0 || histo_count % 2 ||
+ histo_delta <= 0.0 || histo_print <= 0)
+ error->all(FLERR,"Illegal fix hyper/local command");
+ iarg += 5;
+
+ } else if (strcmp(arg[iarg],"lostbond") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal fix hyper/local command");
+ if (strcmp(arg[iarg+1],"error") == 0) lostbond = ERROR;
+ else if (strcmp(arg[iarg+1],"warn") == 0) lostbond = WARN;
+ else if (strcmp(arg[iarg+1],"ignore") == 0) lostbond = IGNORE;
+ else error->all(FLERR,"Illegal fix hyper/local command");
+ iarg += 2;
+
+ } else if (strcmp(arg[iarg],"check/bias") == 0) {
+ if (iarg+3 > narg) error->all(FLERR,"Illegal fix hyper/local command");
+ checkbias = 1;
+ checkbias_every = force->inumeric(FLERR,arg[iarg+1]);
+ if (strcmp(arg[iarg+2],"error") == 0) checkbias_flag = ERROR;
+ else if (strcmp(arg[iarg+2],"warn") == 0) checkbias_flag = WARN;
+ else if (strcmp(arg[iarg+2],"ignore") == 0) checkbias_flag = IGNORE;
+ else error->all(FLERR,"Illegal fix hyper/local command");
+ iarg += 3;
+
+ } else if (strcmp(arg[iarg],"check/coeff") == 0) {
+ if (iarg+3 > narg) error->all(FLERR,"Illegal fix hyper/local command");
+ checkcoeff = 1;
+ checkcoeff_every = force->inumeric(FLERR,arg[iarg+1]);
+ if (strcmp(arg[iarg+2],"error") == 0) checkcoeff_flag = ERROR;
+ else if (strcmp(arg[iarg+2],"warn") == 0) checkcoeff_flag = WARN;
+ else if (strcmp(arg[iarg+2],"ignore") == 0) checkcoeff_flag = IGNORE;
+ else error->all(FLERR,"Illegal fix hyper/local command");
+ iarg += 3;
+
+ } else error->all(FLERR,"Illegal fix hyper/local command");
+ }
+
+ // per-atom data structs
+
+ maxbond = 0;
+ bonds = NULL;
+ numbond = NULL;
+ maxstrain = NULL;
+ maxstrain_region = NULL;
+ maxstrain_bondindex = NULL;
+ biasflag = NULL;
+
+ maxold = old_nall = 0;
+ old2now = NULL;
+ xold = NULL; // NOTE: don't really need except to monitor drift
+ tagold = NULL;
+
+ nboost = maxboost = 0;
+ boost = NULL;
+
+ // maxbondperatom = max # of bonds any atom is part of
+ // will be reset in bond_build()
+ // set comm size needed by this fix
+
+ maxbondperatom = 1;
+ comm_forward = 1;
+ comm_reverse = 1;
+
+ // perform initial allocation of atom-based arrays
+ // register with Atom class
+
+ grow_arrays(atom->nmax);
+ atom->add_callback(0);
+
+ me = comm->me;
+ firstflag = 1;
+
+ allbonds = 0;
+ allboost = 0.0;
+
+ starttime = update->ntimestep;
+ nostrainyet = 1;
+ nnewbond = 0;
+ nevent = 0;
+ nevent_atom = 0;
+ mybias = 0.0;
+
+ histo = allhisto = NULL;
+ if (histoflag) {
+ invhisto_delta = 1.0 / histo_delta;
+ histo_lo = 1.0 - (histo_count/2 * histo_delta);
+ histo = new bigint[histo_count+2];
+ allhisto = new bigint[histo_count+2];
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixHyperLocal::~FixHyperLocal()
+{
+ memory->destroy(bonds);
+ memory->destroy(numbond);
+
+ memory->destroy(maxstrain);
+ memory->destroy(maxstrain_region);
+ memory->destroy(maxstrain_bondindex);
+ memory->destroy(biasflag);
+
+ memory->destroy(old2now);
+ memory->destroy(xold);
+ memory->destroy(tagold);
+ memory->destroy(boost);
+ delete [] histo;
+ delete [] allhisto;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixHyperLocal::setmask()
+{
+ int mask = 0;
+ mask |= PRE_NEIGHBOR;
+ mask |= PRE_REVERSE;
+ mask |= MIN_PRE_NEIGHBOR;
+ mask |= THERMO_ENERGY;
+ return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperLocal::init_hyper()
+{
+ ghost_toofar = 0;
+ lostbond_partner = 0;
+ lostbond_coeff = 0.0;
+ checkbias_count = 0;
+ checkcoeff_count = 0;
+ maxdriftsq = 0.0;
+ maxbondlen = 0.0;
+ maxboostcoeff = 0.0;
+ minboostcoeff = BIG;
+ sumboostcoeff = 0.0;
+ nboost_running = 0;
+ nobias_running = 0;
+ rmaxever = 0.0;
+ rmaxeverbig = 0.0;
+
+ nbondbuild = 0;
+ time_bondbuild = 0.0;
+
+ if (histoflag) histo_steps = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperLocal::init()
+{
+ alpha = update->dt / alpha_user;
+
+ // cutghost = communication cutoff as calculated by Neighbor and Comm
+ // error if cutghost is smaller than Dcut
+ // warn if no drift distance added to cutghost
+
+ if (firstflag) {
+ double cutghost;
+ if (force->pair)
+ cutghost = MAX(force->pair->cutforce+neighbor->skin,comm->cutghostuser);
+ else
+ cutghost = comm->cutghostuser;
+
+ if (cutghost < dcut)
+ error->all(FLERR,"Fix hyper/local bond cutoff exceeds ghost atom range - "
+ "use comm_modify cutoff command");
+ if (cutghost < dcut+cutbond && me == 0)
+ error->warning(FLERR,"Fix hyper/local ghost atom range "
+ "may not allow for atom drift between events");
+ }
+
+ // NOTE: error if not newton,
+ // since bond force bias will not be applied correctly to straddle bonds?
+
+ // NOTE: if exclusions are enabled, some near-neighbors won't appear
+ // will be missed in bonds - WARN if molecular system?
+
+ // need an occasional full neighbor list with cutoff = Dcut
+ // do not need to include neigh skin in cutoff,
+ // b/c this list will be built every time bond_build() is called
+ // NOTE: what if pair style list cutoff > Dcut
+ // or what if neigh skin is huge?
+
+ int irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->pair = 0;
+ neighbor->requests[irequest]->fix = 1;
+ neighbor->requests[irequest]->half = 0;
+ neighbor->requests[irequest]->full = 1;
+ neighbor->requests[irequest]->cut = 1;
+ neighbor->requests[irequest]->cutoff = dcut;
+ neighbor->requests[irequest]->occasional = 1;
+
+ // DEBUG timing output
+
+ timefirst = timesecond = timethird = timefourth = timefifth =
+ timesixth = timeseventh = timetotal = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperLocal::init_list(int id, NeighList *ptr)
+{
+ list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperLocal::setup_pre_neighbor()
+{
+ // called for dynamics and minimization NOTE: check if min is needed?
+
+ pre_neighbor();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperLocal::setup_pre_reverse(int eflag, int vflag)
+{
+ // only called for dynamics, not minimization
+ // setupflag prevents boostostat update of boost coeffs in setup
+ // also prevents increments of nboost_running, nbias_running, sumboostcoeff
+
+ setupflag = 1;
+ pre_reverse(eflag,vflag);
+ setupflag = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperLocal::pre_neighbor()
+{
+ int i,m,n,ilocal,jlocal;
+
+ // convert global ID bond partners back to local indices
+ // need to use closest_image() so can calculate bond lengths
+ // error flag should not happen for a well-behaved system
+ // b/c are only looking up bond partners inside or near my sub-domain
+
+ double **x = atom->x;
+ int nlocal = atom->nlocal;
+
+ int missing = 0;
+ double missing_coeff = 0.0;
+
+ for (i = 0; i < nlocal; i++) {
+ n = numbond[i];
+ for (m = 0; m < n; m++) {
+ jlocal = atom->map(bonds[i][m].jtag);
+ if (jlocal >= 0) bonds[i][m].j = domain->closest_image(i,jlocal);
+ else {
+ bonds[i][m].j = -1;
+ missing++;
+ missing_coeff += bonds[i][m].boostcoeff;
+ if (lostbond != IGNORE) {
+ char str[128];
+ sprintf(str,"Fix hyper/local bond info missing for bond "
+ TAGINT_FORMAT "," TAGINT_FORMAT
+ " with coeff %g at step " BIGINT_FORMAT,
+ atom->tag[i],bonds[i][m].jtag,bonds[i][m].boostcoeff,
+ update->ntimestep);
+ if (lostbond == ERROR) error->one(FLERR,str);
+ if (lostbond == WARN) error->warning(FLERR,str);
+ }
+ }
+ }
+ }
+
+ lostbond_partner += missing;
+ lostbond_coeff += missing_coeff;
+
+ // set old2now to point to current local atom indices
+ // only necessary for tagold entries > 0
+ // because if tagold = 0, atom is not active in Dcut neighbor list
+ // must be done after atoms migrate and ghost atoms setup via comm->borders()
+ // does not matter if there are multiple ghost copies of a global ID
+ // since only need the ghost atom strain, not its coordinates
+ // NOTE: maybe need not use closest image, b/c old2now only used to access
+ // maxstrain values, which will be same for any copy of ghost atom ??
+ // also need it to compute maxdriftsq correctly when a proc spans a PBC
+
+ double distsq;
+
+ for (i = 0; i < old_nall; i++) {
+ if (tagold[i] == 0) continue;
+ ilocal = atom->map(tagold[i]);
+ ilocal = domain->closest_image(xold[i],ilocal);
+ old2now[i] = ilocal;
+
+ if (ilocal >= 0) {
+ distsq = MathExtra::distsq3(x[ilocal],xold[i]);
+ maxdriftsq = MAX(distsq,maxdriftsq);
+ } else ghost_toofar++;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperLocal::pre_reverse(int eflag, int vflag)
+{
+ int i,j,m,ii,jj,inum,jnum,iold,jold,nbond,bondindex;
+ tagint itag,jtag;
+ double xtmp,ytmp,ztmp,delx,dely,delz;
+ double r,r0,estrain,emax,vbias,fbias,fbiasr,boostcoeff;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ //double time1,time2,time3,time4,time5,time6,time7,time8;
+ //time1 = MPI_Wtime();
+
+ // compute current maxstrain and maxstrain_bond for each owned atom
+ // use per-atom full bond list
+ // this is double-calculating for IJ and JI bonds
+ // could compute once, but would have to find/store index of JI bond
+ // order two I,J atoms consistently for IJ and JI calcs
+ // to insure no round-off issue when comparing maxstrain values of I,J
+
+ double **x = atom->x;
+ tagint *tag = atom->tag;
+ int nlocal = atom->nlocal;
+
+ int mybonds = 0;
+ nostrainyet = 0;
+
+ for (i = 0; i < nlocal; i++) {
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itag = tag[i];
+ emax = 0.0;
+ bondindex = -1;
+ jtag = 0;
+
+ nbond = numbond[i];
+ mybonds += nbond;
+ for (m = 0; m < nbond; m++) {
+ j = bonds[i][m].j;
+ if (j < 0) continue;
+ jtag = bonds[i][m].jtag;
+ r0 = bonds[i][m].r0;
+ if (itag < jtag) {
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ } else {
+ delx = x[j][0] - xtmp;
+ dely = x[j][1] - ytmp;
+ delz = x[j][2] - ztmp;;
+ }
+ r = sqrt(delx*delx + dely*dely + delz*delz);
+ maxbondlen = MAX(r,maxbondlen);
+ estrain = fabs(r-r0) / r0;
+ if (estrain > emax) {
+ emax = estrain;
+ bondindex = m;
+ }
+ }
+ maxstrain[i] = emax;
+ maxstrain_bondindex[i] = bondindex;
+ }
+
+ //time2 = MPI_Wtime();
+
+ // forward comm to acquire maxstrain of all ghost atoms
+
+ commflag = STRAIN;
+ comm->forward_comm_fix(this);
+
+ //time3 = MPI_Wtime();
+
+ // use original Dcut neighbor list to check maxstrain of all neighbor atoms
+ // set maxstrain_region of I atoms = maxstrain of I and all J neighs
+ // neighbor list has old indices for IJ b/c reneighboring may have occurred
+ // use old2now[] to convert to current indices
+ // if neighbor is not currently known (too far away),
+ // then assume it was part of an event and its strain = qfactor
+ // this double loop sets maxstrain_region of mostly owned atoms
+ // but possibly some ghost atoms as well
+
+ int nall = nlocal + atom->nghost;
+ for (i = 0; i < nall; i++) maxstrain_region[i] = 0.0;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ // find largest distance from subbox that a ghost atom is with strain < qfactor
+
+ double rmax = rmaxever;
+ double rmaxbig = rmaxeverbig;
+ double *sublo = domain->sublo;
+ double *subhi = domain->subhi;
+
+ for (ii = 0; ii < inum; ii++) {
+ iold = ilist[ii];
+ jlist = firstneigh[iold];
+ jnum = numneigh[iold];
+
+ // I and J may be ghost atoms
+ // will always know I b/c atoms do not drift that far
+ // but may no longer know J if hops outside cutghost
+ // in that case, assume it performed an event, its strain = qfactor
+
+ i = old2now[iold];
+ emax = maxstrain[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ jold = jlist[jj];
+ j = old2now[jold];
+ if (j >= 0) emax = MAX(emax,maxstrain[j]);
+ else {
+ emax = MAX(emax,qfactor);
+ continue;
+ }
+
+ // DEBUG - only good for 2 dims x,y
+
+ if (j >= nlocal) {
+ if (x[j][0] < sublo[0]) rmaxbig = MAX(rmaxbig,sublo[0]-x[j][0]);
+ if (x[j][1] < sublo[1]) rmaxbig = MAX(rmaxbig,sublo[1]-x[j][1]);
+ if (x[j][0] > subhi[0]) rmaxbig = MAX(rmaxbig,x[j][0]-subhi[0]);
+ if (x[j][1] > subhi[1]) rmaxbig = MAX(rmaxbig,x[j][1]-subhi[1]);
+ if (maxstrain[j] < qfactor) {
+ if (x[j][0] < sublo[0]) rmax = MAX(rmax,sublo[0]-x[j][0]);
+ if (x[j][1] < sublo[1]) rmax = MAX(rmax,sublo[1]-x[j][1]);
+ if (x[j][0] > subhi[0]) rmax = MAX(rmax,x[j][0]-subhi[0]);
+ if (x[j][1] > subhi[1]) rmax = MAX(rmax,x[j][1]-subhi[1]);
+ }
+ }
+ }
+
+ maxstrain_region[i] = emax;
+ }
+
+ double rmax2[2],rmax2all[2];
+ rmax2[0] = rmax;
+ rmax2[1] = rmaxbig;
+ MPI_Allreduce(&rmax2,&rmax2all,2,MPI_DOUBLE,MPI_MAX,world);
+ rmaxever = rmax2all[0];
+ rmaxeverbig = rmax2all[1];
+
+ MPI_Allreduce(&mybonds,&allbonds,1,MPI_INT,MPI_SUM,world);
+
+ //time4 = MPI_Wtime();
+
+ // reverse comm to acquire maxstrain_region from ghost atoms
+ // needed b/c neighbor list referred to old owned atoms,
+ // so above loop may set maxstrain_region of ghost atoms
+ // forward comm to acquire maxstrain_region of all ghost atoms
+
+ commflag = STRAINREGION;
+ comm->reverse_comm_fix(this);
+ comm->forward_comm_fix(this);
+
+ //time5 = MPI_Wtime();
+
+ // identify biased bonds and add to boost list
+ // for each max-strain bond IJ of atom I:
+ // bias this bond only if all these conditions hold:
+ // itag < jtag, so bond is only biased once
+ // maxstrain[i] = maxstrain_region[i]
+ // maxstrain[j] = maxstrain_region[j]
+ // NOTE: also need to check that maxstrain[i] = maxstrain[j] ??
+ // don't think so, b/c maxstrain_region[i] includes maxstrain[i]
+
+ nboost = 0;
+
+ for (i = 0; i < nlocal; i++) {
+ if (numbond[i] == 0) continue;
+ itag = tag[i];
+ j = bonds[i][maxstrain_bondindex[i]].j;
+ jtag = tag[j];
+ if (itag > jtag) continue;
+
+ if (maxstrain[i] != maxstrain_region[i]) continue;
+ if (maxstrain[j] != maxstrain_region[j]) continue;
+
+ if (nboost == maxboost) {
+ maxboost += DELTABOOST;
+ memory->grow(boost,maxboost,"hyper/local:boost");
+ }
+ boost[nboost++] = i;
+ }
+
+ //time6 = MPI_Wtime();
+
+ // apply boost force to bonds with locally max strain
+
+ double **f = atom->f;
+
+ int nobias = 0;
+ mybias = 0.0;
+
+ for (int iboost = 0; iboost < nboost; iboost++) {
+ i = boost[iboost];
+ emax = maxstrain[i];
+ if (emax >= qfactor) {
+ nobias++;
+ continue;
+ }
+
+ m = maxstrain_bondindex[i];
+ j = bonds[i][m].j;
+ r0 = bonds[i][m].r0;
+ boostcoeff = bonds[i][m].boostcoeff;
+
+ vbias = boostcoeff * vmax * (1.0 - emax*emax*invqfactorsq);
+ fbias = boostcoeff * 2.0 * vmax * emax / (qfactor*qfactor * r0);
+
+ delx = x[i][0] - x[j][0];
+ dely = x[i][1] - x[j][1];
+ delz = x[i][2] - x[j][2];
+ r = sqrt(delx*delx + dely*dely + delz*delz);
+ fbiasr = fbias / r;
+
+ f[i][0] += delx*fbiasr;
+ f[i][1] += dely*fbiasr;
+ f[i][2] += delz*fbiasr;
+
+ f[j][0] -= delx*fbiasr;
+ f[j][1] -= dely*fbiasr;
+ f[j][2] -= delz*fbiasr;
+
+ mybias += vbias;
+ }
+
+ //time7 = MPI_Wtime();
+
+ // no boostostat update of boost coeffs when pre_reverse called from setup()
+ // nboost_running, nobias_running, sumboostcoeff only incremented on run steps
+ // NOTE: maybe should also not bias any bonds on firststep of this fix
+
+ if (setupflag) return;
+ nboost_running += nboost;
+ nobias_running += nobias;
+
+ // apply boostostat to boost coefficients of all bonds of all owned atoms
+ // use per-atom full bond list
+ // this is double-calculating for IJ and JI bonds
+ // should be identical for both, b/c emax is the same
+ // could compute once, but would have to find/store index of JI bond
+ // delta in boost coeff is function of maxboost_region vs target boost
+ // maxboost_region is function of two maxstrain_regions for I,J
+ // NOTE: if J is lost to I but not vice versa, then biascoeff IJ != JI
+
+
+
+
+ // DEBUG info for Max and Min bias coeffs every step
+
+ // if (update->ntimestep >= 300000) {
+ if (update->ntimestep >= 0) {
+
+ double emaxi,emaxj,maxboost_region;
+ double mincoeff = 1.0e20;
+ double maxcoeff = -1.0e20;
+ double delta;
+ double maxcoeffprev,deltamax,eimax,ejmax,emaximax,emaxjmax,mbrmax;
+ double mincoeffprev,deltamin,eimin,ejmin,emaximin,emaxjmin,mbrmin;
+ int itagmax,jtagmax,jtagmax2;
+ int itagmin,jtagmin,jtagmin2;
+
+ for (i = 0; i < nlocal; i++) {
+ emaxi = maxstrain_region[i];
+ nbond = numbond[i];
+ for (m = 0; m < nbond; m++) {
+ j = bonds[i][m].j;
+ if (j < 0) continue;
+ emaxj = maxstrain_region[j];
+ emax = MAX(emaxi,emaxj);
+ if (emax < qfactor) vbias = vmax * (1.0 - emax*emax*invqfactorsq);
+ else vbias = 0.0;
+ boostcoeff = bonds[i][m].boostcoeff;
+ maxboost_region = exp(beta * boostcoeff*vbias);
+ delta = -alpha * (maxboost_region-boosttarget) / boosttarget;
+ boostcoeff += delta;
+
+ if (boostcoeff > maxcoeff) {
+ maxcoeffprev = boostcoeff - delta;
+ maxcoeff = boostcoeff;
+ itagmax = tag[i];
+ jtagmax = tag[j];
+ jtagmax2 = bonds[i][m].jtag;
+ deltamax = delta;
+ eimax = maxstrain[i];
+ ejmax = maxstrain[j];
+ emaximax = emaxi;
+ emaxjmax = emaxj;
+ mbrmax = maxboost_region;
+ }
+
+ if (boostcoeff < mincoeff) {
+ mincoeffprev = boostcoeff - delta;
+ mincoeff = boostcoeff;
+ itagmin = tag[i];
+ jtagmin = tag[j];
+ jtagmin2 = bonds[i][m].jtag;
+ deltamin = delta;
+ eimin = maxstrain[i];
+ ejmin = maxstrain[j];
+ emaximin = emaxi;
+ emaxjmin = emaxj;
+ mbrmin = maxboost_region;
+ }
+ }
+ }
+
+ // min/max owner = procs that own the Min and Max bias coeffs
+
+ pairme.value = mincoeff;
+ pairme.proc = me;
+ MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MINLOC,world);
+ int minowner = pairall.proc;
+
+ pairme.value = maxcoeff;
+ pairme.proc = me;
+ MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MAXLOC,world);
+ int maxowner = pairall.proc;
+
+ // get info about which atom and bond influenced Max coeff
+ // querystrain = strain of bond that is influencing the max Cij
+ // inbondi,inbondj = ID of I,J atoms in the influence bond
+ // incoeff = bias coeff of influence bond
+ // inbias = 1/0 if inflence bond is biased or not on this step
+ // instraddle = 1/0 if influence bond straddles proc boundary or not
+
+ double querystrain;
+ if (maxowner == me) querystrain = MAX(emaximax,emaxjmax);
+ MPI_Bcast(&querystrain,1,MPI_DOUBLE,maxowner,world);
+
+ int count = 0;
+ for (i = 0; i < nlocal; i++) {
+ m = maxstrain_bondindex[i];
+ if (tag[i] > bonds[i][m].jtag) continue;
+ if (querystrain == maxstrain[i]) count++;
+ }
+
+ int countall;
+ MPI_Allreduce(&count,&countall,1,MPI_INT,MPI_SUM,world);
+
+ int inbias = 0;
+ int instraddle = countall; // output in rare case countall != 1
+ tagint inbondi = 0;
+ tagint inbondj = 0;
+ double incoeff = 0.0;
+
+ if (countall == 1) {
+ int procowner;
+ int proc = 0;
+ for (i = 0; i < nlocal; i++) {
+ m = maxstrain_bondindex[i];
+ if (tag[i] > bonds[i][m].jtag) continue;
+ if (querystrain == maxstrain[i]) {
+ proc = me;
+ inbondi = tag[i];
+ inbondj = bonds[i][m].jtag;
+ incoeff = bonds[i][m].boostcoeff;
+ instraddle = 0;
+ if (bonds[i][m].j >= nlocal) instraddle = 1;
+ inbias = 0;
+ for (int iboost = 0; iboost < nboost; iboost++)
+ if (boost[iboost] == i) inbias = 1;
+ }
+ }
+
+ MPI_Allreduce(&proc,&procowner,1,MPI_INT,MPI_MAX,world);
+ MPI_Bcast(&inbondi,1,MPI_LMP_TAGINT,procowner,world);
+ MPI_Bcast(&inbondj,1,MPI_LMP_TAGINT,procowner,world);
+ MPI_Bcast(&incoeff,1,MPI_DOUBLE,procowner,world);
+ MPI_Bcast(&instraddle,1,MPI_INT,procowner,world);
+ MPI_Bcast(&inbias,1,MPI_INT,procowner,world);
+ }
+
+ // output by procs that own the Min and Max bias coeffs
+
+ if (minowner == me) {
+ printf("MinCoeff %ld %d : %d %d %d : prev/delta/new %g %g %g : "
+ "Eij %g %g %g %g : BR %g\n",
+ update->ntimestep,comm->me,itagmin,jtagmin,jtagmin2,
+ mincoeffprev,deltamin,mincoeff,eimin,ejmin,emaximin,emaxjmin,mbrmin);
+ }
+ if (maxowner == me) {
+ printf("MaxCoeff %ld %d : %d %d %d : prev/delta/new %g %g %g : "
+ "Eij %g %g %g %g : Bond %d %d %g %d %d : BR %g\n",
+ update->ntimestep,comm->me,itagmax,jtagmax,jtagmax2,
+ maxcoeffprev,deltamax,maxcoeff,
+ eimax,ejmax,emaximax,emaxjmax,
+ inbondi,inbondj,incoeff,inbias,instraddle,mbrmax);
+ }
+ }
+
+ // end of DEBUG
+
+
+
+
+
+ double myboost = 0.0;
+ double emaxi,emaxj,maxboost_region;
+
+ for (i = 0; i < nlocal; i++) {
+ emaxi = maxstrain_region[i];
+ nbond = numbond[i];
+ for (m = 0; m < nbond; m++) {
+ j = bonds[i][m].j;
+ if (j < 0) continue;
+ emaxj = maxstrain_region[j];
+ emax = MAX(emaxi,emaxj);
+ if (emax < qfactor) vbias = vmax * (1.0 - emax*emax*invqfactorsq);
+ else vbias = 0.0;
+ boostcoeff = bonds[i][m].boostcoeff;
+ maxboost_region = exp(beta * boostcoeff*vbias);
+ boostcoeff -= alpha * (maxboost_region-boosttarget) / boosttarget;
+ // COMMENT OUT for now
+ //boostcoeff = MIN(boostcoeff,COEFFMAX);
+ myboost += boostcoeff;
+ maxboostcoeff = MAX(maxboostcoeff,boostcoeff);
+ minboostcoeff = MIN(minboostcoeff,boostcoeff);
+ bonds[i][m].boostcoeff = boostcoeff;
+ }
+ }
+
+ // running stats
+
+ MPI_Allreduce(&myboost,&allboost,1,MPI_DOUBLE,MPI_SUM,world);
+ if (allbonds) sumboostcoeff += allboost/allbonds;
+
+ // histogram the bond coeffs and output as requested
+ // do not double count each bond
+
+ if (histoflag && update->ntimestep % histo_every == 0) {
+ if (histo_steps == 0)
+ for (i = 0; i < histo_count+2; i++) histo[i] = 0;
+ histo_steps++;
+
+ int ihisto;
+ for (i = 0; i < nlocal; i++) {
+ nbond = numbond[i];
+ for (m = 0; m < nbond; m++) {
+ if (tag[i] > bonds[i][m].jtag) continue;
+ boostcoeff = bonds[i][m].boostcoeff;
+ if (boostcoeff < histo_lo) ihisto = -1;
+ else ihisto = static_cast<int> ((boostcoeff-histo_lo) * invhisto_delta);
+ if (ihisto >= histo_count) ihisto = histo_count;
+ histo[ihisto+1]++;
+ }
+ }
+
+ if (update->ntimestep % histo_print == 0) {
+ MPI_Allreduce(histo,allhisto,histo_count+2,MPI_LMP_BIGINT,MPI_SUM,world);
+
+ bigint total;
+ for (i = 0; i < histo_count+2; i++) total += allhisto[i];
+
+ if (me == 0) {
+ if (screen) {
+ fprintf(screen,"Histogram of bias coeffs:\n");
+ for (i = 0; i < histo_count+2; i++)
+ fprintf(screen," %g",1.0*allhisto[i]/total);
+ fprintf(screen,"\n");
+ }
+ if (logfile) {
+ fprintf(logfile,"Histogram of bias coeffs:\n");
+ for (i = 0; i < histo_count+2; i++)
+ fprintf(logfile," %g",1.0*allhisto[i]/total);
+ fprintf(logfile,"\n");
+ }
+ }
+ }
+ }
+
+ // check for any biased bonds that are too close to each other
+ // keep a running count for output
+
+ if (checkbias && update->ntimestep % checkbias_every == 0) {
+
+ // mark each atom in a biased bond with ID of partner
+ // nboost loop will mark some ghost atoms
+
+ for (i = 0; i < nall; i++) biasflag[i] = 0;
+
+ for (int iboost = 0; iboost < nboost; iboost++) {
+ i = boost[iboost];
+ m = maxstrain_bondindex[i];
+ j = bonds[i][m].j;
+ biasflag[i] = tag[j];
+ biasflag[j] = tag[i];
+ }
+
+ // reverse comm to acquire biasflag from ghost atoms
+ // needed b/c above loop may set biasflag of ghost atoms
+ // forward comm to acquire biasflag of all ghost atoms
+
+ commflag = BIASFLAG;
+ comm->reverse_comm_fix(this);
+ comm->forward_comm_fix(this);
+
+ // I and J may be ghost atoms
+ // only continue if I is a biased atom
+ // if J is unknonw (drifted ghost) just ignore
+ // if J is biased and is not bonded to I, then flag as too close
+
+ for (ii = 0; ii < inum; ii++) {
+ iold = ilist[ii];
+ i = old2now[iold];
+ if (biasflag[i] == 0) continue;
+
+ jlist = firstneigh[iold];
+ jnum = numneigh[iold];
+
+ for (jj = 0; jj < jnum; jj++) {
+ jold = jlist[jj];
+ j = old2now[jold];
+ if (j < 0) continue;
+ if (biasflag[j] && biasflag[j] != tag[i]) checkbias_count++;
+ }
+ }
+ }
+
+ // check for any bond bias coeffcients that do not match
+ // cannot check unless both atoms IJ are owned by this proc
+ // keep a running count for output
+
+ if (checkcoeff && update->ntimestep % checkcoeff_every == 0) {
+ int jb,jbonds;
+ double ibias,jbias;
+
+ for (i = 0; i < nlocal; i++) {
+ nbond = numbond[i];
+ for (m = 0; m < nbond; m++) {
+ if (tag[i] > bonds[i][m].jtag) continue;
+ j = bonds[i][m].j;
+ if (j < 0) continue;
+ if (j >= nlocal) continue;
+ itag = tag[i];
+ jbonds = numbond[j];
+ for (jb = 0; jb < jbonds; jb++)
+ if (bonds[j][jb].jtag == itag) break;
+ if (jb == jbonds)
+ error->one(FLERR,"Fix hyper/local could not find duplicate bond");
+ if (bonds[i][m].boostcoeff != bonds[j][jb].boostcoeff)
+ checkcoeff_count++;
+ }
+ }
+ }
+
+ // DEBUG timing output
+ // examine these timings for new hyper/local before delete this code
+ // verify that boost cooefs stay in sync for 2 copies of each bond
+ // before delete this code
+
+ /*
+ time8 = MPI_Wtime();
+
+ timefirst += time2-time1;
+ timesecond += time3-time2;
+ timethird += time4-time3;
+ timefourth += time5-time4;
+ timefifth += time6-time5;
+ timesixth += time7-time6;
+ timeseventh += time8-time7;
+ timetotal += time8-time1;
+
+ double bpera = compute_vector(6);
+ double nbperb = compute_vector(7);
+
+ if (update->ntimestep == 10000 && me == 0) {
+ printf("TIME First %g\n",timefirst);
+ printf("TIME Second %g\n",timesecond);
+ printf("TIME Third %g\n",timethird);
+ printf("TIME Fourth %g\n",timefourth);
+ printf("TIME Fifth %g\n",timefifth);
+ printf("TIME Sixth %g\n",timesixth);
+ printf("TIME Seventh %g\n",timeseventh);
+ printf("TIME Total %g\n",timetotal);
+ printf("Bonds/atom, neigh-bonds/bond %g %g\n",bpera,nbperb);
+ also print nbondbuild and tbondbuild
+ }
+ */
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperLocal::min_pre_neighbor()
+{
+ pre_neighbor();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperLocal::build_bond_list(int natom)
+{
+ int i,j,ii,jj,m,n,inum,jnum,nbond;
+ tagint itag,jtag;
+ double xtmp,ytmp,ztmp,delx,dely,delz,rsq,oldcoeff;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ double time1,time2;
+ time1 = MPI_Wtime();
+
+ if (natom) {
+ nevent++;
+ nevent_atom += natom;
+ }
+
+ // trigger Dcut neighbor list build
+ // NOTE: turn off special bonds in this Dcut neigh list?
+
+ neighbor->build_one(list);
+
+ // make copy of old bonds to preserve boostcoeffs for bonds that persist
+ // allocate new numbond
+
+ OneBond **old_bonds = bonds;
+ int *old_numbond = numbond;
+
+ int nmax = atom->nmax;
+ memory->create(numbond,nmax,"hyper/local:numbond");
+
+ // old_nall = value of nall at time bonds are built
+ // reallocate new xold and tagold if necessary
+ // initialize xold to current coords
+ // initialize tagold to zero, so atoms not in neighbor list will remain zero
+
+ old_nall = atom->nlocal + atom->nghost;
+
+ if (old_nall > maxold) {
+ memory->destroy(xold);
+ memory->destroy(tagold);
+ memory->destroy(old2now);
+ maxold = atom->nmax;
+ memory->create(xold,maxold,3,"hyper/local:xold");
+ memory->create(tagold,maxold,"hyper/local:tagold");
+ memory->create(old2now,maxold,"hyper/local:old2now");
+ }
+
+ double **x = atom->x;
+
+ memcpy(&xold[0][0],&x[0][0],3*old_nall*sizeof(double));
+ for (i = 0; i < old_nall; i++) tagold[i] = 0;
+
+ // create and populate new bonds data struct
+ // while loop allows maxbondperatom to increase once if necessary
+ // don't know new maxbondperatom value until end of loop
+ // in practice maxbondperatom will hardly ever increase
+ // since there is a physical max value
+
+ tagint *tag = atom->tag;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ while (1) {
+ bonds = (OneBond **) memory->create(bonds,nmax,maxbondperatom,
+ "hyper/local:bonds");
+ if (bonds) memset(bonds[0],0,nmax*sizeof(OneBond));
+ for (i = 0; i < nlocal; i++) numbond[i] = 0;
+
+ // identify bonds assigned to each owned atom
+ // do not create a bond between two non-group atoms
+ // set tagold = global ID for all I,J atoms used in neighbor list
+ // tagold remains 0 for unused atoms, skipped in pre_neighbor
+
+ int nbondmax = 0;
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ itag = tag[i];
+ tagold[i] = tag[i];
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+ nbond = 0;
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+ jtag = tag[j];
+ tagold[j] = jtag;
+
+ // skip if neither atom I or J are in fix group
+ // order IJ to insure IJ and JI bonds are stored consistently
+
+ if (!(mask[i] & groupbit) && !(mask[j] & groupbit)) continue;
+
+ if (itag < jtag) {
+ delx = xtmp - x[j][0];
+ dely = ytmp - x[j][1];
+ delz = ztmp - x[j][2];
+ } else {
+ delx = x[j][0] - xtmp;
+ dely = x[j][1] - ytmp;
+ delz = x[j][2] - ztmp;
+ }
+
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ // NOTE: could create two bonds for IJ both owned from one calc?
+ // have to skip one of 2 bonds in that case
+
+ if (rsq < cutbondsq) {
+ if (nbond >= maxbondperatom) {
+ nbond++;
+ continue;
+ }
+
+ bonds[i][nbond].r0 = sqrt(rsq);
+ bonds[i][nbond].jtag = tag[j];
+ bonds[i][nbond].j = j;
+
+ if (firstflag) oldcoeff = 0.0;
+ else {
+ oldcoeff = 0.0;
+ jtag = tag[j];
+ n = old_numbond[i];
+ for (m = 0; m < n; m++) {
+ if (old_bonds[i][m].jtag == jtag) {
+ oldcoeff = old_bonds[i][m].boostcoeff;
+ break;
+ }
+ }
+ }
+
+ if (oldcoeff > 0.0) bonds[i][nbond].boostcoeff = oldcoeff;
+ else {
+ bonds[i][nbond].boostcoeff = BOOSTINIT;
+ nnewbond++;
+ }
+ nbond++;
+ }
+ }
+ numbond[i] = nbond;
+ nbondmax = MAX(nbondmax,nbond);
+ }
+
+ // maxbondperatom must increase uniformly on all procs
+ // since bonds are comunicated when atoms migrate
+
+ int allnbondmax;
+ MPI_Allreduce(&nbondmax,&allnbondmax,1,MPI_INT,MPI_MAX,world);
+ if (allnbondmax <= maxbondperatom) break;
+
+ maxbondperatom = allnbondmax;
+ memory->destroy(bonds);
+ }
+
+ // deallocate old_bonds and old_numbond
+
+ memory->destroy(old_bonds);
+ memory->destroy(old_numbond);
+
+ time2 = MPI_Wtime();
+ if (firstflag) nnewbond = 0;
+ else {
+ time_bondbuild += time2-time1;
+ nbondbuild++;
+ }
+ firstflag = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixHyperLocal::pack_forward_comm(int n, int *list, double *buf,
+ int pbc_flag, int *pbc)
+{
+ int i,j,k,m,ibond,nbond,start;
+
+ m = 0;
+
+ // STRAIN
+ // pack maxstrain vector
+
+ if (commflag == STRAIN) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = maxstrain[j];
+ }
+
+ // STRAINREGION
+ // pack maxstrain_region vector
+
+ } else if (commflag == STRAINREGION) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = maxstrain_region[j];
+ }
+
+ // BIASFLAG
+ // pack biasflag vector
+
+ } else if (commflag == BIASFLAG) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ buf[m++] = ubuf(biasflag[j]).d;
+ }
+ }
+
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperLocal::unpack_forward_comm(int n, int first, double *buf)
+{
+ int i,j,k,m,last,ibond,nbond,jlocal,flag;
+
+ m = 0;
+ last = first + n;
+
+ // STRAIN
+ // unpack maxstrain vector
+
+ if (commflag == STRAIN) {
+ for (i = first; i < last; i++) {
+ maxstrain[i] = buf[m++];
+ }
+
+ // STRAINREGION
+ // unpack maxstrain_region vector
+
+ } else if (commflag == STRAINREGION) {
+ for (i = first; i < last; i++) {
+ maxstrain_region[i] = buf[m++];
+ }
+
+ // BIASFLAG
+ // unpack biasflag vector
+
+ } else if (commflag == BIASFLAG) {
+ for (i = first; i < last; i++) {
+ biasflag[i] = (tagint) ubuf(buf[m++]).i;
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixHyperLocal::pack_reverse_comm(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+
+ // STRAINREGION
+ // pack maxstrain_region vector
+
+ if (commflag == STRAINREGION) {
+ for (i = first; i < last; i++) {
+ buf[m++] = maxstrain_region[i];
+ }
+
+ // BIASFLAG
+ // pack biasflag vector
+
+ } else if (commflag == BIASFLAG) {
+ for (i = first; i < last; i++) {
+ buf[m++] = ubuf(biasflag[i]).d;
+ }
+ }
+
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixHyperLocal::unpack_reverse_comm(int n, int *list, double *buf)
+{
+ int i,j,m;
+
+ m = 0;
+
+ // STRAINREGION
+ // unpack maxstrain_region vector
+
+ if (commflag == STRAINREGION) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ maxstrain_region[j] += buf[m++];
+ }
+
+ // BIASFLAG
+ // unpack biasflag vector
+
+ } else if (commflag == BIASFLAG) {
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ biasflag[j] = (tagint) ubuf(buf[m++]).i;
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ allocate atom-based arrays
+------------------------------------------------------------------------- */
+
+void FixHyperLocal::grow_arrays(int nmax)
+{
+ // NOTE: not all of these need to be Nmax in length, could allocate elsewhere
+
+ memory->grow(maxstrain,nmax,"hyper/local:maxstrain");
+ memory->grow(maxstrain_bondindex,nmax,"hyper/local:maxstrain_bondindex");
+ memory->grow(maxstrain_region,nmax,"hyper/local:maxstrain_region");
+ if (checkbias) memory->grow(biasflag,nmax,"hyper/local:biasflag");
+
+ memory->grow(numbond,nmax,"hyper/local:numbond");
+ memory->grow(bonds,nmax,maxbondperatom,"hyper/local:bonds");
+
+ // zero so valgrind does not complain about memcpy() in copy()
+ // also so loops in pre_neighbor() are OK before
+ // bonds are setup for the first time
+
+ if (bonds) {
+ memset(bonds[maxbond],0,(nmax-maxbond)*maxbondperatom*sizeof(OneBond));
+ memset(&numbond[maxbond],0,(nmax-maxbond)*sizeof(int));
+ maxbond = nmax;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ copy values within local atom-based arrays
+------------------------------------------------------------------------- */
+
+void FixHyperLocal::copy_arrays(int i, int j, int delflag)
+{
+ // avoid valgrind copy-to-self warning
+
+ if (i != j) memcpy(bonds[j],bonds[i],numbond[i]*sizeof(OneBond));
+ numbond[j] = numbond[i];
+}
+
+/* ----------------------------------------------------------------------
+ pack values in local atom-based array for exchange with another proc
+------------------------------------------------------------------------- */
+
+int FixHyperLocal::pack_exchange(int i, double *buf)
+{
+ int m = 1;
+ int n = numbond[i];
+ buf[m++] = ubuf(n).d;
+ for (int j = 0; j < n; j++) {
+ buf[m++] = bonds[i][j].r0;
+ buf[m++] = bonds[i][j].boostcoeff;
+ buf[m++] = ubuf(bonds[i][j].jtag).d;
+ buf[m++] = ubuf(bonds[i][j].j).d;
+ }
+
+ buf[0] = m;
+ return m;
+}
+
+/* ----------------------------------------------------------------------
+ unpack values in local atom-based array from exchange with another proc
+------------------------------------------------------------------------- */
+
+int FixHyperLocal::unpack_exchange(int nlocal, double *buf)
+{
+ int m = 1;
+ int n = numbond[nlocal] = (int) ubuf(buf[m++]).i;
+ for (int j = 0; j < n; j++) {
+ bonds[nlocal][j].r0 = buf[m++];
+ bonds[nlocal][j].boostcoeff = buf[m++];
+ bonds[nlocal][j].jtag = (tagint) ubuf(buf[m++]).i;
+ bonds[nlocal][j].j = (int) ubuf(buf[m++]).i;
+ }
+
+ return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixHyperLocal::compute_scalar()
+{
+ double allbias;
+ MPI_Allreduce(&mybias,&allbias,1,MPI_DOUBLE,MPI_SUM,world);
+ return allbias;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixHyperLocal::compute_vector(int i)
+{
+ // 23 vector outputs returned for i = 0-22
+
+ // i = 0 = # of boosted bonds on this step
+ // i = 1 = max strain of any bond on this step
+ // i = 2 = average bias potential for all bonds on this step
+ // i = 3 = ave bonds/atom on this step
+ // i = 4 = ave neighbor bonds/bond on this step
+
+ // i = 5 = fraction of steps and bonds with no bias during this run
+ // i = 6 = max drift distance of any atom during this run
+ // i = 7 = max bond length during this run
+ // i = 8 = average # of boosted bonds/step during this run
+ // i = 9 = average bias potential for all bonds during this run
+ // i = 10 = max bias potential for any bond during this run
+ // i = 11 = min bias potential for any bond during this run
+ // i = 12 = max dist from my box of any ghost atom with
+ // maxstain < qfactor during this run
+ // i = 13 = max dist from my box of any ghost atom with
+ // any maxstrain during this run
+ // i = 14 = count of ghost atoms that could not be found
+ // by any proc at any reneighbor step during this run
+ // i = 15 = count of lost bond partners during this run
+ // i = 16 = average bias coeff for lost bond partners during this run
+ // i = 17 = count of bias overlaps found during this run
+ // i = 18 = count of non-matching bias coefficients found during this run
+
+ // i = 19 = cummulative hyper time
+ // i = 20 = cummulative # of event timesteps since fix created
+ // i = 21 = cummulative # of atoms in events since fix created
+ // i = 22 = cummulative # of new bonds formed since fix created
+
+ if (i == 0) {
+ int nboostall;
+ MPI_Allreduce(&nboost,&nboostall,1,MPI_INT,MPI_SUM,world);
+ return (double) nboostall;
+ }
+
+ if (i == 1) {
+ if (nostrainyet) return 0.0;
+ int nlocal = atom->nlocal;
+ double emax = 0.0;
+ for (int j = 0; j < nlocal; j++)
+ emax = MAX(emax,maxstrain[j]);
+ double eall;
+ MPI_Allreduce(&emax,&eall,1,MPI_DOUBLE,MPI_MAX,world);
+ return eall;
+ }
+
+ if (i == 2) {
+ if (allboost && allbonds) return allboost/allbonds * vmax;
+ return 1.0;
+ }
+
+ if (i == 3) return 1.0*allbonds/atom->natoms;
+
+ if (i == 4) {
+ int nlocal = atom->nlocal;
+ int nbonds = 0; // BIGINT?
+ for (int j = 0; j < nlocal; j++)
+ nbonds += numbond[j];
+ int allbonds;
+ MPI_Allreduce(&nbonds,&allbonds,1,MPI_INT,MPI_SUM,world);
+ int allneigh; // BIGINT?
+ MPI_Allreduce(&list->ipage->ndatum,&allneigh,1,MPI_INT,MPI_SUM,world);
+ double neighsperatom = allneigh/atom->natoms;
+ double bondsperatom = 0.5*allbonds/atom->natoms;
+ return neighsperatom * bondsperatom;
+ }
+
+ if (i == 5) {
+ int allboost_running,allnobias_running;
+ MPI_Allreduce(&nboost_running,&allboost_running,1,MPI_INT,MPI_SUM,world);
+ MPI_Allreduce(&nobias_running,&allnobias_running,1,MPI_INT,MPI_SUM,world);
+ if (allboost_running) return 1.0*allnobias_running / allboost_running;
+ return 0.0;
+ }
+
+ if (i == 6) {
+ double alldriftsq;
+ MPI_Allreduce(&maxdriftsq,&alldriftsq,1,MPI_DOUBLE,MPI_MAX,world);
+ return (double) sqrt(alldriftsq);
+ }
+
+ if (i == 7) {
+ double allbondlen;
+ MPI_Allreduce(&maxbondlen,&allbondlen,1,MPI_DOUBLE,MPI_MAX,world);
+ return allbondlen;
+ }
+
+ if (i == 8) {
+ if (update->ntimestep == update->firststep) return 0.0;
+ int allboost_running;
+ MPI_Allreduce(&nboost_running,&allboost_running,1,MPI_INT,MPI_SUM,world);
+ return 1.0*allboost_running / (update->ntimestep - update->firststep);
+ }
+
+ if (i == 9) {
+ if (update->ntimestep == update->firststep) return 0.0;
+ return sumboostcoeff * vmax / (update->ntimestep - update->firststep);
+ }
+
+ if (i == 10) {
+ double allboostcoeff;
+ MPI_Allreduce(&maxboostcoeff,&allboostcoeff,1,MPI_DOUBLE,MPI_MAX,world);
+ return allboostcoeff * vmax;
+ }
+
+ if (i == 11) {
+ double allboostcoeff;
+ MPI_Allreduce(&minboostcoeff,&allboostcoeff,1,MPI_DOUBLE,MPI_MAX,world);
+ return allboostcoeff * vmax;
+ }
+
+ if (i == 12) return rmaxever;
+ if (i == 13) return rmaxeverbig;
+
+ if (i == 14) {
+ int allghost_toofar;
+ MPI_Allreduce(&ghost_toofar,&allghost_toofar,1,MPI_INT,MPI_SUM,world);
+ return 1.0*allghost_toofar;
+ }
+
+ if (i == 15) {
+ int alllost;
+ MPI_Allreduce(&lostbond_partner,&alllost,1,MPI_INT,MPI_SUM,world);
+ return 1.0*alllost;
+ }
+
+ if (i == 16) {
+ int alllost;
+ MPI_Allreduce(&lostbond_partner,&alllost,1,MPI_INT,MPI_SUM,world);
+ double allcoeff;
+ MPI_Allreduce(&lostbond_coeff,&allcoeff,1,MPI_DOUBLE,MPI_SUM,world);
+ if (alllost == 0) return 0;
+ return allcoeff/alllost;
+ }
+
+ if (i == 17) {
+ int allclose;
+ MPI_Allreduce(&checkbias_count,&allclose,1,MPI_INT,MPI_SUM,world);
+ return 1.0*allclose;
+ }
+
+ if (i == 18) {
+ int allcoeff;
+ MPI_Allreduce(&checkcoeff_count,&allcoeff,1,MPI_INT,MPI_SUM,world);
+ return 1.0*allcoeff;
+ }
+
+ if (i == 19) {
+ return boosttarget * update->dt * (update->ntimestep - starttime);
+ }
+
+ if (i == 20) return (double) nevent;
+ if (i == 21) return (double) nevent_atom;
+
+ if (i == 22) {
+ int allnew;
+ MPI_Allreduce(&nnewbond,&allnew,1,MPI_INT,MPI_SUM,world);
+ return (double) 0.5*allnew;
+ }
+}
+
+/* ----------------------------------------------------------------------
+ wrapper on compute_vector()
+ used by hyper.cpp to call FixHyper
+------------------------------------------------------------------------- */
+
+double FixHyperLocal::query(int i)
+{
+ if (i == 1) return compute_vector(19); // cummulative hyper time
+ if (i == 2) return compute_vector(20); // nevent
+ if (i == 3) return compute_vector(21); // nevent_atom
+ if (i == 4) return compute_vector(3); // ave bonds/atom
+ if (i == 5) return compute_vector(6); // maxdrift
+ if (i == 6) return compute_vector(7); // maxbondlen
+ if (i == 7) return compute_vector(5); // fraction with zero bias
+
+ // unique to local hyper
+
+ if (i == 8) return compute_vector(22); // number of new bonds
+ if (i == 9) return 1.0*maxbondperatom; // max bonds/atom
+ if (i == 10) return compute_vector(8); // ave # of boosted bonds/step
+ if (i == 11) return compute_vector(9); // ave boost coeff over all bonds
+ if (i == 12) return compute_vector(10); // max boost cooef for any bond
+ if (i == 13) return compute_vector(11); // max boost cooef for any bond
+ if (i == 14) return compute_vector(4); // neighbor bonds/bond
+ if (i == 15) return compute_vector(2); // ave boost cooef now
+ if (i == 16) return time_bondbuild; // CPU time for bond_build calls
+ if (i == 17) return rmaxever; // ghost atom distance for < maxstrain
+ if (i == 18) return rmaxeverbig; // ghost atom distance for any strain
+ if (i == 19) return compute_vector(14); // count of ghost atoms not found
+ if (i == 20) return compute_vector(15); // count of lost bond partners
+ if (i == 21) return compute_vector(16); // ave bias coeff of long bonds
+ if (i == 22) return compute_vector(17); // count of bias overlaps
+ if (i == 23) return compute_vector(18); // count of non-matching bias coeffs
+
+ error->all(FLERR,"Invalid query to fix hyper/local");
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of per-atom and per-bond arrays
+------------------------------------------------------------------------- */
+
+double FixHyperLocal::memory_usage()
+{
+ int nmax = atom->nmax;
+ double bytes = 2*nmax * sizeof(int); // numbond, maxstrain_bondindex
+ bytes = 2*nmax * sizeof(double); // maxstrain, maxstrain_region
+ bytes += maxbondperatom*nmax * sizeof(OneBond); // bonds
+ bytes += 3*maxold * sizeof(double); // xold
+ bytes += maxold * sizeof(tagint); // tagold
+ bytes += maxold * sizeof(int); // old2now
+ return bytes;
+}
diff --git a/src/REPLICA/fix_hyper_local.h b/src/REPLICA/fix_hyper_local.h
new file mode 100644
index 0000000000..5150386929
--- /dev/null
+++ b/src/REPLICA/fix_hyper_local.h
@@ -0,0 +1,175 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(hyper/local,FixHyperLocal)
+
+#else
+
+#ifndef LMP_FIX_HYPER_LOCAL_H
+#define LMP_FIX_HYPER_LOCAL_H
+
+#include "fix_hyper.h"
+
+namespace LAMMPS_NS {
+
+class FixHyperLocal : public FixHyper {
+ public:
+ FixHyperLocal(class LAMMPS *, int, char **);
+ ~FixHyperLocal();
+ int setmask();
+ void init();
+ void init_list(int, class NeighList *);
+ void setup_pre_neighbor();
+ void setup_pre_reverse(int, int);
+ void pre_neighbor();
+ void pre_reverse(int, int);
+ void min_pre_neighbor();
+ double compute_scalar();
+ double compute_vector(int);
+ double query(int);
+
+ int pack_forward_comm(int, int *, double *, int, int *);
+ void unpack_forward_comm(int, int, double *);
+ int pack_reverse_comm(int, int, double *);
+ void unpack_reverse_comm(int, int *, double *);
+
+ void grow_arrays(int);
+ void copy_arrays(int, int, int);
+ int pack_exchange(int, double *);
+ int unpack_exchange(int, double *);
+
+ double memory_usage();
+
+ // extra methods visible to callers
+
+ void init_hyper();
+ void build_bond_list(int);
+
+ private:
+ int me;
+ double cutbond,qfactor,vmax,tequil,dcut;
+ double alpha_user; // timescale to apply boostostat (time units)
+ double alpha; // unitless dt/alpha_user
+ double boosttarget; // target value of boost
+ int histoflag;
+ int lostbond,lostbond_partner;
+ double lostbond_coeff;
+ int checkbias,checkbias_every,checkbias_flag,checkbias_count;
+ int checkcoeff,checkcoeff_every,checkcoeff_flag,checkcoeff_count;
+
+ int setupflag; // 1 during setup, 0 during run
+ int firstflag; // set for first time bond_build takes place
+ int nostrainyet; // 1 until maxstrain is first computed
+
+ int nboost_running,nobias_running;
+ int nbondbuild;
+ double time_bondbuild;
+ bigint starttime;
+ double sumboostcoeff; // sum of aveboost at every timestep
+ int allbonds; // sum of bond count on this step
+ double allboost; // sum of boostcoeff on all bonds on this step
+
+ int nnewbond; // running tally of number of new bonds created
+ int maxbondperatom; // max # of bonds any atom ever has
+ int commflag; // flag for communication mode
+ int nevent; // # of events that trigger bond rebuild
+ int nevent_atom; // # of atoms that experienced an event
+ double cutbondsq,dcutsq;
+ double beta,t_hyper,invqfactorsq;
+ double mybias;
+ double maxbondlen; // cummulative max length of any bond
+ double maxdriftsq; // max distance any atom drifts from original pos
+ double maxboostcoeff; // cummulative max boost coeff for any bond
+ double minboostcoeff; // cummulative min boost coeff for any bond
+ double rmaxever,rmaxeverbig;
+ int ghost_toofar;
+
+ double timefirst,timesecond,timethird,timefourth;
+ double timefifth,timesixth,timeseventh,timetotal;
+
+ // data structs for per-atom and per-bond info
+ // all of these are for current owned and ghost atoms
+ // except list and old2now are for atom indices at time of last bond build
+
+ class NeighList *list; // full neigh list up to Dcut distance
+ // created only when bonds are rebuilt
+
+ int *old2now; // o2n[i] = current local index of old atom i
+ // stored for old owned and ghost atoms
+ // I = old index when bonds were last created
+ // old indices are stored in old neighbor list
+
+ double **xold; // coords of owned+ghost atoms when bonds created
+ tagint *tagold; // global IDs of owned+ghost atoms when b created
+
+ int maxold; // allocated size of old2now
+ int maxbond; // allocated size of bonds
+ int old_nall; // nlocal+nghost when old2now was last setup
+
+ struct OneBond { // single IJ bond, atom I is owner
+ double r0; // original relaxed bond length
+ double boostcoeff; // boost coefficient
+ tagint jtag; // global index of J atom in bond IJ
+ int j; // local index of J atom in bond IJ
+ };
+
+ struct OneBond **bonds; // 2d array of bonds for owned atoms
+ int *numbond; // number of bonds for each owned atom
+
+ double *maxstrain; // max-strain of any bond atom I is part of
+ // for owned and ghost atoms
+ double *maxstrain_region; // max-strain of any neighbor atom J of atom I
+ // for owned and ghost atoms
+ int *maxstrain_bondindex; // index of max-strain bond of each atom I
+ // just for owned atoms
+ tagint *biasflag; // atoms in biased bonds marked with bond partner
+ // for owned and ghost atoms
+
+ // list of boosted bonds that this proc will bias
+
+ int maxboost; // allocated size of boost list
+ int nboost; // # of boosted bonds I own
+ int *boost; // index of atom I in each boosted bond
+
+ // histogramming of bond boost cooeficients
+
+ int histo_flag,histo_every,histo_count,histo_print,histo_steps;
+ double histo_delta,invhisto_delta,histo_lo;
+ bigint *histo,*allhisto;
+
+ // DEBUG: MPI data struct for finding min/max bias coeffs via Allreduce
+
+ struct Two {
+ double value;
+ int proc;
+ };
+ Two pairme,pairall;
+
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+*/
diff --git a/src/REPLICA/hyper.cpp b/src/REPLICA/hyper.cpp
new file mode 100644
index 0000000000..0449943215
--- /dev/null
+++ b/src/REPLICA/hyper.cpp
@@ -0,0 +1,548 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include "hyper.h"
+#include "update.h"
+#include "atom.h"
+#include "domain.h"
+#include "region.h"
+#include "integrate.h"
+#include "min.h"
+#include "force.h"
+#include "neighbor.h"
+#include "modify.h"
+#include "compute_event_displace.h"
+#include "fix_hyper.h"
+#include "fix_event_hyper.h"
+#include "output.h"
+#include "dump.h"
+#include "finish.h"
+#include "timer.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{NOHYPER,GLOBAL,LOCAL};
+
+/* ---------------------------------------------------------------------- */
+
+Hyper::Hyper(LAMMPS *lmp) : Pointers(lmp) {}
+
+/* ----------------------------------------------------------------------
+ perform hyperdynamics simulation
+------------------------------------------------------------------------- */
+
+void Hyper::command(int narg, char **arg)
+{
+ MPI_Comm_rank(world,&me);
+ MPI_Comm_size(world,&nprocs);
+
+ // error checks
+
+ if (domain->box_exist == 0)
+ error->all(FLERR,"Hyper command before simulation box is defined");
+
+ if (narg < 4) error->all(FLERR,"Illegal hyper command");
+
+ int nsteps = force->inumeric(FLERR,arg[0]);
+ t_event = force->inumeric(FLERR,arg[1]);
+
+ char *id_fix = new char[strlen(arg[2])+1];
+ strcpy(id_fix,arg[2]);
+
+ char *id_compute = new char[strlen(arg[3])+1];
+ strcpy(id_compute,arg[3]);
+
+ options(narg-4,&arg[4]);
+
+ // total # of timesteps must be multiple of t_event
+
+ if (t_event <= 0)
+ error->all(FLERR,"Invalid t_event in hyper command");
+ if (nsteps % t_event)
+ error->all(FLERR,"Hyper nsteps must be multiple of t_event");
+ if (rebond < 0)
+ error->all(FLERR,"Invalid rebond in hyper command");
+ if (rebond && rebond % t_event)
+ error->all(FLERR,"Hyper rebond must be multiple of t_event");
+
+ // FixHyper class performs global or local hyperdynamics
+
+ int hyperenable,hyperstyle;
+
+ if (strcmp(id_fix,"NULL") == 0) {
+ hyperenable = 0;
+ hyperstyle = NOHYPER;
+ } else {
+ int ifix = modify->find_fix(id_fix);
+ if (ifix < 0) error->all(FLERR,"Could not find fix ID for hyper");
+ fix_hyper = (FixHyper *) modify->fix[ifix];
+ int dim;
+ int *hyperflag = (int *) fix_hyper->extract("hyperflag",dim);
+ if (hyperflag == NULL || *hyperflag == 0)
+ error->all(FLERR,"Hyper fix is not a valid hyperdynamics fix");
+ if (*hyperflag == 1) hyperstyle = GLOBAL;
+ if (*hyperflag == 2) hyperstyle = LOCAL;
+ hyperenable = 1;
+ }
+
+ // create FixEventHyper class to store event and pre-quench states
+
+ char **args = new char*[3];
+ args[0] = (char *) "hyper_event";
+ args[1] = (char *) "all";
+ args[2] = (char *) "EVENT/HYPER";
+ modify->add_fix(3,args);
+ fix_event = (FixEventHyper *) modify->fix[modify->nfix-1];
+ delete [] args;
+
+ // create Finish for timing output
+
+ finish = new Finish(lmp);
+
+ // assign FixEventHyper to event-detection compute
+ // necessary so it will know atom coords at last event
+
+ int icompute = modify->find_compute(id_compute);
+ if (icompute < 0) error->all(FLERR,"Could not find compute ID for hyper");
+ compute_event = (ComputeEventDisplace *) modify->compute[icompute];
+ compute_event->reset_extra_compute_fix("hyper_event");
+
+ // reset reneighboring criteria since will perform minimizations
+
+ neigh_every = neighbor->every;
+ neigh_delay = neighbor->delay;
+ neigh_dist_check = neighbor->dist_check;
+
+ if (neigh_every != 1 || neigh_delay != 0 || neigh_dist_check != 1) {
+ if (me == 0)
+ error->warning(FLERR,"Resetting reneighboring criteria during hyper");
+ }
+
+ neighbor->every = 1;
+ neighbor->delay = 0;
+ neighbor->dist_check = 1;
+
+ // initialize hyper as if one long dynamics run
+
+ update->whichflag = 1;
+ update->nsteps = nsteps;
+ update->beginstep = update->firststep = update->ntimestep;
+ update->endstep = update->laststep = update->firststep + nsteps;
+ if (update->laststep < 0)
+ error->all(FLERR,"Too many timesteps");
+
+ lmp->init();
+
+ // init minimizer settings and minimizer itself
+
+ update->etol = etol;
+ update->ftol = ftol;
+ update->max_eval = maxeval;
+
+ // cannot use hyper with a changing box
+ // removing this restriction would require saving/restoring box params
+
+ if (domain->box_change)
+ error->all(FLERR,"Cannot use hyper with a changing box");
+
+ // cannot use hyper with time-dependent fixes or regions
+
+ for (int i = 0; i < modify->nfix; i++)
+ if (modify->fix[i]->time_depend)
+ error->all(FLERR,"Cannot use hyper with a time-dependent fix defined");
+
+ for (int i = 0; i < domain->nregion; i++)
+ if (domain->regions[i]->dynamic_check())
+ error->all(FLERR,"Cannot use hyper with a time-dependent region defined");
+
+ // perform hyperdynamics simulation
+
+ timer->init();
+ timer->barrier_start();
+ time_start = timer->get_wall(Timer::TOTAL);
+
+ nbuild = ndanger = 0;
+ time_dynamics = time_quench = 0.0;
+ double clock = 0.0;
+
+ if (hyperenable) fix_hyper->init_hyper();
+
+ timer->barrier_start();
+ time_start = timer->get_wall(Timer::TOTAL);
+
+ // perform initial minimization and bond list creation
+
+ int nevent = 0;
+ int nevent_atoms = 0;
+
+ fix_event->store_state_quench();
+ quench(1);
+ if (dumpflag)
+ for (int idump = 0; idump < ndump; idump++)
+ output->dump[dumplist[idump]]->write();
+ fix_event->store_event();
+ if (hyperenable) fix_hyper->build_bond_list(0);
+ fix_event->restore_state_quench();
+
+ // main loop: dynamics, store state, quench, check event, restore state
+
+ int ecount;
+ int istep = 0;
+
+ while (istep < nsteps) {
+ dynamics(t_event,time_dynamics);
+ fix_event->store_state_quench();
+ quench(0);
+
+ ecount = compute_event->all_events();
+
+ if (ecount) {
+ nevent++;
+ nevent_atoms += ecount;
+
+ if (dumpflag)
+ for (int idump = 0; idump < ndump; idump++)
+ output->dump[dumplist[idump]]->write();
+ fix_event->store_event();
+ if (hyperenable) fix_hyper->build_bond_list(ecount);
+
+ } else if (rebond && update->ntimestep % rebond == 0) {
+ fix_event->store_event();
+ if (hyperenable) fix_hyper->build_bond_list(ecount);
+ }
+
+ fix_event->restore_state_quench();
+
+ // NOTE: replace clock with hypertime
+ //clock = clock + t_event*universe->nworlds;
+ if (stepmode == 0) istep = update->ntimestep - update->beginstep;
+ else istep = clock;
+ }
+
+ if (stepmode) nsteps = update->ntimestep - update->beginstep;
+
+ // set total timers and counters so Finish() will process them
+
+ timer->set_wall(Timer::TOTAL,time_start);
+ timer->barrier_stop();
+
+ timer->set_wall(Timer::DYNAMICS,time_dynamics);
+ timer->set_wall(Timer::QUENCH,time_quench);
+
+ neighbor->ncalls = nbuild;
+ neighbor->ndanger = ndanger;
+
+ update->nsteps = nsteps;
+
+ if (me == 0) {
+ if (screen) fprintf(screen,"Final hyper stats ...\n\n");
+ if (logfile) fprintf(logfile,"Final hyper stats ...\n\n");
+ }
+
+ // subset of quantities also available in fix hyper output
+
+ int nevent_running = 0;
+ int nevent_atoms_running = 0;
+ double t_hyper = 0.0;
+ double avebonds = 0.0;
+ double maxdrift = 0.0;
+ double maxbondlen = 0.0;
+ double fraczero = 1.0;
+
+ double nnewbond,avenboost,aveboostcoeff,maxboostcoeff,minboostcoeff;
+ double maxbondperatom,neighbondperbond,aveboostnow;
+ double tbondbuild,rmaxever,rmaxeverbig,allghost_toofar;
+ double lostbond,lostbondcoeff,biasoverlap,nonmatchbiascoeff;
+
+ if (hyperenable) {
+ t_hyper = fix_hyper->query(1);
+ nevent_running = fix_hyper->query(2);
+ nevent_atoms_running = fix_hyper->query(3);
+ avebonds = fix_hyper->query(4);
+ maxdrift = fix_hyper->query(5);
+ maxbondlen = fix_hyper->query(6);
+ fraczero = fix_hyper->query(7);
+
+ if (hyperstyle == LOCAL) {
+ nnewbond = fix_hyper->query(8);
+ maxbondperatom = fix_hyper->query(9);
+ avenboost = fix_hyper->query(10);
+ aveboostcoeff = fix_hyper->query(11);
+ maxboostcoeff = fix_hyper->query(12);
+ minboostcoeff = fix_hyper->query(13);
+ neighbondperbond = fix_hyper->query(14);
+ aveboostnow = fix_hyper->query(15);
+ tbondbuild = fix_hyper->query(16);
+ rmaxever = fix_hyper->query(17);
+ rmaxeverbig = fix_hyper->query(18);
+ allghost_toofar = fix_hyper->query(19);
+ lostbond = fix_hyper->query(20);
+ lostbondcoeff = fix_hyper->query(21);
+ biasoverlap = fix_hyper->query(22);
+ nonmatchbiascoeff = fix_hyper->query(23);
+ }
+ }
+
+ if (me == 0) {
+ if (screen) {
+ fprintf(screen,"Cummulative quantities for fix hyper:\n");
+ fprintf(screen," hyper time = %g\n",t_hyper);
+ fprintf(screen," event timesteps = %d\n",nevent_running);
+ fprintf(screen," # of atoms in events = %d\n",nevent_atoms_running);
+ fprintf(screen,"Quantities for this hyper run:\n");
+ fprintf(screen," event timesteps = %d\n",nevent);
+ fprintf(screen," # of atoms in events = %d\n",nevent_atoms);
+ fprintf(screen," max length of any bond = %g\n",maxbondlen);
+ fprintf(screen," max drift distance of any atom = %g\n",maxdrift);
+ fprintf(screen," fraction of steps & bonds with zero bias = %g\n",
+ fraczero);
+ fprintf(screen,"Current quantities:\n");
+ fprintf(screen," ave bonds/atom = %g\n",avebonds);
+
+ if (hyperstyle == LOCAL) {
+ fprintf(screen,"Cummulative quantities specific to fix hyper/local:\n");
+ fprintf(screen," # of new bonds formed = %g\n",nnewbond);
+ fprintf(screen," max bonds/atom = %g\n",maxbondperatom);
+ fprintf(screen,"Quantities for this hyper run specific to "
+ "fix hyper/local:\n");
+ fprintf(screen," ave boosted bonds/step = %g\n",avenboost);
+ fprintf(screen," ave boost coeff of all bonds = %g\n",aveboostcoeff);
+ fprintf(screen," max boost coeff of any bond = %g\n",maxboostcoeff);
+ fprintf(screen," min boost coeff of any bond = %g\n",minboostcoeff);
+ fprintf(screen," max dist from my box of any "
+ "non-maxstrain bond ghost atom = %g\n",rmaxever);
+ fprintf(screen," max dist from my box of any bond ghost atom = %g\n",
+ rmaxeverbig);
+ fprintf(screen," count of bond ghost neighbors "
+ "not found on reneighbor steps = %g\n",allghost_toofar);
+ fprintf(screen," lost bond partners = %g\n",lostbond);
+ fprintf(screen," ave bias coeff for lost bond partners = %g\n",
+ lostbondcoeff);
+ fprintf(screen," bias overlaps = %g\n",biasoverlap);
+ fprintf(screen," non-matching bias coeffs = %g\n",nonmatchbiascoeff);
+ fprintf(screen," CPU time for bond builds = %g\n",tbondbuild);
+ fprintf(screen,"Current quantities specific to fix hyper/local:\n");
+ fprintf(screen," neighbor bonds/bond = %g\n",neighbondperbond);
+ fprintf(screen," ave boost coeff for all bonds = %g\n",aveboostnow);
+ }
+ fprintf(screen,"\n");
+ }
+
+ if (logfile) {
+ fprintf(logfile,"Cummulative quantities for fix hyper:\n");
+ fprintf(logfile," hyper time = %g\n",t_hyper);
+ fprintf(logfile," event timesteps = %d\n",nevent_running);
+ fprintf(logfile," # of atoms in events = %d\n",nevent_atoms_running);
+ fprintf(logfile,"Quantities for this hyper run:\n");
+ fprintf(logfile," event timesteps = %d\n",nevent);
+ fprintf(logfile," # of atoms in events = %d\n",nevent_atoms);
+ fprintf(logfile," max length of any bond = %g\n",maxbondlen);
+ fprintf(logfile," max drift distance of any atom = %g\n",maxdrift);
+ fprintf(logfile," fraction of steps & bonds with zero bias = %g\n",
+ fraczero);
+ fprintf(logfile,"Current quantities:\n");
+ fprintf(logfile," ave bonds/atom = %g\n",avebonds);
+
+ if (hyperstyle == LOCAL) {
+ fprintf(logfile,"Cummulative quantities specific tofix hyper/local:\n");
+ fprintf(logfile," # of new bonds formed = %g\n",nnewbond);
+ fprintf(logfile," max bonds/atom = %g\n",maxbondperatom);
+ fprintf(logfile,"Quantities for this hyper run specific to "
+ "fix hyper/local:\n");
+ fprintf(logfile," ave boosted bonds/step = %g\n",avenboost);
+ fprintf(logfile," ave boost coeff of all bonds = %g\n",aveboostcoeff);
+ fprintf(logfile," max boost coeff of any bond = %g\n",maxboostcoeff);
+ fprintf(logfile," min boost coeff of any bond = %g\n",minboostcoeff);
+ fprintf(logfile," max dist from my box of any "
+ "non-maxstrain bond ghost atom = %g\n",rmaxever);
+ fprintf(logfile," max dist from my box of any bond ghost atom = %g\n",
+ rmaxeverbig);
+ fprintf(logfile," count of ghost bond neighbors "
+ "not found on reneighbor steps = %g\n",allghost_toofar);
+ fprintf(logfile," lost bond partners = %g\n",lostbond);
+ fprintf(logfile," ave bias coeff for lost bond partners = %g\n",
+ lostbondcoeff);
+ fprintf(logfile," bias overlaps = %g\n",biasoverlap);
+ fprintf(logfile," non-matching bias coeffs = %g\n",nonmatchbiascoeff);
+ fprintf(logfile," CPU time for bond builds = %g\n",tbondbuild);
+ fprintf(logfile,"Current quantities specific to fix hyper/local:\n");
+ fprintf(logfile," neighbor bonds/bond = %g\n",neighbondperbond);
+ fprintf(logfile," ave boost coeff for all bonds = %g\n",aveboostnow);
+ }
+ fprintf(logfile,"\n");
+ }
+ }
+
+ // timing stats
+
+ finish->end(4);
+
+ update->whichflag = 0;
+ update->firststep = update->laststep = 0;
+ update->beginstep = update->endstep = 0;
+
+ // reset reneighboring criteria
+
+ neighbor->every = neigh_every;
+ neighbor->delay = neigh_delay;
+ neighbor->dist_check = neigh_dist_check;
+
+ delete [] id_fix;
+ delete [] id_compute;
+ memory->destroy(dumplist);
+ delete finish;
+ modify->delete_fix("hyper_event");
+
+ compute_event->reset_extra_compute_fix(NULL);
+}
+
+/* ----------------------------------------------------------------------
+ short dynamics run
+------------------------------------------------------------------------- */
+
+void Hyper::dynamics(int nsteps, double &time_category)
+{
+ update->whichflag = 1;
+ update->nsteps = nsteps;
+
+ // full init works
+ // need to try partial init or setup
+
+ 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;
+
+ timer->barrier_start();
+ update->integrate->run(nsteps);
+ timer->barrier_stop();
+ time_dynamics += timer->get_wall(Timer::TOTAL);
+
+ nbuild += neighbor->ncalls - ncalls;
+ ndanger += neighbor->ndanger;
+
+ update->integrate->cleanup();
+ finish->end(0);
+}
+
+/* ----------------------------------------------------------------------
+ quench minimization
+ flag = 1 to trigger output of memory in setup() call
+------------------------------------------------------------------------- */
+
+void Hyper::quench(int flag)
+{
+ bigint ntimestep_hold = update->ntimestep;
+ bigint endstep_hold = update->endstep;
+
+ // need to change whichflag so that minimize->setup() calling
+ // modify->setup() will call fix->min_setup()
+
+ update->whichflag = 2;
+ update->nsteps = maxiter;
+ update->endstep = update->laststep = update->firststep + maxiter;
+ if (update->laststep < 0)
+ error->all(FLERR,"Too many iterations");
+ update->restrict_output = 1;
+
+ // full init works
+
+ lmp->init();
+ update->minimize->setup(flag);
+
+ // partial init does not work
+
+ //modify->addstep_compute_all(update->ntimestep);
+ //update->minimize->setup_minimal(1);
+
+ // NOTE: what doing with ncalls?
+ int ncalls = neighbor->ncalls;
+
+ timer->barrier_start();
+ update->minimize->run(maxiter);
+ timer->barrier_stop();
+ time_quench += timer->get_wall(Timer::TOTAL);
+
+ update->minimize->cleanup();
+ finish->end(0);
+
+ // reset timestep as if quench did not occur
+ // clear timestep storage from computes, since now invalid
+
+ update->restrict_output = 0;
+ update->ntimestep = ntimestep_hold;
+ update->endstep = update->laststep = endstep_hold;
+ for (int i = 0; i < modify->ncompute; i++)
+ if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
+}
+
+/* ----------------------------------------------------------------------
+ parse optional parameters at end of hyper input line
+------------------------------------------------------------------------- */
+
+void Hyper::options(int narg, char **arg)
+{
+ // set defaults
+
+ etol = 1.0e-4;
+ ftol = 1.0e-4;
+ maxiter = 40;
+ maxeval = 50;
+ stepmode = 0;
+ dumpflag = 0;
+ ndump = 0;
+ dumplist = NULL;
+ rebond = 0;
+
+ int iarg = 0;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"min") == 0) {
+ if (iarg+5 > narg) error->all(FLERR,"Illegal hyper command");
+ etol = force->numeric(FLERR,arg[iarg+1]);
+ ftol = force->numeric(FLERR,arg[iarg+2]);
+ maxiter = force->inumeric(FLERR,arg[iarg+3]);
+ maxeval = force->inumeric(FLERR,arg[iarg+4]);
+ if (maxiter < 0) error->all(FLERR,"Illegal hyper command");
+ iarg += 5;
+
+ } else if (strcmp(arg[iarg],"time") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal hyper command");
+ if (strcmp(arg[iarg+1],"steps") == 0) stepmode = 0;
+ else if (strcmp(arg[iarg+1],"clock") == 0) stepmode = 1;
+ else error->all(FLERR,"Illegal hyper command");
+ iarg += 2;
+
+ } else if (strcmp(arg[iarg],"dump") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal hyper command");
+ dumpflag = 1;
+ int idump = output->find_dump(arg[iarg+1]);
+ if (idump < 0)
+ error->all(FLERR,"Dump ID in hyper command does not exist");
+ memory->grow(dumplist,ndump+1,"hyper:dumplist");
+ dumplist[ndump++] = idump;
+ iarg += 2;
+
+ } else if (strcmp(arg[iarg],"rebond") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal hyper command");
+ rebond = force->inumeric(FLERR,arg[iarg+1]);
+ iarg += 2;
+
+ } else error->all(FLERR,"Illegal hyper command");
+ }
+}
diff --git a/src/REPLICA/hyper.h b/src/REPLICA/hyper.h
new file mode 100644
index 0000000000..1b05172bcf
--- /dev/null
+++ b/src/REPLICA/hyper.h
@@ -0,0 +1,65 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMMAND_CLASS
+
+CommandStyle(hyper,Hyper)
+
+#else
+
+#ifndef LMP_HYPER_H
+#define LMP_HYPER_H
+
+#include "pointers.h"
+
+namespace LAMMPS_NS {
+
+class Hyper : protected Pointers {
+ public:
+ Hyper(class LAMMPS *);
+ ~Hyper() {}
+ void command(int, char **);
+
+ private:
+ int me,nprocs;
+ int t_event;
+ double etol,ftol;
+ int maxiter,maxeval;
+ int stepmode,dumpflag,ndump,rebond;
+ int *dumplist;
+
+ int neigh_every,neigh_delay,neigh_dist_check;
+ int quench_reneighbor;
+ bigint nbuild,ndanger;
+
+ double time_dynamics,time_quench;
+ double time_start;
+
+ class FixHyper *fix_hyper;
+ class FixEventHyper *fix_event;
+ class ComputeEventDisplace *compute_event;
+ class Finish *finish;
+
+ void dynamics(int, double &);
+ void quench(int flag);
+ void options(int, char **);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+*/
--
GitLab