NEB Upgrade (Nov. 2013)

Hi Steve/ LAMMPs community,

I have been trying to use the neb technique with the newly added “each” syntax, with the very simple case of the formation of an adatom and a subsequent “hop” of that atom to an adjacent position. The simulation runs, but produces a rubbish output. I do not understand what would cause this - perhaps you can elaborate Steve?

Firstly, I modified the in.neb.hop1 script so that the final position is final.hop2. Then, I modified the dump commands by simply modifying the dump output lines:

$ dump 1 nebatoms custom 200 dump.neb.* id x y z
$dump 2 nonneb custom 200 dump.noneb.* id x y z

To produce dump.X.0 - dump.X.1141 dumpfiles. I then modified these files by deleting the headers (a requirement that I find rather annoying! - why can it not automatically identify the id, and correlate the adjacent numerical values with x, y and z?


Then, I create a new file, based on a very minor modification of in.neb.hop1 to create the file in.neb.Largehop1:

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 i index 200 400 600 800 1000 1200 1400 1600 1800 2000

create geometry with flat surface

log log.My_neb_1

lattice hex 0.9
region box block 0 20 0 10 -0.25 0.25
create_box 3 box
create_atoms 1 box

mass * 1.0

LJ potentials

pair_style lj/cut 2.5
pair_coeff * * 1.0 1.0 2.5
pair_modify shift yes

initial minimization to relax surface

minimize 1.0e-15 1.0e-15 50000 500000
reset_timestep 0

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
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.01 50000 50000 100 each dump.noneb.$i

This converged with a reverse energy barrier of 0!?! Furthermore, the forwards energy barrier was ~30% less % of the energy barrier obtained for hop1 alone!.

I re-attempted this exact same methodology, just for hop1… I got even more rubbish results - an energy barrier that was approaching E-05 eV!?!

Does anyone know if the each command can even actually be used to produce a viable result?

In response to queries a) and b) - As I am still not obtaing a consistent output even for the basic examples provided, I do not believe that I am in any manner a reliable source. However, when I run the simulation mentioned above with a timestep:

timestep 0.0001

I obtained exactly the same result as with timestep 0.05.

in response to part c, it would be extremely valuable to enable investigations of strained samples with PBCs (avoiding free surface effects)…? But that is essentially just theoretical.

Thanks - I would really appreciate advice on running neb with the each method… If it is possible?


As per my previous email, I am wondering about the accuracy of the NEB command. Referring again to the Example of 2 “adatom hops”, I have created a sequence of dump files which list the id x y z information of the NEB atoms, every 200 timesteps. I then appropriately modified these files in accordance with the format required (as specified: These dumpfiles were also renamed so they could be accessed via partition variables (i.e. they were dumped at every 200 timesteps as " .* ").

When I use the NEB command to attempt to find the MEP from the initial state to the final (2-hops) state, one would anticipate that the energy barrier should be somewhat LESS than what would be required to overcome two sequential steps (noting that the MEP involves 11 intermediate “steps”). However, I find that the energy barrier provides an substantial over-estimation (i.e., 0.01924708 compared with 0.0071048579+0.0016018802)…

Please inform me if it is simply impossible to use the NEB command except for very small/ single-step atomic transitions… Please, also check out my methods. Here is the input file and I have attached the necessary dumpfiles:

You will need to run this simulation with 11 partitions… mpirun -np 11 ./lmp_NeSI_26.02.2014 -partition 11x1 -in in.neb.hop1


dump.neb.11 (484 Bytes)

dump.neb.10 (482 Bytes)

dump.neb.9 (477 Bytes)

dump.neb.8 (475 Bytes)

dump.neb.7 (483 Bytes)

dump.neb.6 (479 Bytes)

dump.neb.5 (484 Bytes)

dump.neb.4 (484 Bytes)

dump.neb.3 (484 Bytes)

dump.neb.2 (483 Bytes)

dump.neb.2 (483 Bytes)

dump.neb.0 (616 Bytes)

I’ll check it out. I saw your original email,
but haven’t got to this issue yet. Thanks
for posting the files.


Hi Steve,

Would you please re-consider this issue, as an NEB analysis that tracks the evolution as a sequence of dumpfiles based on an original LAMMPs simulation could enable particularly interesting energetic analysis…

Kind regards,

yes, will hopefully get to it this week

sorry for the delay,