-enable specific accelerator support via '-k on' "command-line switch"_Section_start.html#start_7, |
+enable specific accelerator support via '-k on' "command-line switch"_Section_start.html#start_6, |
only needed for KOKKOS package |
-set any needed options for the package via "-pk" "command-line switch"_Section_start.html#start_7 or "package"_package.html command, |
+set any needed options for the package via "-pk" "command-line switch"_Section_start.html#start_6 or "package"_package.html command, |
only if defaults need to be changed |
-use accelerated styles in your input via "-sf" "command-line switch"_Section_start.html#start_7 or "suffix"_suffix.html command | lmp_machine -in in.script -sf gpu
+use accelerated styles in your input via "-sf" "command-line switch"_Section_start.html#start_6 or "suffix"_suffix.html command | lmp_machine -in in.script -sf gpu
:tb(c=2,s=|)
Note that the first 4 steps can be done as a single command, using the
src/Make.py tool. This tool is discussed in "Section
-2.4"_Section_start.html#start_4 of the manual, and its use is
+4"_Section_packages.html of the manual, and its use is
illustrated in the individual accelerator sections. Typically these
steps only need to be done once, to create an executable that uses one
or more accelerator packages.
The last 4 steps can all be done from the command-line when LAMMPS is
launched, without changing your input script, as illustrated in the
individual accelerator sections. Or you can add
"package"_package.html and "suffix"_suffix.html commands to your input
script.
NOTE: With a few exceptions, you can build a single LAMMPS executable
with all its accelerator packages installed. Note however that the
USER-INTEL and KOKKOS packages require you to choose one of their
hardware options when building for a specific platform. I.e. CPU or
Phi option for the USER-INTEL package. Or the OpenMP, Cuda, or Phi
option for the KOKKOS package.
These are the exceptions. You cannot build a single executable with:
both the USER-INTEL Phi and KOKKOS Phi options
the USER-INTEL Phi or Kokkos Phi option, and the GPU package :ul
See the examples/accelerate/README and make.list files for sample
Make.py commands that build LAMMPS with any or all of the accelerator
packages. As an example, here is a command that builds with all the
GPU related packages installed (GPU, KOKKOS with Cuda), including
settings to build the needed auxiliary GPU libraries for Kepler GPUs:
This is not allowed if the kspace_modify overlap setting is no. :dd
{PPPM order < minimum allowed order} :dt
The default minimum order is 2. This can be reset by the
kspace_modify minorder command. :dd
{PPPM order cannot be < 2 or > than %d} :dt
This is a limitation of the PPPM implementation in LAMMPS. :dd
{PPPMDisp Coulomb grid is too large} :dt
The global PPPM grid is larger than OFFSET in one or more dimensions.
OFFSET is currently set to 4096. You likely need to decrease the
requested accuracy. :dd
{PPPMDisp Dispersion grid is too large} :dt
The global PPPM grid is larger than OFFSET in one or more dimensions.
OFFSET is currently set to 4096. You likely need to decrease the
requested accuracy. :dd
{PPPMDisp can only currently be used with comm_style brick} :dt
This is a current restriction in LAMMPS. :dd
{PPPMDisp coulomb order cannot be greater than %d} :dt
This is a limitation of the PPPM implementation in LAMMPS. :dd
{PPPMDisp used but no parameters set, for further information please see the pppm/disp documentation} :dt
An efficient and accurate usage of the pppm/disp requires settings via the kspace_modify command. Please see the pppm/disp documentation for further instructions. :dd
{PRD command before simulation box is defined} :dt
The prd command cannot be used before a read_data,
read_restart, or create_box command. :dd
{PRD nsteps must be multiple of t_event} :dt
Self-explanatory. :dd
{PRD t_corr must be multiple of t_event} :dt
Self-explanatory. :dd
{Package command after simulation box is defined} :dt
The package command cannot be used afer a read_data, read_restart, or
create_box command. :dd
{Package cuda command without USER-CUDA package enabled} :dt
The USER-CUDA package must be installed via "make yes-user-cuda"
before LAMMPS is built, and the "-c on" must be used to enable the
package. :dd
{Package gpu command without GPU package installed} :dt
The GPU package must be installed via "make yes-gpu" before LAMMPS is
built. :dd
{Package intel command without USER-INTEL package installed} :dt
The USER-INTEL package must be installed via "make yes-user-intel"
before LAMMPS is built. :dd
{Package kokkos command without KOKKOS package enabled} :dt
The KOKKOS package must be installed via "make yes-kokkos" before
LAMMPS is built, and the "-k on" must be used to enable the package. :dd
{Package omp command without USER-OMP package installed} :dt
The USER-OMP package must be installed via "make yes-user-omp" before
LAMMPS is built. :dd
{Pair body requires atom style body} :dt
Self-explanatory. :dd
{Pair body requires body style nparticle} :dt
This pair style is specific to the nparticle body style. :dd
This may lead to a larger grid than desired. See the kspace_modify overlap
command to prevent changing of the PPPM order. :dd
{Replacing a fix, but new group != old group} :dt
The ID and style of a fix match for a fix you are changing with a fix
command, but the new group you are specifying does not match the old
group. :dd
{Replicating in a non-periodic dimension} :dt
The parameters for a replicate command will cause a non-periodic
dimension to be replicated; this may cause unwanted behavior. :dd
{Resetting reneighboring criteria during PRD} :dt
A PRD simulation requires that neigh_modify settings be delay = 0,
every = 1, check = yes. Since these settings were not in place,
LAMMPS changed them and will restore them to their original values
after the PRD simulation. :dd
{Resetting reneighboring criteria during TAD} :dt
A TAD simulation requires that neigh_modify settings be delay = 0,
every = 1, check = yes. Since these settings were not in place,
LAMMPS changed them and will restore them to their original values
after the PRD simulation. :dd
{Resetting reneighboring criteria during minimization} :dt
Minimization requires that neigh_modify settings be delay = 0, every =
1, check = yes. Since these settings were not in place, LAMMPS
changed them and will restore them to their original values after the
minimization. :dd
{Restart file used different # of processors} :dt
The restart file was written out by a LAMMPS simulation running on a
different number of processors. Due to round-off, the trajectories of
your restarted simulation may diverge a little more quickly than if
you ran on the same # of processors. :dd
{Restart file used different 3d processor grid} :dt
The restart file was written out by a LAMMPS simulation running on a
different 3d grid of processors. Due to round-off, the trajectories
of your restarted simulation may diverge a little more quickly than if
you ran on the same # of processors. :dd
{Restart file used different boundary settings, using restart file values} :dt
Your input script cannot change these restart file settings. :dd
{Restart file used different newton bond setting, using restart file value} :dt
The restart file value will override the setting in the input script. :dd
{Restart file used different newton pair setting, using input script value} :dt
The input script value will override the setting in the restart file. :dd
{Restrain problem: %d %ld %d %d %d %d} :dt
Conformation of the 4 listed dihedral atoms is extreme; you may want
to check your simulation geometry. :dd
{Running PRD with only one replica} :dt
This is allowed, but you will get no parallel speed-up. :dd
{SRD bin shifting turned on due to small lamda} :dt
This is done to try to preserve accuracy. :dd
{SRD bin size for fix srd differs from user request} :dt
Fix SRD had to adjust the bin size to fit the simulation box. See the
cubic keyword if you want this message to be an error vs warning. :dd
{SRD bins for fix srd are not cubic enough} :dt
The bin shape is not within tolerance of cubic. See the cubic
keyword if you want this message to be an error vs warning. :dd
{SRD particle %d started inside big particle %d on step %ld bounce %d} :dt
See the inside keyword if you want this message to be an error vs
warning. :dd
{SRD particle %d started inside wall %d on step %ld bounce %d} :dt
See the inside keyword if you want this message to be an error vs
warning. :dd
{Shake determinant < 0.0} :dt
The determinant of the quadratic equation being solved for a single
cluster specified by the fix shake command is numerically suspect. LAMMPS
will set it to 0.0 and continue. :dd
{Shell command '%s' failed with error '%s'} :dt
Self-explanatory. :dd
{Shell command returned with non-zero status} :dt
This may indicate the shell command did not operate as expected. :dd
{Should not allow rigid bodies to bounce off relecting walls} :dt
LAMMPS allows this, but their dynamics are not computed correctly. :dd
{Should not use fix nve/limit with fix shake or fix rattle} :dt
This will lead to invalid constraint forces in the SHAKE/RATTLE
computation. :dd
{Simulations might be very slow because of large number of structure factors} :dt
Self-explanatory. :dd
{Slab correction not needed for MSM} :dt
Slab correction is intended to be used with Ewald or PPPM and is not needed by MSM. :dd
{System is not charge neutral, net charge = %g} :dt
The total charge on all atoms on the system is not 0.0.
For some KSpace solvers this is only a warning. :dd
{Table inner cutoff >= outer cutoff} :dt
You specified an inner cutoff for a Coulombic table that is longer
than the global cutoff. Probably not what you wanted. :dd
{Temperature for MSST is not for group all} :dt
User-assigned temperature to MSST fix does not compute temperature for
all atoms. Since MSST computes a global pressure, the kinetic energy
contribution from the temperature is assumed to also be for all atoms.
Thus the pressure used by MSST could be inaccurate. :dd
{Temperature for NPT is not for group all} :dt
User-assigned temperature to NPT fix does not compute temperature for
all atoms. Since NPT computes a global pressure, the kinetic energy
contribution from the temperature is assumed to also be for all atoms.
Thus the pressure used by NPT could be inaccurate. :dd
{Temperature for fix modify is not for group all} :dt
The temperature compute is being used with a pressure calculation
which does operate on group all, so this may be inconsistent. :dd
{Temperature for thermo pressure is not for group all} :dt
User-assigned temperature to thermo via the thermo_modify command does
not compute temperature for all atoms. Since thermo computes a global
pressure, the kinetic energy contribution from the temperature is
assumed to also be for all atoms. Thus the pressure printed by thermo
could be inaccurate. :dd
{The fix ave/spatial command has been replaced by the more flexible fix ave/chunk and compute chunk/atom commands -- fix ave/spatial will be removed in the summer of 2015} :dt
Self-explanatory. :dd
{The minimizer does not re-orient dipoles when using fix efield} :dt
This means that only the atom coordinates will be minimized,
not the orientation of the dipoles. :dd
{Too many common neighbors in CNA %d times} :dt
More than the maximum # of neighbors was found multiple times. This
was unexpected. :dd
{Too many inner timesteps in fix ttm} :dt
Self-explanatory. :dd
{Too many neighbors in CNA for %d atoms} :dt
More than the maximum # of neighbors was found multiple times. This
was unexpected. :dd
{Triclinic box skew is large} :dt
The displacement in a skewed direction is normally required to be less
than half the box length in that dimension. E.g. the xy tilt must be
between -half and +half of the x box length. You have relaxed the
constraint using the box tilt command, but the warning means that a
LAMMPS simulation may be inefficient as a result. :dd
{Use special bonds = 0,1,1 with bond style fene} :dt
Most FENE models need this setting for the special_bonds command. :dd
{Use special bonds = 0,1,1 with bond style fene/expand} :dt
Most FENE models need this setting for the special_bonds command. :dd
{Using a manybody potential with bonds/angles/dihedrals and special_bond exclusions} :dt
This is likely not what you want to do. The exclusion settings will
eliminate neighbors in the neighbor list, which the manybody potential
needs to calculated its terms correctly. :dd
{Using compute temp/deform with inconsistent fix deform remap option} :dt
Fix nvt/sllod assumes deforming atoms have a velocity profile provided
by "remap v" or "remap none" as a fix deform option. :dd
{Using compute temp/deform with no fix deform defined} :dt
This is probably an error, since it makes little sense to use
compute temp/deform in this case. :dd
{Using fix srd with box deformation but no SRD thermostat} :dt
The deformation will heat the SRD particles so this can
be dangerous. :dd
{Using kspace solver on system with no charge} :dt
Self-explanatory. :dd
{Using largest cut-off for lj/long/dipole/long long long} :dt
Self-explanatory. :dd
{Using largest cutoff for buck/long/coul/long} :dt
Self-explanatory. :dd
{Using largest cutoff for lj/long/coul/long} :dt
Self-explanatory. :dd
{Using largest cutoff for pair_style lj/long/tip4p/long} :dt
Self-explanatory. :dd
{Using package gpu without any pair style defined} :dt
Self-explanatory. :dd
{Using pair potential shift with pair_modify compute no} :dt
The shift effects will thus not be computed. :dd
{Using pair tail corrections with nonperiodic system} :dt
This is probably a bogus thing to do, since tail corrections are
computed by integrating the density of a periodic system out to
infinity. :dd
{Using pair tail corrections with pair_modify compute no} :dt
The tail corrections will thus not be computed. :dd
{pair style reax is now deprecated and will soon be retired. Users should switch to pair_style reax/c} :dt
6.22 Calculating a diffusion coefficient :link(howto_22),h4
The diffusion coefficient D of a material can be measured in at least
2 ways using various options in LAMMPS. See the examples/DIFFUSE
directory for scripts that implement the 2 methods discussed here for
a simple Lennard-Jones fluid model.
The first method is to measure the mean-squared displacement (MSD) of
the system, via the "compute msd"_compute_msd.html command. The slope
of the MSD versus time is proportional to the diffusion coefficient.
The instantaneous MSD values can be accumulated in a vector via the
"fix vector"_fix_vector.html command, and a line fit to the vector to
compute its slope via the "variable slope"_variable.html function, and
thus extract D.
The second method is to measure the velocity auto-correlation function
(VACF) of the system, via the "compute vacf"_compute_vacf.html
command. The time-integral of the VACF is proportional to the
diffusion coefficient. The instantaneous VACF values can be
accumulated in a vector via the "fix vector"_fix_vector.html command,
and time integrated via the "variable trap"_variable.html function,
and thus extract D.
:line
6.23 Using chunks to calculate system properties :link(howto_23),h4
In LAMMS, "chunks" are collections of atoms, as defined by the
"compute chunk/atom"_compute_chunk_atom.html command, which assigns
each atom to a chunk ID (or to no chunk at all). The number of chunks
and the assignment of chunk IDs to atoms can be static or change over
time. Examples of "chunks" are molecules or spatial bins or atoms
with similar values (e.g. coordination number or potential energy).
The per-atom chunk IDs can be used as input to two other kinds of
commands, to calculate various properties of a system:
"fix ave/chunk"_fix_ave_chunk.html
any of the "compute */chunk"_compute.html commands :ul
Here, each of the 3 kinds of chunk-related commands is briefly
overviewed. Then some examples are given of how to compute different
properties with chunk commands.
Compute chunk/atom command: :h5
This compute can assign atoms to chunks of various styles. Only atoms
in the specified group and optional specified region are assigned to a
chunk. Here are some possible chunk definitions:
atoms in same molecule | chunk ID = molecule ID |
atoms of same atom type | chunk ID = atom type |
all atoms with same atom property (charge, radius, etc) | chunk ID = output of compute property/atom |
atoms in same cluster | chunk ID = output of "compute cluster/atom"_compute_cluster_atom.html command |
atoms in same spatial bin | chunk ID = bin ID |
atoms in same rigid body | chunk ID = molecule ID used to define rigid bodies |
atoms with similar potential energy | chunk ID = output of "compute pe/atom"_compute_pe_atom.html |
atoms with same local defect structure | chunk ID = output of "compute centro/atom"_compute_centro_atom.html or "compute coord/atom"_compute_coord_atom.html command :tb(s=|,c=2)
Note that chunk IDs are integer values, so for atom properties or
computes that produce a floating point value, they will be truncated
to an integer. You could also use the compute in a variable that
scales the floating point value to spread it across multiple integers.
Spatial bins can be of various kinds, e.g. 1d bins = slabs, 2d bins =
pencils, 3d bins = boxes, spherical bins, cylindrical bins.
This compute also calculates the number of chunks {Nchunk}, which is
used by other commands to tally per-chunk data. {Nchunk} can be a
static value or change over time (e.g. the number of clusters). The
chunk ID for an individual atom can also be static (e.g. a molecule
ID), or dynamic (e.g. what spatial bin an atom is in as it moves).
Note that this compute allows the per-atom output of other
"computes"_compute.html, "fixes"_fix.html, and
"variables"_variable.html to be used to define chunk IDs for each
atom. This means you can write your own compute or fix to output a
per-atom quantity to use as chunk ID. See
"Section 10"_Section_modify.html of the documentation for how to
do this. You can also define a "per-atom variable"_variable.html in
the input script that uses a formula to generate a chunk ID for each
atom.
Fix ave/chunk command: :h5
This fix takes the ID of a "compute
chunk/atom"_compute_chunk_atom.html command as input. For each chunk,
it then sums one or more specified per-atom values over the atoms in
each chunk. The per-atom values can be any atom property, such as
mpirun -np 72 -ppn 36 lmp_machine -sf intel -in in.script -pk intel 0 omp 2 mode double # Don't use any coprocessors that might be available, use 2 OpenMP threads for each task, use double precision :pre
[Or run with the USER-INTEL package by editing an input script:]
As an alternative to adding command-line arguments, the input script
can be edited to enable the USER-INTEL package. This requires adding
the "package intel"_package.html command to the top of the input
script. For the second example above, this would be:
package intel 0 omp 2 mode double :pre
To enable the USER-INTEL package only for individual styles, you can
add an "intel" suffix to the individual style, e.g.:
pair_style lj/cut/intel 2.5 :pre
Alternatively, the "suffix intel"_suffix.html command can be added to
the input script to enable USER-INTEL styles for the commands that
follow in the input script.
[Tuning for Performance:]
NOTE: The USER-INTEL package will perform better with modifications
to the input script when "PPPM"_kspace_style.html is used:
"kspace_modify diff ad"_kspace_modify.html should be added to the
input script.
Long-Range Thread (LRT) mode is an option to the "package
intel"_package.html command that can improve performance when using
"PPPM"_kspace_style.html for long-range electrostatics on processors
with SMT. It generates an extra pthread for each MPI task. The thread
is dedicated to performing some of the PPPM calculations and MPI
communications. On Intel Xeon Phi x200 series CPUs, this will likely
always improve performance, even on a single node. On Intel Xeon
processors, using this mode might result in better performance when
using multiple nodes, depending on the machine. To use this mode,
specify that the number of OpenMP threads is one less than would
normally be used for the run and add the "lrt yes" option to the "-pk"
command-line suffix or "package intel" command. For example, if a run
would normally perform best with "-pk intel 0 omp 4", instead use
"-pk intel 0 omp 3 lrt yes". When using LRT, you should set the
environment variable "KMP_AFFINITY=none". LRT mode is not supported
when using offload.
NOTE: Changing the "newton"_newton.html setting to off can improve
performance and/or scalability for simple 2-body potentials such as
lj/cut or when using LRT mode on processors supporting AVX-512.
Not all styles are supported in the USER-INTEL package. You can mix
the USER-INTEL package with styles from the "OPT"_accelerate_opt.html
package or the "USER-OMP package"_accelerate_omp.html. Of course,
this requires that these packages were installed at build time. This
can performed automatically by using "-sf hybrid intel opt" or
"-sf hybrid intel omp" command-line options. Alternatively, the "opt"
and "omp" suffixes can be appended manually in the input script. For
the latter, the "package omp"_package.html command must be in the
input script or the "-pk omp Nt" "command-line
-switch"_Section_start.html#start_7 must be used where Nt is the
+switch"_Section_start.html#start_6 must be used where Nt is the
number of OpenMP threads. The number of OpenMP threads should not be
set differently for the different packages. Note that the "suffix
hybrid intel omp"_suffix.html command can also be used within the
input script to automatically append the "omp" suffix to styles when
USER-INTEL styles are not available.
NOTE: For simulations on higher node counts, add "processors * * *
grid numa"_processors.html" to the beginning of the input script for
better scalability.
When running on many nodes, performance might be better when using
fewer OpenMP threads and more MPI tasks. This will depend on the
simulation and the machine. Using the "verlet/split"_run_style.html
run style might also give better performance for simulations with
"PPPM"_kspace_style.html electrostatics. Note that this is an
alternative to LRT mode and the two cannot be used together.
Currently, when using Intel MPI with Intel Xeon Phi x200 series
CPUs, better performance might be obtained by setting the
environment variable "I_MPI_SHM_LMT=shm" for Linux kernels that do
not yet have full support for AVX-512. Runs on Intel Xeon Phi x200
series processors will always perform better using MCDRAM. Please
consult your system documentation for the best approach to specify
that MPI runs are performed in MCDRAM.
[Tuning for Offload Performance:]
The default settings for offload should give good performance.
When using LAMMPS with offload to Intel coprocessors, best performance
will typically be achieved with concurrent calculations performed on
both the CPU and the coprocessor. This is achieved by offloading only
a fraction of the neighbor and pair computations to the coprocessor or
using "hybrid"_pair_hybrid.html pair styles where only one style uses
the "intel" suffix. For simulations with long-range electrostatics or
bond, angle, dihedral, improper calculations, computation and data
transfer to the coprocessor will run concurrently with computations
and MPI communications for these calculations on the host CPU. This
is illustrated in the figure below for the rhodopsin protein benchmark
running on E5-2697v2 processors with a Intel Xeon Phi 7120p
coprocessor. In this plot, the vertical access is time and routines
running at the same time are running concurrently on both the host and
the coprocessor.
:c,image(JPG/offload_knc.png)
The fraction of the offloaded work is controlled by the {balance}
keyword in the "package intel"_package.html command. A balance of 0
runs all calculations on the CPU. A balance of 1 runs all
supported calculations on the coprocessor. A balance of 0.5 runs half
of the calculations on the coprocessor. Setting the balance to -1
(the default) will enable dynamic load balancing that continously
adjusts the fraction of offloaded work throughout the simulation.
Because data transfer cannot be timed, this option typically produces
results within 5 to 10 percent of the optimal fixed balance.
If running short benchmark runs with dynamic load balancing, adding a
short warm-up run (10-20 steps) will allow the load-balancer to find a
near-optimal setting that will carry over to additional runs.
The default for the "package intel"_package.html command is to have
all the MPI tasks on a given compute node use a single Xeon Phi
coprocessor. In general, running with a large number of MPI tasks on
each node will perform best with offload. Each MPI task will
automatically get affinity to a subset of the hardware threads
available on the coprocessor. For example, if your card has 61 cores,
with 60 cores available for offload and 4 hardware threads per core
(240 total threads), running with 24 MPI tasks per node will cause
each MPI task to use a subset of 10 threads on the coprocessor. Fine
tuning of the number of threads to use per MPI task or the number of
threads to use per core can be accomplished with keyword settings of
the "package intel"_package.html command.
The USER-INTEL package has two modes for deciding which atoms will be
handled by the coprocessor. This choice is controlled with the {ghost}
keyword of the "package intel"_package.html command. When set to 0,
ghost atoms (atoms at the borders between MPI tasks) are not offloaded
to the card. This allows for overlap of MPI communication of forces
with computation on the coprocessor when the "newton"_newton.html
setting is "on". The default is dependent on the style being used,
however, better performance may be achieved by setting this option
explicitly.
When using offload with CPU Hyper-Threading disabled, it may help
performance to use fewer MPI tasks and OpenMP threads than available
cores. This is due to the fact that additional threads are generated
internally to handle the asynchronous offload tasks.
If pair computations are being offloaded to an Intel Xeon Phi
coprocessor, a diagnostic line is printed to the screen (not to the
log file), during the setup phase of a run, indicating that offload
mode is being used and indicating the number of coprocessor threads
per MPI task. Additionally, an offload timing summary is printed at
the end of each run. When offloading, the frequency for "atom
sorting"_atom_modify.html is changed to 1 so that the per-atom data is
effectively sorted at every rebuild of the neighbor lists. All the
available coprocessor threads on each Phi will be divided among MPI
tasks, unless the {tptask} option of the "-pk intel" "command-line
-switch"_Section_start.html#start_7 is used to limit the coprocessor
+switch"_Section_start.html#start_6 is used to limit the coprocessor
threads per MPI task.
[Restrictions:]
When offloading to a coprocessor, "hybrid"_pair_hybrid.html styles
that require skip lists for neighbor builds cannot be offloaded.
Using "hybrid/overlay"_pair_hybrid.html is allowed. Only one intel
accelerated style may be used with hybrid styles.
"Special_bonds"_special_bonds.html exclusion lists are not currently
supported with offload, however, the same effect can often be
accomplished by setting cutoffs for excluded atom types to 0. None of
the pair styles in the USER-INTEL package currently support the
"inner", "middle", "outer" options for rRESPA integration via the
"run_style respa"_run_style.html command; only the "pair" option is
supported.
[References:]
Brown, W.M., Carrillo, J.-M.Y., Mishra, B., Gavhane, N., Thakker, F.M., De Kraker, A.R., Yamada, M., Ang, J.A., Plimpton, S.J., "Optimizing Classical Molecular Dynamics in LAMMPS," in Intel Xeon Phi Processor High Performance Programming: Knights Landing Edition, J. Jeffers, J. Reinders, A. Sodani, Eds. Morgan Kaufmann. :ulb,l
Brown, W. M., Semin, A., Hebenstreit, M., Khvostov, S., Raman, K., Plimpton, S.J. "Increasing Molecular Dynamics Simulation Rates with an 8-Fold Increase in Electrical Power Efficiency."_http://dl.acm.org/citation.cfm?id=3014915 2016 High Performance Computing, Networking, Storage and Analysis, SC16: International Conference (pp. 82-95). :l
{dhugoniot} is the departure from the Hugoniot (temperature units).
{drayleigh} is the departure from the Rayleigh line (pressure units).
{lagrangian_speed} is the laboratory-frame Lagrangian speed (particle velocity) of the computational cell (velocity units).
{lagrangian_position} is the computational cell position in the reference frame moving at the shock speed. This is usually a good estimate of distance of the computational cell behind the shock front. :ol
To print these quantities to the log file with descriptive column
headers, the following LAMMPS commands are suggested:
fix msst all msst z
fix_modify msst energy yes
variable dhug equal f_msst\[1\]
variable dray equal f_msst\[2\]
variable lgr_vel equal f_msst\[3\]
variable lgr_pos equal f_msst\[4\]
thermo_style custom step temp ke pe lz pzz etotal v_dhug v_dray v_lgr_vel v_lgr_pos f_msst :pre
These fixes compute a global scalar and a global vector of 4
quantities, which can be accessed by various "output
commands"_Section_howto.html#howto_15. The scalar values calculated
by this fix are "extensive"; the vector values are "intensive".
[Restrictions:]
This fix style is part of the SHOCK 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.
All cell dimensions must be periodic. This fix can not be used with a
triclinic cell. The MSST fix has been tested only for the group-ID
ID, group-ID are documented in "fix"_fix.html command :ulb,l
style_name = {nvt} or {npt} or {nph} :l
one or more keyword/value pairs may be appended :l
keyword = {temp} or {iso} or {aniso} or {tri} or {x} or {y} or {z} or {xy} or {yz} or {xz} or {couple} or {tchain} or {pchain} or {mtk} or {tloop} or {ploop} or {nreset} or {drag} or {dilate} or {scalexy} or {scaleyz} or {scalexz} or {flip} or {fixedpoint} or {update}
{temp} values = Tstart Tstop Tdamp
Tstart,Tstop = external temperature at start/end of run
Tdamp = temperature damping parameter (time units)
{iso} or {aniso} or {tri} values = Pstart Pstop Pdamp
Pstart,Pstop = scalar external pressure at start/end of run (pressure units)
Pdamp = pressure damping parameter (time units)
{x} or {y} or {z} or {xy} or {yz} or {xz} values = Pstart Pstop Pdamp
Pstart,Pstop = external stress tensor component at start/end of run (pressure units)
Pdamp = stress damping parameter (time units)
{couple} = {none} or {xyz} or {xy} or {yz} or {xz}
{tchain} value = N
N = length of thermostat chain (1 = single thermostat)
{pchain} values = N
N length of thermostat chain on barostat (0 = no thermostat)
{mtk} value = {yes} or {no} = add in MTK adjustment term or not
{tloop} value = M
M = number of sub-cycles to perform on thermostat
{ploop} value = M
M = number of sub-cycles to perform on barostat thermostat
{nreset} value = reset reference cell every this many timesteps
{drag} value = Df
Df = drag factor added to barostat/thermostat (0.0 = no drag)
{dilate} value = dilate-group-ID
dilate-group-ID = only dilate atoms in this group due to barostat volume changes
{scalexy} value = {yes} or {no} = scale xy with ly
{scaleyz} value = {yes} or {no} = scale yz with lz
{scalexz} value = {yes} or {no} = scale xz with lz
{flip} value = {yes} or {no} = allow or disallow box flips when it becomes highly skewed
{fixedpoint} values = x y z
x,y,z = perform barostat dilation/contraction around this point (distance units)
{update} value = {dipole} or {dipole/dlm}
dipole = update dipole orientation (only for sphere variants)
dipole/dlm = use DLM integrator to update dipole orientation (only for sphere variants) :pre
:ule
[Examples:]
fix 1 all nvt temp 300.0 300.0 100.0
fix 1 water npt temp 300.0 300.0 100.0 iso 0.0 0.0 1000.0
ID, group-ID are documented in "fix"_fix.html command :ulb,l
one or more keyword value pairs may be appended
keyword = {temp} or {iso} or {aniso} or {tri} or {x} or {y} or {z} or {couple} or {tchain} or {pchain} or {mtk} or {tloop} or {ploop} or {nreset} or {drag} or {dilate} or {scaleyz} or {scalexz} or {scalexy}
{temp} values = Value1 Value2 Tdamp
Value1, Value2 = Nose-Hoover target temperatures, ignored by Hugoniostat
Tdamp = temperature damping parameter (time units)
{iso} or {aniso} or {tri} values = Pstart Pstop Pdamp
Pstart,Pstop = scalar external pressures, must be equal (pressure units)
Pdamp = pressure damping parameter (time units)
{x} or {y} or {z} or {xy} or {yz} or {xz} values = Pstart Pstop Pdamp
Pstart,Pstop = external stress tensor components, must be equal (pressure units)
Pdamp = stress damping parameter (time units)
{couple} = {none} or {xyz} or {xy} or {yz} or {xz}
{tchain} value = length of thermostat chain (1 = single thermostat)
{pchain} values = length of thermostat chain on barostat (0 = no thermostat)
{mtk} value = {yes} or {no} = add in MTK adjustment term or not
{tloop} value = number of sub-cycles to perform on thermostat
{ploop} value = number of sub-cycles to perform on barostat thermostat
{nreset} value = reset reference cell every this many timesteps
{drag} value = drag factor added to barostat/thermostat (0.0 = no drag)
{dilate} value = {all} or {partial}
{scaleyz} value = {yes} or {no} = scale yz with lz
{scalexz} value = {yes} or {no} = scale xz with lz
{scalexy} value = {yes} or {no} = scale xy with ly :pre
:ule
[Examples:]
fix myhug all nphug temp 1.0 1.0 10.0 z 40.0 40.0 70.0
fix myhug all nphug temp 1.0 1.0 10.0 iso 40.0 40.0 70.0 drag 200.0 tchain 1 pchain 0 :pre
[Description:]
This command is a variant of the Nose-Hoover
"fix npt"_fix_nh.html fix style.
It performs time integration of the Hugoniostat equations
of motion developed by Ravelo et al. "(Ravelo)"_#Ravelo1.
These equations compress the system to a state with average
axial stress or pressure equal to the specified target value
and that satisfies the Rankine-Hugoniot (RH)
jump conditions for steady shocks.
The compression can be performed
either
hydrostatically (using keyword {iso}, {aniso}, or {tri}) or uniaxially
(using keywords {x}, {y}, or {z}). In the hydrostatic case,
the cell dimensions change dynamically so that the average axial stress
in all three directions converges towards the specified target value.
In the uniaxial case, the chosen cell dimension changes dynamically
so that the average
axial stress in that direction converges towards the target value. The
other two cell dimensions are kept fixed (zero lateral strain).
This leads to the following additional restrictions on the keywords:
One and only one of the following keywords should be used: {iso}, {aniso}, {tri}, {x}, {y}, {z}
The specified initial and final target pressures must be the same.
The keywords {xy}, {xz}, {yz} may not be used.
The only admissible value for the couple keyword is {xyz}, which has the same effect as keyword {iso}
The {temp} keyword must be used to specify the time constant for kinetic energy relaxation, but initial and final target temperature values are ignored. :ul
Essentially, a Hugoniostat simulation is an NPT simulation in which the
user-specified target temperature is replaced with a time-dependent
target temperature Tt obtained from the following equation:
:c,image(Eqs/fix_nphug.jpg)
where T and Tt are the instantaneous and target temperatures,
P and P0 are the instantaneous and reference pressures or axial stresses,
depending on whether hydrostatic or uniaxial compression is being
performed, V and V0 are the instantaneous and reference volumes,
E and E0 are the instantaneous and reference internal energy (potential
plus kinetic), Ndof is the number of degrees of freedom used in the
definition of temperature, and kB is the Boltzmann constant. Delta is the
negative deviation of the instantaneous temperature from the target temperature.
When the system reaches a stable equilibrium, the value of Delta should
fluctuate about zero.
The values of E0, V0, and P0 are the instantaneous values at the start of
the simulation. These can be overridden using the fix_modify keywords {e0},
{v0}, and {p0} described below.
:line
NOTE: Unlike the "fix temp/berendsen"_fix_temp_berendsen.html command
which performs thermostatting but NO time integration, this fix
performs thermostatting/barostatting AND time integration. Thus you
should not use any other time integration fix, such as "fix
nve"_fix_nve.html on atoms to which this fix is applied. Likewise,
this fix should not be used on atoms that have their temperature
controlled by another fix - e.g. by "fix langevin"_fix_nh.html or "fix
temp/rescale"_fix_temp_rescale.html commands.
:line
This fix computes a temperature and pressure at each timestep. To do
this, the fix creates its own computes of style "temp" and "pressure",
as if one of these two sets of commands had been issued:
fix ID group-ID style bodystyle args keyword values ... :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
style = {rigid} or {rigid/nve} or {rigid/nvt} or {rigid/npt} or {rigid/nph} or {rigid/small} or {rigid/nve/small} or {rigid/nvt/small} or {rigid/npt/small} or {rigid/nph/small} :l
bodystyle = {single} or {molecule} or {group} :l
{single} args = none
{molecule} args = none
{group} args = N groupID1 groupID2 ...
N = # of groups
groupID1, groupID2, ... = list of N group IDs :pre
zero or more keyword/value pairs may be appended :l
keyword = {langevin} or {reinit} or {temp} or {iso} or {aniso} or {x} or {y} or {z} or {couple} or {tparam} or {pchain} or {dilate} or {force} or {torque} or {infile} :l
{langevin} values = Tstart Tstop Tperiod seed
Tstart,Tstop = desired temperature at start/stop of run (temperature units)
Tdamp = temperature damping parameter (time units)
seed = random number seed to use for white noise (positive integer)
{reinit} = {yes} or {no}
{temp} values = Tstart Tstop Tdamp
Tstart,Tstop = desired temperature at start/stop of run (temperature units)
Tdamp = temperature damping parameter (time units)
{iso} or {aniso} values = Pstart Pstop Pdamp
Pstart,Pstop = scalar external pressure at start/end of run (pressure units)
Pdamp = pressure damping parameter (time units)
{x} or {y} or {z} values = Pstart Pstop Pdamp
Pstart,Pstop = external stress tensor component at start/end of run (pressure units)
Pdamp = stress damping parameter (time units)
{couple} = {none} or {xyz} or {xy} or {yz} or {xz}
{tparam} values = Tchain Titer Torder
Tchain = length of Nose/Hoover thermostat chain
Titer = number of thermostat iterations performed
Torder = 3 or 5 = Yoshida-Suzuki integration parameters
{pchain} values = Pchain
Pchain = length of the Nose/Hoover thermostat chain coupled with the barostat
{dilate} value = dilate-group-ID
dilate-group-ID = only dilate atoms in this group due to barostat volume changes
{force} values = M xflag yflag zflag
M = which rigid body from 1-Nbody (see asterisk form below)
xflag,yflag,zflag = off/on if component of center-of-mass force is active
{torque} values = M xflag yflag zflag
M = which rigid body from 1-Nbody (see asterisk form below)
xflag,yflag,zflag = off/on if component of center-of-mass torque is active
{infile} filename
filename = file with per-body values of mass, center-of-mass, moments of inertia
{mol} value = template-ID
template-ID = ID of molecule template specified in a separate "molecule"_molecule.html command :pre
:ule
[Examples:]
fix 1 clump rigid single reinit yes
fix 1 clump rigid/small molecule
fix 1 clump rigid single force 1 off off on langevin 1.0 1.0 1.0 428984
style = {lj/charmm/coul/charmm} or {lj/charmm/coul/charmm/implicit} or {lj/charmm/coul/long} or {lj/charmm/coul/msm} or {lj/charmmfsw/coul/charmmfsh} or {lj/charmmfsw/coul/long}
args = list of arguments for a particular style :ul
style = {lj/cut} or {lj/cut/coul/cut} or {lj/cut/coul/debye} or {lj/cut/coul/dsf} or {lj/cut/coul/long} or {lj/cut/coul/long/cs} or {lj/cut/coul/msm} or {lj/cut/tip4p/long}
args = list of arguments for a particular style :ul
{lj/cut} args = cutoff
cutoff = global cutoff for Lennard Jones interactions (distance units)
{lj/cut/coul/cut} args = cutoff (cutoff2)
cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units)
{lj/cut/coul/debye} args = kappa cutoff (cutoff2)
kappa = inverse of the Debye length (inverse distance units)
cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units)
style = {lj/cut/soft} or {lj/cut/coul/cut/soft} or {lj/cut/coul/long/soft} or {lj/cut/tip4p/long/soft} or {lj/charmm/coul/long/soft} or {coul/cut/soft} or {coul/long/soft} or {tip4p/long/soft}
args = list of arguments for a particular style :ul
{lj/cut/soft} args = n alpha_lj cutoff
n, alpha_LJ = parameters of soft-core potential
cutoff = global cutoff for Lennard-Jones interactions (distance units)
{lj/cut/coul/cut/soft} args = n alpha_LJ alpha_C cutoff (cutoff2)
n, alpha_LJ, alpha_C = parameters of soft-core potential
cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units)
{lj/cut/coul/long/soft} args = n alpha_LJ alpha_C cutoff
n, alpha_LJ, alpha_C = parameters of the soft-core potential
cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units)
keyword = {add} or {offset} or {shift} or {extra/atom/types} or {extra/bond/types} or {extra/angle/types} or {extra/dihedral/types} or {extra/improper/types} or {extra/bond/per/atom} or {extra/angle/per/atom} or {extra/dihedral/per/atom} or {extra/improper/per/atom} or {group} or {nocoeff} or {fix} :l
{add} arg = {append} or {Nstart} or {merge}
append = add new atoms with IDs appended to current IDs
Nstart = add new atoms with IDs starting with Nstart
merge = add new atoms with their IDs unchanged
{offset} args = toff boff aoff doff ioff
toff = offset to add to atom types
boff = offset to add to bond types
aoff = offset to add to angle types
doff = offset to add to dihedral types
ioff = offset to add to improper types
{shift} args = Sx Sy Sz
Sx,Sy,Sz = distance to shift atoms when adding to system (distance units)
{extra/atom/types} arg = # of extra atom types
{extra/bond/types} arg = # of extra bond types
{extra/angle/types} arg = # of extra angle types
{extra/dihedral/types} arg = # of extra dihedral types
{extra/improper/types} arg = # of extra improper types
{extra/bond/per/atom} arg = leave space for this many new bonds per atom
{extra/angle/per/atom} arg = leave space for this many new angles per atom
{extra/dihedral/per/atom} arg = leave space for this many new dihedrals per atom
{extra/improper/per/atom} arg = leave space for this many new impropers per atom
{extra/special/per/atom} arg = leave space for extra 1-2,1-3,1-4 interactions per atom
{group} args = groupID
groupID = add atoms in data file to this group
{nocoeff} = ignore force field parameters
{fix} args = fix-ID header-string section-string
fix-ID = ID of fix to process header lines and sections of data file
header-string = header lines containing this string will be passed to fix
section-string = section names with this string will be passed to fix :pre
style = {delete} or {index} or {loop} or {world} or {universe} or {uloop} or {string} or {format} or {getenv} or {file} or {atomfile} or {python} or {internal} or {equal} or {vector} or {atom} :l
{delete} = no args
{index} args = one or more strings
{loop} args = N
N = integer size of loop, loop from 1 to N inclusive
{loop} args = N pad
N = integer size of loop, loop from 1 to N inclusive
pad = all values will be same length, e.g. 001, 002, ..., 100
{loop} args = N1 N2
N1,N2 = loop from N1 to N2 inclusive
{loop} args = N1 N2 pad
N1,N2 = loop from N1 to N2 inclusive
pad = all values will be same length, e.g. 050, 051, ..., 100
{world} args = one string for each partition of processors
{universe} args = one or more strings
{uloop} args = N
N = integer size of loop
{uloop} args = N pad
N = integer size of loop
pad = all values will be same length, e.g. 001, 002, ..., 100
{string} arg = one string
{format} args = vname fstr
vname = name of equal-style variable to evaluate
fstr = C-style format string
{getenv} arg = one string
{file} arg = filename
{atomfile} arg = filename
{python} arg = function
{internal} arg = numeric value
{equal} or {vector} or {atom} args = one formula containing numbers, thermo keywords, math operations, group functions, atom values and vectors, compute/fix/variable references
numbers = 0.0, 100, -5.4, 2.8e-4, etc
constants = PI, version, on, off, true, false, yes, no
thermo keywords = vol, ke, press, etc from "thermo_style"_thermo_style.html
Group functions are specified as keywords followed by one or two
parenthesized arguments. The first argument {ID} is the group-ID.
The {dim} argument, if it exists, is {x} or {y} or {z}. The {dir}
argument, if it exists, is {xmin}, {xmax}, {ymin}, {ymax}, {zmin}, or
{zmax}. The {dimdim} argument, if it exists, is {xx} or {yy} or {zz}
or {xy} or {yz} or {xz}.
The group function count() is the number of atoms in the group. The
group functions mass() and charge() are the total mass and charge of
the group. Xcm() and vcm() return components of the position and
velocity of the center of mass of the group. Fcm() returns a
component of the total force on the group of atoms. Bound() returns
the min/max of a particular coordinate for all atoms in the group.
Gyration() computes the radius-of-gyration of the group of atoms. See
the "compute gyration"_compute_gyration.html command for a definition
of the formula. Angmom() returns components of the angular momentum
of the group of atoms around its center of mass. Torque() returns
components of the torque on the group of atoms around its center of
mass, based on current forces on the atoms. Inertia() returns one of
6 components of the symmetric inertia tensor of the group of atoms
around its center of mass, ordered as Ixx,Iyy,Izz,Ixy,Iyz,Ixz.
Omega() returns components of the angular velocity of the group of
atoms around its center of mass.
Region functions are specified exactly the same way as group functions
except they take an extra final argument {IDR} which is the region ID.
The function is computed for all atoms that are in both the group and
the region. If the group is "all", then the only criteria for atom
inclusion is that it be in the region.
:line
Special Functions :h4
Special functions take specific kinds of arguments, meaning their
arguments cannot be formulas themselves.
The sum(x), min(x), max(x), ave(x), trap(x), and slope(x) functions
each take 1 argument which is of the form "c_ID" or "c_ID\[N\]" or
"f_ID" or "f_ID\[N\]" or "v_name". The first two are computes and the
second two are fixes; the ID in the reference should be replaced by
the ID of a compute or fix defined elsewhere in the input script. The
compute or fix must produce either a global vector or array. If it
produces a global vector, then the notation without "\[N\]" should be
used. If it produces a global array, then the notation with "\[N\]"
should be used, when N is an integer, to specify which column of the
global array is being referenced. The last form of argument "v_name"
is for a vector-style variable where "name" is replaced by the name of
the variable.
These functions operate on a global vector of inputs and reduce it to
a single scalar value. This is analogous to the operation of the
"compute reduce"_compute_reduce.html command, which performs similar
operations on per-atom and local vectors.
The sum() function calculates the sum of all the vector elements. The
min() and max() functions find the minimum and maximum element
respectively. The ave() function is the same as sum() except that it
divides the result by the length of the vector.
The trap() function is the same as sum() except the first and last
elements are multiplied by a weighting factor of 1/2 when performing
the sum. This effectively implements an integration via the
trapezoidal rule on the global vector of data. I.e. consider a set of
points, equally spaced by 1 in their x coordinate: (1,V1), (2,V2),
..., (N,VN), where the Vi are the values in the global vector of
length N. The integral from 1 to N of these points is trap(). When
appropriately normalized by the timestep size, this function is useful
for calculating integrals of time-series data, like that generated by
the "fix ave/correlate"_fix_ave_correlate.html command.
The slope() function uses linear regression to fit a line to the set
of points, equally spaced by 1 in their x coordinate: (1,V1), (2,V2),
..., (N,VN), where the Vi are the values in the global vector of
length N. The returned value is the slope of the line. If the line
has a single point or is vertical, it returns 1.0e20.
The gmask(x) function takes 1 argument which is a group ID. It
can only be used in atom-style variables. It returns a 1 for
atoms that are in the group, and a 0 for atoms that are not.
The rmask(x) function takes 1 argument which is a region ID. It can
only be used in atom-style variables. It returns a 1 for atoms that
are in the geometric region, and a 0 for atoms that are not.
The grmask(x,y) function takes 2 arguments. The first is a group ID,
and the second is a region ID. It can only be used in atom-style
variables. It returns a 1 for atoms that are in both the group and
region, and a 0 for atoms that are not in both.
The next(x) function takes 1 argument which is a variable ID (not
"v_foo", just "foo"). It must be for a file-style or atomfile-style
variable. Each time the next() function is invoked (i.e. each time
the equal-style or atom-style variable is evaluated), the following
steps occur.
For file-style variables, the current string value stored by the
file-style variable is converted to a numeric value and returned by
the function. And the next string value in the file is read and
stored. Note that if the line previously read from the file was not a
numeric string, then it will typically evaluate to 0.0, which is
likely not what you want.
For atomfile-style variables, the current per-atom values stored by
the atomfile-style variable are returned by the function. And the
next set of per-atom values in the file is read and stored.
Since file-style and atomfile-style variables read and store the first
line of the file or first set of per-atoms values when they are
defined in the input script, these are the value(s) that will be
returned the first time the next() function is invoked. If next() is
invoked more times than there are lines or sets of lines in the file,
the variable is deleted, similar to how the "next"_next.html command
operates.
:line
Feature Functions :h4
Feature functions allow to probe the running LAMMPS executable for
whether specific features are either active, defined, or available.
The functions take two arguments, a {category} and a corresponding
{argument}. The arguments are strings thus cannot be formulas
themselves (only $-style immediate variable expansion is possible).
Return value is either 1.0 or 0.0 depending on whether the function
evaluates to true or false, respectively.
The {is_active()} function allows to query for active settings which
are grouped by categories. Currently supported categories and
arguments are:
{package} (argument = {cuda} or {gpu} or {intel} or {kokkos} or {omp})
{newton} (argument = {pair} or {bond} or {any})
{pair} (argument = {single} or {respa} or {manybody} or {tail} or {shift})
{comm_style} (argument = {brick} or {tiled})
{min_style} (argument = any of the compiled in minimizer styles)
{run_style} (argument = any of the compiled in run styles)
{atom_style} (argument = any of the compiled in atom styles)
{pair_style} (argument = any of the compiled in pair styles)
{bond_style} (argument = any of the compiled in bond styles)
{angle_style} (argument = any of the compiled in angle styles)
{dihedral_style} (argument = any of the compiled in dihedral styles)
{improper_style} (argument = any of the compiled in improper styles)
{kspace_style} (argument = any of the compiled in kspace styles) :ul
Most of the settings are self-explanatory, the {single} argument in the
{pair} category allows to check whether a pair style supports a
Pair::single() function as needed by compute group/group and others
features or LAMMPS, {respa} allows to check whether the inner/middle/outer
mode of r-RESPA is supported. In the various style categories,
the checking is also done using suffix flags, if available and enabled.
Example 1: disable use of suffix for pppm when using GPU package (i.e. run it on the CPU concurrently to running the pair style on the GPU), but do use the suffix otherwise (e.g. with USER-OMP).
pair_style lj/cut/coul/long 14.0
if $(is_active(package,gpu)) then "suffix off"
kspace_style pppm :pre
Example 2: use r-RESPA with inner/outer cutoff, if supported by pair style, otherwise fall back to using pair and reducing the outer time step
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/colvars/colvars
// Please update all Colvars source files before making any changes.
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#ifndef COLVARDEPS_H
#define COLVARDEPS_H
#include "colvarmodule.h"
#include "colvarparse.h"
/// \brief Parent class for a member object of a bias, cv or cvc etc. containing features and
/// their dependencies, and handling dependency resolution
///
/// There are 3 kinds of features:
/// 1. Dynamic features are under the control of the dependency resolution
/// system. They may be enabled or disabled depending on dependencies.
/// 2. User features may be enabled based on user input (they may trigger a failure upon dependency resolution, though)
/// 3. Static features are static properties of the object, determined
/// programatically at initialization time.
///
+/// In all classes, feature 0 is active. When an object is inactivated
+/// all its children dependencies are dereferenced (free_children_deps)
+/// While the object is inactive, no dependency solving is done on children
+/// it is done when the object is activated back (restore_children_deps)
class colvardeps {
public:
- colvardeps() {}
+ colvardeps();
virtual ~colvardeps();
// Subclasses should initialize the following members:
std::string description; // reference to object name (cv, cvc etc.)
/// This contains the current state of each feature for each object
+ // since the feature class only contains static properties
struct feature_state {
feature_state(bool a, bool e)
- : available(a), enabled(e) {}
+ : available(a), enabled(e), ref_count(0) {}
/// Feature may be enabled, subject to possible dependencies
bool available;
/// Currently enabled - this flag is subject to change dynamically
/// TODO consider implications for dependency solving: anyone who disables
/// it should trigger a refresh of parent objects
bool enabled; // see if this should be private depending on implementation
+
// bool enabledOnce; // this should trigger an update when object is evaluated
+
+ /// Number of features requiring this one as a dependency
+ /// When it falls to zero:
+ /// - a dynamic feature is disabled automatically
+ /// - other features may be disabled statically
+ int ref_count;
+ /// List of features that were enabled by this one
+ /// as part of an alternate requirement (for ref counting purposes)
+ /// This is necessary because we don't know which feature in the list
+ /// we enabled, otherwise
+ std::vector<int> alternate_refs;
};
+protected:
+ /// Time step multiplier (for coarse-timestep biases & colvars)
+ /// Biases and colvars will only be calculated at those times
+ /// (f_cvb_awake and f_cv_awake); a
+ /// Biases use this to apply "impulse" biasing forces at the outer timestep
+ /// Unused by lower-level objects (cvcs and atom groups)
+ int time_step_factor;
private:
/// List of the states of all features
std::vector<feature_state> feature_states;
/// Enum of possible feature types
enum feature_type {
f_type_not_set,
f_type_dynamic,
f_type_user,
f_type_static
};
public:
+ /// \brief returns time_step_factor
+ inline int get_time_step_factor() const {return time_step_factor;}
+
/// Pair a numerical feature ID with a description and type
void init_feature(int feature_id, const char *description, feature_type type = f_type_not_set);
- /// Describes a feature and its dependecies
+ /// Describes a feature and its dependencies
/// used in a static array within each subclass
class feature {
public:
feature() {}
~feature() {}
std::string description; // Set by derived object initializer
// features that this feature requires in the same object
// NOTE: we have no safety mechanism against circular dependencies, however, they would have to be internal to an object (ie. requires_self or requires_alt)
std::vector<int> requires_self;
// Features that are incompatible, ie. required disabled
// if enabled, they will cause a dependency failure (they will not be disabled)
// To enforce these dependencies regardless of the order in which they
// are enabled, they are always set in a symmetric way, so whichever is enabled
// second will cause the dependency to fail
std::vector<int> requires_exclude;
// sets of features that are required in an alternate way
// when parent feature is enabled, if none are enabled, the first one listed that is available will be enabled
std::vector<std::vector<int> > requires_alt;
// features that this feature requires in children
std::vector<int> requires_children;
inline bool is_dynamic() { return type == f_type_dynamic; }
inline bool is_static() { return type == f_type_static; }
inline bool is_user() { return type == f_type_user; }
/// Type of this feature, from the enum feature_type
// we skip dist2(), because periodicities and the like should
// matter: boundaries should be EXACTLY the same (otherwise,
// map_grid() should be used)
if ( (std::fabs(other_grid.lower_boundaries[i] -
lower_boundaries[i]) > 1.0E-10) ||
(std::fabs(other_grid.upper_boundaries[i] -
upper_boundaries[i]) > 1.0E-10) ||
(std::fabs(other_grid.widths[i] -
widths[i]) > 1.0E-10) ||
(data.size() != other_grid.data.size()) ) {
cvm::error("Error: inconsistency between "
"two grids that are supposed to be equal, "
"aside from the data stored.\n");
return;
}
}
}
/// \brief Read grid entry in restart file
std::istream & read_restart(std::istream &is)
{
size_t const start_pos = is.tellg();
std::string key, conf;
if ((is >> key) && (key == std::string("grid_parameters"))) {
is.seekg(start_pos, std::ios::beg);
is >> colvarparse::read_block("grid_parameters", conf);
parse_params(conf, colvarparse::parse_silent);
} else {
cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n");
is.seekg(start_pos, std::ios::beg);
}
read_raw(is);
return is;
}
/// \brief Write grid entry in restart file
std::ostream & write_restart(std::ostream &os)
{
write_params(os);
write_raw(os);
return os;
}
/// \brief Write the grid data without labels, as they are
/// represented in memory
/// \param buf_size Number of values per line
std::ostream & write_raw(std::ostream &os,
size_t const buf_size = 3)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
std::vector<int> ix = new_index();
size_t count = 0;
for ( ; index_ok(ix); incr(ix)) {
for (size_t imult = 0; imult < mult; imult++) {
os << " "
<< std::setw(w) << std::setprecision(p)
<< value_output(ix, imult);
if (((++count) % buf_size) == 0)
os << "\n";
}
}
// write a final newline only if buffer is not empty
if ((count % buf_size) != 0)
os << "\n";
return os;
}
/// \brief Read data written by colvar_grid::write_raw()
std::istream & read_raw(std::istream &is)
{
size_t const start_pos = is.tellg();
for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix)) {
for (size_t imult = 0; imult < mult; imult++) {
T new_value;
if (is >> new_value) {
value_input(ix, new_value, imult);
} else {
is.clear();
is.seekg(start_pos, std::ios::beg);
is.setstate(std::ios::failbit);
cvm::error("Error: failed to read all of the grid points from file. Possible explanations: grid parameters in the configuration (lowerBoundary, upperBoundary, width) are different from those in the file, or the file is corrupt/incomplete.\n");
return is;
}
}
}
has_data = true;
return is;
}
/// \brief Write the grid in a format which is both human readable
/// and suitable for visualization e.g. with gnuplot
void write_multicol(std::ostream &os)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
// Data in the header: nColvars, then for each
// xiMin, dXi, nPoints, periodic
os << std::setw(2) << "# " << nd << "\n";
for (size_t i = 0; i < nd; i++) {
os << "# "
<< std::setw(10) << lower_boundaries[i]
<< std::setw(10) << widths[i]
<< std::setw(10) << nx[i] << " "
<< periodic[i] << "\n";
}
for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix) ) {
if (ix.back() == 0) {
// if the last index is 0, add a new line to mark the new record
os << "\n";
}
for (size_t i = 0; i < nd; i++) {
os << " "
<< std::setw(w) << std::setprecision(p)
<< bin_to_value_scalar(ix[i], i);
}
os << " ";
for (size_t imult = 0; imult < mult; imult++) {
os << " "
<< std::setw(w) << std::setprecision(p)
<< value_output(ix, imult);
}
os << "\n";
}
}
/// \brief Read a grid written by colvar_grid::write_multicol()
/// Adding data if add is true, replacing if false