Merge pull request #490 from EmileMaras/NEB-Change

[Syntax:]

fix ID group-ID neb Kspring :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

fix ID group-ID neb Kspring keyword value :pre

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

neb = style name of this fix command :l

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

zero or more keyword/value pairs may be appended :l

keyword = {nudg_style} or {perp} or {freend} or {freend_k_spring} :l

  {nudg_style} value = {neigh} or {idealpos}

    {neigh} = the parallel nudging force is calculated from the distances to neighbouring replicas (in this case, Kspring is in force/distance units)

    {idealpos} = the parallel nudging force is proportional to the distance between the replica and its interpolated ideal position (in this case Kspring is in force units)

  {perp} value {none} or kspring2

    {none} = no perpendicular spring force is applied

    {kspring2} = spring constant for the perpendicular nudging force (in force/distance units)

  {freeend} value = {none} or {ini} or {final} or {finaleini} or {final2eini}

    {none} = no nudging force is applied to the first and last replicas

    {ini} = set the first replica to be a free end 

    {final} = set the last replica to be a free end

    {finaleini} = set the last replica to be a free end and set its target energy as that of the first replica

    {final2eini} = same as {finaleini} plus prevent intermediate replicas to have a lower energy than the first replica

  {freeend_kspring}  value = kspring3

    kspring3 = spring constant of the perpendicular spring force (per distance units)    

   :pre

[Examples:]

fix 1 active neb 10.0 :pre

fix 1 active neb 10.0

fix 2 all neb 1.0 perp 1.0 freeend final

fix 1 all neb 1.0 nudg_style idealpos freeend final2eini freend_kspring 1:pre

[Description:]

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

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 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:

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) but the climbing replica the force vector

becomes:

Fi = -Grad(V) + (Grad(V) dot That) That + Fnudgparallel + 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.  Fnudgparallel is an

artificial nudging 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).

The keyword {nudg_style} allow to specify how to parallel

nudging force is computed. With a value of idealpos, the spring 

force is computed as suggested in "(E)"_#E :

Fnudgparallel=-{Kspring}* (RD-RDideal)/(2 meanDist) :pre

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.

When {nudg_style} has a value of neigh (or by default), the parallel 

nudging force is computed as in "(Henkelman1)"_#Henkelman1 by 

connecting each intermediate replica with the previous and the next 

image:

Fnudgparallel= {Kspring}* (|Ri+1 - Ri| - |Ri - Ri-1|) :pre

The parallel nudging force associated with the key word idealpos should

usually be more efficient at keeping the images equally spaced.

:line

The keyword {perp} 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)"_#Maras1).

The perpendicular spring force is given by

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

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

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)*kspring3) That,  {when} Grad(V) dot That < 0

Fi = -Grad(V)+ (Grad(V) dot That + (ETarget- E)*kspring3) That, {when} Grad(V) dot That > 0

:pre

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 {finaleini} or {final2eini} 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 {finaleini}, 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 {final2eini}

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 {final2eini} 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 {finaleini} or

{final2eini} cannot reach the convergence criteria.

:line

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.

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

"neb"_neb.html

[Default:] none

[Default:]

The option defaults are nudg_style = neigh, perp = none, freeend = none and freend_kspring = 1.

:line

:link(Henkelman1)

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

:link(Henkelman2)

[(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(Maras1)

[(Maras)] Maras, Trushin, Stukowski, Ala-Nissila, Jonsson,

Comp Phys Comm, 205, 13-21 (2016)

mpirun -np 4 lmp_g++ -partition 4x1 -in in.neb.hop1

mpirun -np 4 lmp_g++ -partition 4x1 -in in.neb.hop2

mpirun -np 4 lmp_g++ -partition 4x1 -in in.neb.hop1freeend

mpirun -np 3 lmp_g++ -partition 3x1 -in in.neb.sivac

mpirun -np 8 lmp_g++ -partition 4x2 -in in.neb.hop1

mpirun -np 8 lmp_g++ -partition 4x2 -in in.neb.hop2

mpirun -np 8 lmp_g++ -partition 4x2 -in in.neb.hop1freeend

mpirun -np 6 lmp_g++ -partition 3x2 -in in.neb.sivac

mpirun -np 9 lmp_g++ -partition 3x3 -in in.neb.sivac

Note that more than 4 replicas should be used for a precise estimate 

of the activation energy corresponding to a transition.

If you uncomment the dump command lines in the input scripts, you can

create dump files to do visualization from via Python tools: (see

lammps/tools/README and lammps/tools/python/README for more info on

group		nonneb subtract all nebatoms

fix		1 lower setforce 0.0 0.0 0.0

fix		2 nebatoms neb 1.0

fix		2 nebatoms neb 1.0 nudg_style idealpos

fix		3 all enforce2d

thermo		100

# 2d NEB surface simulation, hop from surface to become adatom

dimension	2

boundary	p s p

atom_style	atomic

neighbor	0.3 bin

neigh_modify	delay 5

atom_modify	map array sort 0 0.0

variable	u uloop 20

# create geometry with flat surface

lattice		hex 0.9

region		box block 0 20 0 10 -0.25 0.25

read_data        initial.hop1freeend

# LJ potentials

pair_style	lj/cut 2.5

pair_coeff	* * 1.0 1.0 2.5

pair_modify	shift yes

# define groups

region	        1 block INF INF INF 1.25 INF INF

group		lower region 1

group		mobile subtract all lower

set		group lower type 2

timestep	0.05

# group of NEB atoms - either block or single atom ID 412

region		surround block 10 18 17 20 0 0 units box

group		nebatoms region surround

#group		nebatoms id 412

set		group nebatoms type 3

group		nonneb subtract all nebatoms

fix		1 lower setforce 0.0 0.0 0.0

fix		2 nebatoms neb 1.0 nudg_style idealpos freeend ini

fix		3 all enforce2d

thermo		100

#dump		1 nebatoms atom 10 dump.neb.$u

#dump		2 nonneb atom 10 dump.nonneb.$u

# run NEB for 2000 steps or to force tolerance

min_style	quickmin

neb		0.0 0.1 1000 1000 100 final final.hop1