2010-10-01 04:38:58 +08:00
|
|
|
<HTML>
|
|
|
|
|
<CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A>
|
|
|
|
|
</CENTER>
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
<HR>
|
|
|
|
|
|
|
|
|
|
<H3>neb command
|
|
|
|
|
</H3>
|
|
|
|
|
<P><B>Syntax:</B>
|
|
|
|
|
</P>
|
2011-01-05 07:39:13 +08:00
|
|
|
<PRE>neb etol ftol N1 N2 Nevery filename
|
2010-10-01 04:38:58 +08:00
|
|
|
</PRE>
|
2011-01-05 07:39:13 +08:00
|
|
|
<UL><LI>etol = stopping tolerance for energy (energy units)
|
|
|
|
|
<LI>ftol = stopping tolerance for force (force units)
|
2010-10-01 04:38:58 +08:00
|
|
|
<LI>N1 = max # of iterations (timesteps) to run initial NEB
|
|
|
|
|
<LI>N2 = max # of iterations (timesteps) to run barrier-climbing NEB
|
|
|
|
|
<LI>Nevery = print replica energies and reaction coordinates every this many timesteps
|
|
|
|
|
<LI>filename = file specifying final atom coordinates on other side of barrier
|
|
|
|
|
</UL>
|
|
|
|
|
<P><B>Examples:</B>
|
|
|
|
|
</P>
|
2011-01-05 07:39:13 +08:00
|
|
|
<PRE>neb 0.1 0.0 1000 500 50 coords.final
|
|
|
|
|
neb 0.0 0.001 1000 500 50 coords.final
|
2010-10-01 04:38:58 +08:00
|
|
|
</PRE>
|
|
|
|
|
<P><B>Description:</B>
|
|
|
|
|
</P>
|
|
|
|
|
<P>Perform a nudged elastic band (NEB) calculation using multiple
|
2010-10-01 06:30:55 +08:00
|
|
|
replicas of a system. Two or more replicas must be used, two of which
|
|
|
|
|
are the end points of the transition path.
|
2010-10-01 04:38:58 +08:00
|
|
|
</P>
|
|
|
|
|
<P>NEB is a method for finding both the atomic configurations and height
|
|
|
|
|
of the energy barrier associated with a transition state, e.g. for an
|
|
|
|
|
atom to perform a diffusive hop from one energy basin to another in a
|
|
|
|
|
coordinated fashion with its neighbors. The implementation in LAMMPS
|
|
|
|
|
follows the discussion in these 3 papers: <A HREF = "#Henkelman1">(Henkelman1)</A>,
|
|
|
|
|
<A HREF = "#Henkelman2">(Henkelman2)</A>, and <A HREF = "#Nakano">(Nakano)</A>.
|
|
|
|
|
</P>
|
|
|
|
|
<P>Each replica runs on a partition of one or more processors. Processor
|
|
|
|
|
partitions are defined at run-time using the -partition command-line
|
|
|
|
|
switch; see <A HREF = "Section_start.html#2_6">this section</A> of the manual. Note
|
|
|
|
|
that if you have MPI installed, you can run a multi-replica simulation
|
|
|
|
|
with more replicas (partitions) than you have physical processors, e.g
|
|
|
|
|
you can run a 10-replica simulation on one or two processors. You
|
|
|
|
|
will simply not get the performance speed-up you would see with one or
|
|
|
|
|
more physical processors per replica. See <A HREF = "Section_howto.html#4_5">this
|
|
|
|
|
section</A> of the manual for further discussion.
|
|
|
|
|
</P>
|
|
|
|
|
<P>NOTE: The current NEB implementation in LAMMPS restricts you to having
|
|
|
|
|
exactly one processor per replica.
|
|
|
|
|
</P>
|
|
|
|
|
<P>When a NEB calculation is performed, it is assumed that each replica
|
|
|
|
|
is running the same model, though LAMMPS does not check for this.
|
|
|
|
|
I.e. the simulation domain, the number of atoms, the interaction
|
|
|
|
|
potentials, and the starting configuration when the neb command is
|
|
|
|
|
issued should be the same for every replica.
|
|
|
|
|
</P>
|
|
|
|
|
<P>In a NEB calculation each atom in a replica is connected to the same
|
|
|
|
|
atom in adjacent replicas by springs, which induce inter-replica
|
|
|
|
|
forces. These forces are imposed by the <A HREF = "fix_neb.html">fix neb</A>
|
|
|
|
|
command, which must be used in conjunction with the neb command. The
|
|
|
|
|
group used to define the fix neb command specifies which atoms the
|
|
|
|
|
inter-replica springs are applied to. These are the NEB atoms.
|
|
|
|
|
Additional atoms can be present in your system, e.g. to provide a
|
|
|
|
|
background force field or simply to hold fixed during the NEB
|
|
|
|
|
procedure, but they will not be part of the barrier finding procedure.
|
|
|
|
|
</P>
|
|
|
|
|
<P>The "starting configuration" for NEB should be a state with the NEB
|
|
|
|
|
atoms (and all other atoms) having coordinates on one side of the
|
|
|
|
|
energy barrier. These coordinates will be assigned to the first
|
2010-10-01 06:30:55 +08:00
|
|
|
replica #1. The coordinates should be close to a local energy
|
|
|
|
|
minimum. A perfect energy minimum is not required, since NEB runs via
|
|
|
|
|
damped dynamics which will tend to drive the configuration of replica
|
|
|
|
|
#1 to a true energy minimum, but you will typically get better
|
|
|
|
|
convergence if the initial state is already at a minimum. For
|
|
|
|
|
example, for a system with a free surface, the surface should be fully
|
|
|
|
|
relaxed before attempting a NEB calculation.
|
2010-10-01 04:38:58 +08:00
|
|
|
</P>
|
|
|
|
|
<P>The final configuration is specified in the input <I>filename</I>, which is
|
|
|
|
|
formatted as described below. Only coordinates for NEB atoms or a
|
|
|
|
|
subset of them should be listed in the file; they represent the state
|
|
|
|
|
of the system on the other side of the barrier, at or near an energy
|
|
|
|
|
minimum. These coordinates will be assigned to the last replica #M.
|
|
|
|
|
The final coordinates of atoms not listed in <I>filename</I> are set equal
|
|
|
|
|
to their initial coordinates. Again, a perfect energy minimum is not
|
|
|
|
|
required for the final configuration, since the atoms in replica #M
|
|
|
|
|
will tend to move during the NEB procedure to the nearest energy
|
|
|
|
|
minimum. Also note that a final coordinate does not need to be
|
|
|
|
|
specified for a NEB atom if you expect it to only displace slightly
|
|
|
|
|
during the NEB procedure. For example, only the final coordinate of
|
|
|
|
|
the single atom diffusing into a vacancy need be specified if the
|
|
|
|
|
surrounding atoms will only relax slightly in the final configuration.
|
|
|
|
|
</P>
|
|
|
|
|
<P>The initial coordinates of all atoms (not just NEB atoms) in the
|
|
|
|
|
intermediate replicas #2,#3,...,#M-1 are set to values linearly
|
|
|
|
|
interpolated between the corresponding atoms in replicas #1 and #M.
|
|
|
|
|
</P>
|
|
|
|
|
<P>A NEB calculation has two stages, each of which is a minimization
|
|
|
|
|
procedure, performed via damped dynamics. To enable this, you must
|
|
|
|
|
first define an appropriate <A HREF = "min_style.html">min_style</A>, such as
|
|
|
|
|
<I>quickmin</I> or <I>fire</I>. The <I>cg</I>, <I>sd</I>, and <I>hftn</I> styles cannot be
|
|
|
|
|
used, since they perform iterative line searches in their inner loop,
|
|
|
|
|
which cannot be easily synchronized across multiple replicas.
|
|
|
|
|
</P>
|
2011-01-05 07:39:13 +08:00
|
|
|
<P>The minimizer tolerances for energy and force are set by <I>etol</I> and <I>ftol</I>,
|
|
|
|
|
the same as for
|
|
|
|
|
the <A HREF = "minimize.html">minimize</A> command.
|
|
|
|
|
</P>
|
|
|
|
|
<P>A non-zero <I>etol</I>
|
|
|
|
|
means that the NEB calculation will terminate if the energy criterion is met
|
|
|
|
|
by every replica. The energies being compared to
|
|
|
|
|
<I>etol</I> do not include any contribution from the inter-replica forces, since
|
|
|
|
|
these are non-conservative.
|
|
|
|
|
A non-zero <I>ftol</I>
|
|
|
|
|
means that the NEB calculation will terminate if the force criterion is met
|
|
|
|
|
by every replica. The forces being compared to
|
2010-10-05 07:31:54 +08:00
|
|
|
<I>ftol</I> include the inter-replica forces between an atom and its images
|
|
|
|
|
in adjacent replicas.
|
|
|
|
|
</P>
|
|
|
|
|
<P>The maximum number of iterations in each stage is set by <I>N1</I> and
|
|
|
|
|
<I>N2</I>. These are effectively timestep counts since each iteration of
|
|
|
|
|
damped dynamics is like a single timestep in a dynamics
|
|
|
|
|
<A HREF = "run.html">run</A>. During both stages, the potential energy of each
|
|
|
|
|
replica and its normalized distance along the reaction path (reaction
|
|
|
|
|
coordinate RD) will be printed to the screen and log file every
|
|
|
|
|
<I>Nevery</I> timesteps. The RD is 0 and 1 for the first and last replica.
|
2011-01-05 07:39:13 +08:00
|
|
|
For intermediate replicas, it is the cumulative distance (normalized
|
|
|
|
|
by the total cumulative distance) between adjacent replicas, where
|
2010-10-05 07:31:54 +08:00
|
|
|
"distance" is defined as the length of the 3N-vector of differences in
|
|
|
|
|
atomic coordinates, where N is the number of NEB atoms involved in the
|
2010-10-01 04:38:58 +08:00
|
|
|
transition. These outputs allow you to monitor NEB's progress in
|
|
|
|
|
finding a good energy barrier. <I>N1</I> and <I>N2</I> must both be multiples
|
|
|
|
|
of <I>Nevery</I>.
|
|
|
|
|
</P>
|
|
|
|
|
<P>In the first stage of NEB, the set of replicas should converge toward
|
|
|
|
|
the minimum energy path (MEP) of conformational states that transition
|
|
|
|
|
over the barrier. The MEP for a barrier is defined as a sequence of
|
|
|
|
|
3N-dimensional states that cross the barrier at its saddle point, each
|
|
|
|
|
of which has a potential energy gradient parallel to the MEP itself.
|
|
|
|
|
The replica states will also be roughly equally spaced along the MEP
|
|
|
|
|
due to the inter-replica spring force added by the <A HREF = "fix_neb.html">fix
|
|
|
|
|
neb</A> command.
|
|
|
|
|
</P>
|
2010-11-11 01:29:32 +08:00
|
|
|
<P>In the second stage of NEB, the replica with the highest energy
|
2010-10-01 04:38:58 +08:00
|
|
|
is selected and the inter-replica forces on it are converted to a
|
|
|
|
|
force that drives its atom coordinates to the top or saddle point of
|
2010-11-11 01:29:32 +08:00
|
|
|
the barrier, via the barrier-climbing calculation described in
|
2010-10-01 04:38:58 +08:00
|
|
|
<A HREF = "#Hinkelman2">(Henkelman2)</A>. As before, the other replicas rearrange
|
|
|
|
|
themselves along the MEP so as to be roughly equally spaced.
|
|
|
|
|
</P>
|
|
|
|
|
<P>When both stages are complete, if the NEB calculation was successful,
|
|
|
|
|
one of the replicas should be an atomic configuration at the top or
|
|
|
|
|
saddle point of the barrier, the potential energies for the set of
|
|
|
|
|
replicas should represent the energy profile of the barrier along the
|
2011-01-05 07:39:13 +08:00
|
|
|
MEP, and the configurations of the replicas should be a sequence of
|
2010-10-01 04:38:58 +08:00
|
|
|
configurations along the MEP.
|
|
|
|
|
</P>
|
|
|
|
|
<HR>
|
|
|
|
|
|
|
|
|
|
<P>A few other settings in your input script are required or advised to
|
|
|
|
|
perform a NEB calculation.
|
|
|
|
|
</P>
|
|
|
|
|
<P>The "atom_modify sort 0 0.0" command should be used to turn off atom
|
|
|
|
|
sorting.
|
|
|
|
|
</P>
|
|
|
|
|
<P>NOTE: This sorting restriction will be removed in a future version of
|
|
|
|
|
NEB in LAMMPS.
|
|
|
|
|
</P>
|
|
|
|
|
<P>The minimizers in LAMMPS operate on all atoms in your system, even
|
|
|
|
|
non-NEB atoms, as defined above. To prevent non-NEB atoms from moving
|
|
|
|
|
during the minimization, you should use the <A HREF = "fix_setforce.html">fix
|
|
|
|
|
setforce</A> command to set the force on each of those
|
|
|
|
|
atoms to 0.0. This is not required, and may not even be desired in
|
|
|
|
|
some cases, but if those atoms move too far (e.g. because the initial
|
|
|
|
|
state of your system was not well-minimized), it can cause problems
|
|
|
|
|
for the NEB procedure.
|
|
|
|
|
</P>
|
|
|
|
|
<P>The damped dynamics <A HREF = "min_style.html">minimizers</A>, such as <I>quickmin</I>
|
|
|
|
|
and <I>fire</I>), adjust the position and velocity of the atoms via an
|
|
|
|
|
Euler integration step. Thus you must define an appropriate
|
|
|
|
|
<A HREF = "timestep.html">timestep</A> to use with NEB. Using the same timestep
|
|
|
|
|
that would be used for a dynamics <A HREF = "run.html">run</A> of your system is
|
|
|
|
|
advised.
|
|
|
|
|
</P>
|
|
|
|
|
<HR>
|
|
|
|
|
|
|
|
|
|
<P>The specified <I>filename</I> contains atom coordinates for the final
|
|
|
|
|
configuration. Only atoms with coordinates different than the initial
|
|
|
|
|
configuration need to be specified, i.e. those geometrically near the
|
|
|
|
|
barrier.
|
|
|
|
|
</P>
|
|
|
|
|
<P>The file can be ASCII text or a gzipped text file (detected by a .gz
|
|
|
|
|
suffix). The file should contain one line per atom, formatted
|
|
|
|
|
with the atom ID, followed by the final x,y,z coordinates:
|
|
|
|
|
</P>
|
|
|
|
|
<PRE>125 24.97311 1.69005 23.46956
|
|
|
|
|
126 1.94691 2.79640 1.92799
|
|
|
|
|
127 0.15906 3.46099 0.79121
|
|
|
|
|
...
|
|
|
|
|
</PRE>
|
|
|
|
|
<P>The lines can be listed in any order.
|
|
|
|
|
</P>
|
|
|
|
|
<HR>
|
|
|
|
|
|
|
|
|
|
<P>Four kinds of output can be generated during a NEB calculation: energy
|
|
|
|
|
barrier statistics, thermodynamic output by each replica, dump files,
|
|
|
|
|
and restart files.
|
|
|
|
|
</P>
|
|
|
|
|
<P>When running with multiple partitions (each of which is a replica in
|
2010-10-07 03:43:46 +08:00
|
|
|
this case), the print-out to the screen and master log.lammps file
|
|
|
|
|
contains a line of output, printed once every <I>Nevery</I> timesteps. It
|
|
|
|
|
contains the timestep, the maximum force per replica, the maximum
|
2010-11-11 01:29:32 +08:00
|
|
|
force per atom (in any replica), potential gradients in the initial,
|
2011-01-05 07:39:13 +08:00
|
|
|
final, and climbing replicas,
|
|
|
|
|
the forward and backward energy barriers,
|
|
|
|
|
the total reaction coordinate (RDT), and
|
|
|
|
|
the normalized reaction coordinate and potential energy of each replica.
|
|
|
|
|
</P>
|
|
|
|
|
<P>The "maximum force per replica" is
|
2010-10-07 03:43:46 +08:00
|
|
|
the two-norm of the 3N-length force vector for the atoms in each
|
|
|
|
|
replica, maximized across replicas, which is what the <I>ftol</I> setting
|
|
|
|
|
is checking against. In this case, N is all the atoms in each
|
|
|
|
|
replica. The "maximum force per atom" is the maximum force component
|
2010-11-11 01:29:32 +08:00
|
|
|
of any atom in any replica. The potential gradients are the two-norm
|
|
|
|
|
of the 3N-length force vector solely due to the interaction potential i.e.
|
|
|
|
|
without adding in inter-replica forces. Note that inter-replica forces
|
|
|
|
|
are zero in the initial and final replicas, and only affect
|
|
|
|
|
the direction in the climbing replica. For this reason, the "maximum
|
|
|
|
|
force per replica" is often equal to the potential gradient in the
|
|
|
|
|
climbing replica. In the first stage of NEB, there is no climbing
|
|
|
|
|
replica, and so the potential gradient in the highest energy replica
|
|
|
|
|
is reported, since this replica will become the climbing replica
|
|
|
|
|
in the second stage of NEB.
|
|
|
|
|
</P>
|
2011-01-05 07:39:13 +08:00
|
|
|
<P>The "reaction coordinate" (RD) for each
|
2010-11-11 01:29:32 +08:00
|
|
|
replica is the two-norm of the 3N-length vector of distances between
|
2011-01-05 07:39:13 +08:00
|
|
|
its atoms and the preceding replica's atoms, added to the RD of the
|
|
|
|
|
preceding replica. The RD of the first replica RD1 = 0.0;
|
|
|
|
|
the RD of the final replica RDN = RDT, the total reaction coordinate.
|
|
|
|
|
The normalized RDs are divided by RDT,
|
|
|
|
|
so that they form a monotonically increasing sequence
|
|
|
|
|
from zero to one. When computing RD, N only includes the atoms
|
|
|
|
|
being operated on by the fix neb command.
|
|
|
|
|
</P>
|
|
|
|
|
<P>The forward (reverse) energy barrier is the potential energy of the highest
|
|
|
|
|
replica minus the energy of the first (last) replica.
|
2010-10-01 04:38:58 +08:00
|
|
|
</P>
|
|
|
|
|
<P>When running on multiple partitions, LAMMPS produces additional log
|
|
|
|
|
files for each partition, e.g. log.lammps.0, log.lammps.1, etc. For a
|
|
|
|
|
NEB calculation, these contain the thermodynamic output for each
|
|
|
|
|
replica.
|
|
|
|
|
</P>
|
|
|
|
|
<P>If <A HREF = "dump.html">dump</A> commands in the input script define a filename
|
|
|
|
|
that includes a <I>universe</I> or <I>uloop</I> style <A HREF = "variable.html">variable</A>,
|
|
|
|
|
then one dump file (per dump command) will be created for each
|
|
|
|
|
replica. At the end of the NEB calculation, the final snapshot in
|
|
|
|
|
each file will contain the sequence of snapshots that transition the
|
|
|
|
|
system over the energy barrier. Earlier snapshots will show the
|
|
|
|
|
convergence of the replicas to the MEP.
|
|
|
|
|
</P>
|
|
|
|
|
<P>Likewise, <A HREF = "restart.html">restart</A> filenames can be specified with a
|
|
|
|
|
<I>universe</I> or <I>uloop</I> style <A HREF = "variable.html">variable</A>, to generate
|
|
|
|
|
restart files for each replica. These may be useful if the NEB
|
|
|
|
|
calculation fails to converge properly to the MEP, and you wish to
|
|
|
|
|
restart the calculation from an intermediate point with altered
|
|
|
|
|
parameters.
|
|
|
|
|
</P>
|
|
|
|
|
<P>There are 2 Python scripts provided in the tools/python directory,
|
2010-11-05 22:29:46 +08:00
|
|
|
neb_combine.py and neb_final.py, which are useful in analyzing output
|
|
|
|
|
from a NEB calculation. Assume a NEB simulation with M replicas, and
|
|
|
|
|
the NEB atoms labelled with a specific atom type.
|
|
|
|
|
</P>
|
|
|
|
|
<P>The neb_combine.py script extracts atom coords for the NEB atoms from
|
|
|
|
|
all M dump files and creates a single dump file where each snapshot
|
|
|
|
|
contains the NEB atoms from all the replicas and one copy of non-NEB
|
|
|
|
|
atoms from the first replica (presumed to be identical in other
|
|
|
|
|
replicas). This can be visualized/animated to see how the NEB atoms
|
|
|
|
|
relax as the NEB calculation proceeds.
|
|
|
|
|
</P>
|
|
|
|
|
<P>The neb_final.py script extracts the final snapshot from each of the M
|
|
|
|
|
dump files to create a single dump file with M snapshots. This can be
|
2010-10-01 04:38:58 +08:00
|
|
|
visualized to watch the system make its transition over the energy
|
|
|
|
|
barrier.
|
|
|
|
|
</P>
|
2010-11-05 22:29:46 +08:00
|
|
|
<P>To illustrate, here are images from the final snapshot produced by the
|
|
|
|
|
neb_combine.py script run on the dump files produced by the two
|
|
|
|
|
example input scripts in examples/neb. Click on them to see a larger
|
|
|
|
|
image.
|
|
|
|
|
</P>
|
|
|
|
|
<A HREF = "JPG/hop1.jpg"><IMG SRC = "JPG/hop1_small.jpg"></A>
|
|
|
|
|
|
|
|
|
|
<A HREF = "JPG/hop2.jpg"><IMG SRC = "JPG/hop2_small.jpg"></A>
|
|
|
|
|
|
2010-10-01 04:38:58 +08:00
|
|
|
<HR>
|
|
|
|
|
|
|
|
|
|
<P><B>Restrictions:</B>
|
|
|
|
|
</P>
|
|
|
|
|
<P>This command can only be used if LAMMPS was built with the "replica"
|
|
|
|
|
package. See the <A HREF = "Section_start.html#2_3">Making LAMMPS</A> section for
|
|
|
|
|
more info on packages.
|
|
|
|
|
</P>
|
|
|
|
|
<P><B>Related commands:</B>
|
|
|
|
|
</P>
|
|
|
|
|
<P><A HREF = "prd.html">prd</A>, <A HREF = "temper.html">temper</A>, <A HREF = "fix_langevin.html">fix
|
|
|
|
|
langevin</A>, <A HREF = "fix_viscous.html">fix viscous</A>
|
|
|
|
|
</P>
|
|
|
|
|
<P><B>Default:</B> none
|
|
|
|
|
</P>
|
|
|
|
|
<HR>
|
|
|
|
|
|
|
|
|
|
<A NAME = "Henkelman1"></A>
|
|
|
|
|
|
2010-10-01 05:22:06 +08:00
|
|
|
<P><B>(Henkelman1)</B> Henkelman and Jonsson, J Chem Phys, 113, 9978-9985 (2000).
|
2010-10-01 04:38:58 +08:00
|
|
|
</P>
|
|
|
|
|
<A NAME = "Henkelman2"></A>
|
|
|
|
|
|
2010-10-01 05:22:06 +08:00
|
|
|
<P><B>(Henkelman2)</B> Henkelman, Uberuaga, Jonsson, J Chem Phys, 113,
|
2010-10-01 04:38:58 +08:00
|
|
|
9901-9904 (2000).
|
|
|
|
|
</P>
|
|
|
|
|
<A NAME = "Nakano"></A>
|
|
|
|
|
|
|
|
|
|
<P><B>(Nakano)</B> Nakano, Comp Phys Comm, 178, 280-289 (2008).
|
|
|
|
|
</P>
|
|
|
|
|
</HTML>
|
