added several features to the NEB

# SOME DESCRIPTIVE TITLE.

# Copyright (C) YEAR Free Software Foundation, Inc.

# FIRST AUTHOR <EMAIL@ADDRESS>, YEAR.

#

#, fuzzy

msgid ""

msgstr ""

"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c

:link(lws,http://lammps.sandia.gov)

[Syntax:]

fix ID group-ID neb Kspring :pre

fix ID group-ID neb Kspring keyword value :pre

ID, group-ID are documented in "fix"_fix.html command

neb = style name of this fix command

Kspring = inter-replica spring constant (force/distance units) :ul

Kspring = parallel spring constant (force/distance units) :ul

keyword = {idealpos} or {nearestneigh} or {PerpSpring} or {freeend}

  {idealpos} value = none = each replica is attached with a spring to its interpolated ideal position (default value)

  {nearestneigh} value = none = each replica is attached with a spring with the previous and next replica.

  {PerpSpring} value = KspringPerpend

  {freeend} value = ini or final or finalWithRespToIni or finalAndInterWithRespToIni

[Examples:]

fix 1 active neb 10.0 :pre

fix 2 all neb 1.0 PerpSpring 1.0 freeend final :pre

fix 1 all neb 1.0 nearestneigh freeend finalAndInterWithRespToIni :pre

[Description:]

Add inter-replica forces to atoms in the group for a multi-replica

simulation run via the "neb"_neb.html command to perform a nudged

elastic band (NEB) calculation for transition state finding.  Hi-level

explanations of NEB are given with the "neb"_neb.html command and in

"Section 6.5"_Section_howto.html#howto_5 of the manual.  The fix

neb command must be used with the "neb" command to define how

inter-replica forces are computed.

Only the N atoms in the fix group experience inter-replica forces.

Atoms in the two end-point replicas do not experience these forces,

but those in intermediate replicas do.  During the initial stage of

NEB, the 3N-length vector of interatomic forces Fi = -Grad(V) acting

on the atoms of each intermediate replica I is altered, as described

in the "(Henkelman1)"_#Henkelman1 paper, to become:

Fi = -Grad(V) + (Grad(V) dot That) That + Kspring (| Ri+i - Ri | - | Ri - Ri-1 |) That :pre

Ri are the atomic coordinates of replica I; Ri-1 and Ri+1 are the

coordinates of its neighbor replicas.  That (t with a hat over it) is

the unit "tangent" vector for replica I which is a function of Ri,

Ri-1, Ri+1, and the potential energy of the 3 replicas; it points

roughly in the direction of (Ri+i - Ri-1); see the

"(Henkelman1)"_#Henkelman1 paper for details.

The first two terms in the above equation are the component of the

interatomic forces perpendicular to the tangent vector.  The last term

is a spring force between replica I and its neighbors, parallel to the

tangent vector direction with the specified spring constant {Kspring}.

The effect of the first two terms is to push the atoms of each replica

toward the minimum energy path (MEP) of conformational states that

transition over the energy barrier.  The MEP for an energy barrier is

defined as a sequence of 3N-dimensional states which cross the barrier

at its saddle point, each of which has a potential energy gradient

parallel to the MEP itself.

The effect of the last term is to push each replica away from its two

neighbors in a direction along the MEP, so that the final set of

states are equidistant from each other.

During the second stage of NEB, the forces on the N atoms in the

replica nearest the top of the energy barrier are altered so that it

climbs to the top of the barrier and finds the saddle point.  The

forces on atoms in this replica are described in the

"(Henkelman2)"_#Henkelman2 paper, and become:

Add a nudging force to atoms in the group for a multi-replica simulation run via the "neb"_neb.html command to perform a nudged elastic band (NEB) calculation for finding the transition state.  Hi-level

explanations of NEB are given with the "neb"_neb.html command and in "Section_howto 5"_Section_howto.html#howto_5 of the manual.  The fix neb command must be used with the "neb" command and defines how

nudging inter-replica forces are computed.

A NEB calculation is divided in two stages. In the first stage n replicas are relaxed toward a MEP and in a second stage, the climbing image scheme (see "(Henkelman2)"_#Henkelman2) is turned on so that the replica having the highest energy relaxes toward the saddle point (i.e. the point of highest energy along the MEP).

One purpose of the nudging forces is to keep the replicas equally spaced. 

During the NEB, the 3N-length vector of interatomic force Fi = -Grad(V) of replicas i is altered. For all intermediate replicas (i.e. for 1<i<n) except for the climbing replica the force vector becomes:

Fi = -Grad(V) + (Grad(V) dot That) That + Fspringparallel + Fspringperp :pre

That is the unit "tangent" vector for replica i and is a function of Ri, Ri-1,

Ri+1, and the potential energy of the 3 replicas; it points roughly in the direction of (Ri+i - Ri-1) (see the

"(Henkelman1)"_#Henkelman1 paper for details).

Ri are the atomic coordinates of replica i; Ri-1 and Ri+1 are the

coordinates of its neighbor replicas.  

The term (Grad(V) dot That) is used to remove the component of the gradient parallel to the path which would tend to distribute the replica unevenly along the path.

Fspringparallel is an artificial spring force which is applied only in the tangent direction and which maintains the replicas equally spaced (see below for more information). 

Fspringperp is an optinal artificial spring which is applied only perpendicular to the tangent and which prevent the paths from forming too acute kinks (see below for more information). 

In the second stage of the NEB, the interatomic force Fi for the climbing replica (which is the replica of highest energy) becomes :

Fi = -Grad(V) + 2 (Grad(V) dot That) That :pre

The inter-replica forces for the other replicas are unchanged from the

first equation.

By default, the force acting on the first and last replicas is not altered so that during the NEB relaxation, these ending replicas relax toward local minima. However it is possible to use the key word {freeend} to allow either the initial or the final replica to relax toward a MEP while constraining its energy. 

The interatomic force Fi for the free end image becomes : 

Fi = -Grad(V)+ (Grad(V) dot That + E-ETarget) That when Grad(V) dot That < 0

Fi = -Grad(V)+ (Grad(V) dot That + ETarget- E) That when Grad(V) dot That > 0

where E is the energy of the free end replica and ETarget is the target energy.

When the value {ini} ({final}) is used after the keyword {freeend}, the first (last) replica is considered as a free end. The target energy is set to the energy of the replica at starting of the NEB calculation. When the value {finalWithRespToIni} or {finalAndInterWithRespToIni} is used the last image is considered as a free end and the target energy is equal to the energy of the first replica (which can evolve during the NEB relaxation). 

With the value {finalWithRespToIni}, when the initial path is too far from the MEP, an intermediate repilica might relax "faster" and get a lower energy than the last replica. The benefit of the free end is then lost since this intermediate replica will relax toward a local minima. This behavior can be prevented by using the value {finalAndInterWithRespToIni} which remove entirely the contribution of the gradient for all intermediate replica which have a lower energy than the initial one thus preventing these replicae to over-relax.  After converging a NEB with the {finalAndInterWithRespToIni} value it

is recommended to check that all intermediate replica have a larger energy than the initial replica. Finally note that if the last replica converges toward a local minimum with a larger energy than the energy of the first replica, a free end neb calculation with the value {finalWithRespToIni} or {finalAndInterWithRespToIni} cannot reach the convergence criteria.

:line

The keywords {idealpos} and {nearestneigh} allow to specify how to parallel spring force is computed. 

If the keyword {idealpos} is used or by default, the spring force is computed as suggested in "(E)"_#E :

Fspringparallel=-{Kspring}* (RD-RDideal)/(2 meanDist)

where RD is the "reaction coordinate" see "neb"_neb.html section, and RDideal is the ideal RD for which all the images are equally spaced (i.e. RDideal = (i-1)*meanDist when the climbing image is off, where i is the replica number). The meanDist is the average distance between replicas.

If the keyword  {nearestneigh} is used, the parallel spring force is computed as in "(Henkelman1)"_#Henkelman1 by connecting each intermediate replica with the previous and the next image:

Fspringparallel= {Kspring}* (|Ri+1 - Ri| - |Ri - Ri-1|)

The parallel spring force associated with the key word idealpos should usually be more efficient at keeping the images equally spaced.

:line

The keyword {PerpSpring} allows to add a spring force perpendicular to the path in order to prevent the path from becoming too kinky. It can improve significantly the convergence of the NEB when the resolution is poor (i.e. when too few images are used) (see "(Maras)"_#Maras).

The perpendicular spring force is given by 

Fspringperp = {Kspringperp} * f(Ri-1,Ri,Ri+1) (Ri+1 + Ri-1 - 2 Ri)

 f(Ri-1 Ri R+1) is a smooth scalar function of the angle Ri-1 Ri Ri+1. It is equal to 0 when the path is straight and is equal to 1 when the angle Ri-1 Ri Ri+1 is accute. f(Ri-1 Ri R+1) is defined in "(Jonsson)"_#Jonsson 

:line

[Restart, fix_modify, output, run start/stop, minimize info:]

[Default:] none

:link(Henkelman1)

:link(Henkelman)

[(Henkelman1)] Henkelman and Jonsson, J Chem Phys, 113, 9978-9985 (2000).

:link(Henkelman2)

:link(Henkelman)

[(Henkelman2)] Henkelman, Uberuaga, Jonsson, J Chem Phys, 113,

9901-9904 (2000).

:link(E)

[(E)] E, Ren, Vanden-Eijnden, Phys Rev B, 66, 052301 (2002)

:link(Jonsson)

[(Jonsson)] Jonsson, Mills and  Jacobsen, in Classical and Quantum Dynamics in Condensed Phase Simulations, edited by  Berne, Ciccotti, and  Coker ͑World Scientific, Singapore, 1998͒, p. 385

:link(Maras)

[(Maras)] Maras, Trushin, Stukowski, Ala-Nissila, Jonsson, Comp Phys Comm, 205, 13-21 (2016)

[Syntax:]

neb etol ftol N1 N2 Nevery file-style arg :pre

neb etol ftol N1 N2 Nevery file-style arg keyword :pre

etol = stopping tolerance for energy (energy units) :ulb,l

ftol = stopping tolerance for force (force units) :l

    filename = unique filename for each replica (except first) with its initial coords

  {none} arg = no argument

    all replicas assumed to already have their initial coords :pre

keyword = {verbose} :pre

:ule

[Examples:]

neb 0.1 0.0 1000 500 50 final coords.final

neb 0.0 0.001 1000 500 50 each coords.initial.$i

neb 0.0 0.001 1000 500 50 none :pre

neb 0.0 0.001 1000 500 50 none verbose :pre

[Description:]

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: "(HenkelmanA)"_#HenkelmanA,

"(HenkelmanB)"_#HenkelmanB, and "(Nakano)"_#Nakano3.

follows the discussion in these 4 papers: "(HenkelmanA)"_#HenkelmanA,

"(HenkelmanB)"_#HenkelmanB, "(Nakano)"_#Nakano and "(Maras)"_#Maras.

Each replica runs on a partition of one or more processors.  Processor

partitions are defined at run-time using the -partition command-line

discussion.

NOTE: As explained below, a NEB calculation perfoms a damped dynamics

minimization across all the replicas.  The minimizer uses whatever

minimization across all the replicas.  The mimimizer uses whatever

timestep you have defined in your input script, via the

"timestep"_timestep.html command.  Often NEB will converge more

quickly if you use a timestep about 10x larger than you would normally

potentials, and the starting configuration when the neb command is

issued should be the same for every replica.

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 "fix neb"_fix_neb.html

In a NEB calculation each replica is connected to other replicas by inter-replica

nudging forces.  These forces are imposed by the "fix neb"_fix_neb.html

command, which must be used in conjunction with the neb command.  The

group used to define the fix neb command defines the NEB atoms which

are the only ones that inter-replica springs are applied to.  If the

computed in the usual way due only to other atoms within their

replica.  Conceptually, the non-NEB atoms provide a background force

field for the NEB atoms.  They can be allowed to move during the NEB

minimization procedure (which will typically induce different

minimiation procedure (which will typically induce different

coordinates for non-NEB atoms in different replicas), or held fixed

using other LAMMPS commands such as "fix setforce"_fix_setforce.html.

Note that the "partition"_partition.html command can be used to invoke

discussed below.  Only atoms whose initial coordinates should differ

from the current configuration need be specified.

Conceptually, the initial configuration for the first replica should

be a state with all the atoms (NEB and non-NEB) having coordinates on

one side of the energy barrier.  A perfect energy minimum is not

required, since atoms in the first replica experience no spring forces

from the 2nd replica.  Thus the damped dynamics minimization will

drive the first replica to an energy minimum if it is not already

there.  However, 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.

Likewise, the initial configuration of the final replica should be a

state with all the atoms (NEB and non-NEB) on the other side of the

energy barrier.  Again, a perfect energy minimum is not required,

since the atoms in the last replica also experience no spring forces

from the next-to-last replica, and thus the damped dynamics

minimization will drive it to an energy minimum.

Conceptually, the initial and final configurations for the first replica should

be states on either side of an energy barrier. 

As explained below, the initial configurations of intermediate

replicas can be atomic coordinates interpolated in a linear fashion

between the first and last replicas.  This is often adequate state for

between the first and last replicas.  This is often adequate for

simple transitions.  For more complex transitions, it may lead to slow

convergence or even bad results if the minimum energy path (MEP, see

below) of states over the barrier cannot be correctly converged to

from such an initial configuration.  In this case, you will want to

from such an initial path.  In this case, you will want to

generate initial states for the intermediate replicas that are

geometrically closer to the MEP and read them in.

replicas, the 2nd replica will assign a position that is 10% of the

distance along a line between the starting and final point, and the

9th replica will assign a position that is 90% of the distance along

the line.  Note that this procedure to produce consistent coordinates

the line.  Note that for this procedure to produce consistent coordinates

across all the replicas, the current coordinates need to be the same

in all replicas.  LAMMPS does not check for this, but invalid initial

configurations will likely result if it is not the case.

Each replica (except the first replica) will read its file, formatted

as described below, and for any atom that appears in the file, assign

the specified coordinates to its atom.  The various files do not need

the specified coordinates to this atom.  The various files do not need

to contain the same set of atoms.

For a {file-style} setting of {none}, no filename is specified.  Each

A non-zero {etol} means that the NEB calculation will terminate if the

energy criterion is met by every replica.  The energies being compared

to {etol} do not include any contribution from the inter-replica

to {etol} do not include any contribution from the inter-replica nudging

forces, since these are non-conservative.  A non-zero {ftol} means

that the NEB calculation will terminate if the force criterion is met

by every replica.  The forces being compared to {ftol} include the

inter-replica forces between an atom and its images in adjacent

replicas.

inter-replica nudging forces.

The maximum number of iterations in each stage is set by {N1} and

{N2}.  These are effectively timestep counts since each iteration of

of {Nevery}.

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.

a minimum energy path (MEP) of conformational states that transition

over a barrier.  The MEP for a transition is defined as a sequence of

3N-dimensional states, each  of which has a potential energy gradient parallel to the MEP itself.

The configuration of highest energy along a MEP corresponds to a saddle point.

The replica states will also be roughly equally spaced along the MEP

due to the inter-replica spring force added by the "fix

due to the inter-replica nugding force added by the "fix

neb"_fix_neb.html command.

In the second stage of NEB, the replica with the highest energy

"(HenkelmanB)"_#HenkelmanB.  As before, the other replicas rearrange

themselves along the MEP so as to be roughly equally spaced.

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

MEP, and the configurations of the replicas should be a sequence of

configurations along the MEP.

When both stages are complete, if the NEB calculation was successful, the configurations of the replicas should be along (close to) the MEP and the replica with the highest energy should be an atomic configuration at (close to) the saddle point of

the transition. The potential energies for the set of

replicas represents the energy profile of the transition along the

MEP.

:line

replica.  The "maximum force per atom" is the maximum force component

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.

without adding in inter-replica forces.    

The "reaction coordinate" (RD) for each

replica is the two-norm of the 3N-length vector of distances between

The forward (reverse) energy barrier is the potential energy of the highest

replica minus the energy of the first (last) replica.

Supplementary informations for all replicas  can be printed out to the screen and master

log.lammps file by adding the verbose keyword. These informations include the following.

The "path angle" (pathangle) for the replica i which is the angle 

between the 3N-length vectors (Ri-1 - Ri) and  (Ri+1 - Ri) (where Ri is the

atomic coordinates of replica i). A "path angle" of 180 indicates that replicas

i-1, i and i+1 are aligned.

 "angletangrad" is the angle between the 3N-length tangent vector and

the 3N-length force vector at image i. The tangent vector is calculated as in "(HenkelmanA)"_#HenkelmanA for all intermediate replicas and at R2 - R1 and RM - RM-1 for the first and last replica, respectively.

"anglegrad" is the angle between the 3N-length energy gradient vector of replica i and that of

replica i+1. It is not defined for the final replica and reads nan.

gradV is the norm of the energy gradient of image i.

ReplicaForce is the two-norm of the 3N-length force vector (including nudging forces) for replica i.

MaxAtomForce is the maximum force component of any atom in replica i.

When a NEB calculation does not converge properly, these suplementary

informations can help understanding what is going wrong. For instance when the

path angle becomes accute the definition of tangent used in the NEB

calculation is questionable and the NEB cannot may diverge

"(Maras)"_#Maras.

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

There are 2 Python scripts provided in the tools/python directory,

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 labeled with a specific atom type.

the NEB atoms labelled with a specific atom type.

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

package.  See the "Making LAMMPS"_Section_start.html#start_3 section

for more info on packages.

:line

[Related commands:]

"prd"_prd.html, "temper"_temper.html, "fix

[(HenkelmanB)] Henkelman, Uberuaga, Jonsson, J Chem Phys, 113,

9901-9904 (2000).

:link(Nakano3)

:link(Nakano)

[(Nakano)] Nakano, Comp Phys Comm, 178, 280-289 (2008).

:link(Maras)

[(Maras)] Maras, Trushin, Stukowski, Ala-Nissila, Jonsson, Comp Phys Comm, 205, 13-21 (2016)

class FixNEB : public Fix {

 public:

  double veng,plen,nlen;

  double veng,plen,nlen,dotpath, dottangrad,gradlen,dotgrad ;

  int rclimber;

  double gradvnorm;

  FixNEB(class LAMMPS *, int, char **);

  ~FixNEB();

 private:

  int me,nprocs,nprocs_universe;

  double kspring;

  double kspring, kspringPerp,EIniIni,EFinalIni;

  bool StandardNEB, NEBLongRange,PerpSpring, FreeEndIni,FreeEndFinal,FreeEndFinalWithRespToEIni,FinalAndInterWithRespToEIni ;

  int ireplica,nreplica;

  int procnext,procprev;

  int cmode;

  MPI_Comm uworld;

  MPI_Comm rootworld;

  char *id_pe;

  class Compute *pe;

  int nebatoms;                // # of active NEB atoms

  int nebatoms;

  int ntotal;                  // total # of atoms, NEB or not

  int maxlocal;                // size of xprev,xnext,tangent arrays

  double **xprev,**xnext;      // coords of my owned atoms in neighbor replicas

  double **tangent;            // work vector for inter-replica forces

  double *nlenall;

  double **xprev,**xnext,**fnext,**springF;

  double **tangent;

  double **xsend,**xrecv;      // coords to send/recv to/from other replica

  double **fsend,**frecv;      // coords to send/recv to/from other replica

  tagint *tagsend,*tagrecv;    // ditto for atom IDs

                                 // info gathered from all procs in my replica

  double **xsendall,**xrecvall;    // coords to send/recv to/from other replica

  double **fsendall,**frecvall;    // force to send/recv to/from other replica

  tagint *tagsendall,*tagrecvall;  // ditto for atom IDs

  int *counts,*displacements;   // used for MPI_Gather

  void inter_replica_comm();

  void reallocate();

};

}