Regarding the use of the NEB technique

LAMMPS LAMMPS Mailing List Mirror
1

Hi LAMMPs community,

I am somewhat perplexed by a recent observation that I have made using the NEB technique.

For a simple validation of the technique, I have simulated the diffusion of a single atom into a vacancy within bulk Aluminium. Rather oddly, the original estimation for the MEP (i.e., at timestep ~0), is very close to the expected value, being ~1.23eV. This is a good approximation of the theoretical vacancy migration energy of pure aluminium from first principles (i.e., 1.27eV [1]) and from experiments at 360-480K (1.31 eV [2]). However, once the minimization and climbing NEB phases are complete, the energy is dramatically smaller than expected, at ~0.63eV. ALSO, the number of partitions significantly influences the magnitude of the FINAL, minimized energy barrier, but not the INITIAL.

Please refer to the attached txt file, for the simple case of 1 atom diffusing in pure Aluminium (using Al99.eam.alloy potential file). This compares an NEB with 40 replicas and an NEB with 16 replicas.

Perhaps more importantly, in my case, I have found that the NEB technique will form a good MEP at timestep 0, with a clear reverse energy barrier. However, after fully minimized, there is no reverse energy barrier. This result is consistent for several different cases that I have tested. This does not seem likely, because the end state is fully minimized using both the cg and quickmin minimizers to an energy tolerance of 1e-25 (energy IS the stopping criterion)! It may be possible that, because the end state is a higher energy than the initial, it is only a metastable state. However, in light of the results for my “benchmark” simulation using NEB to calculate the MEP for single atom diffusion, maybe it is more accurate to assume the energy path at t=0 to be a better MEP, which does have a clear reverse energy barrier (EBR).

Is it just because the interpolation inherently forms a bell-shape curve? I am not so sure, because the magnitude of the energy barrier is HIGHLY dependent on the initial and final states (and NOT just related to the potential energy difference). You can refer to the attaced Word file for a clear comparison of the different EPs.

ALSO NOTE: The final energy barrier at t=0 is mostly independent of the spring constant used, but is somewhat influenced by the timestep and very much by the number of replicas.

Can anyone with experience using the NEB technique comment on this observation? Is it best to just assume the initial energy path (at t=0) is the best!?!

Thanks!!!

16_Partition_1atomDiffuse_Al99.eam.alloy.txt (2.37 KB)

40_Partition_1atomDiffuse_Al99.eam.alloy.txt (4.81 KB)

NEB_Weird_Results_For_LAMMPs_Community.docx (266 KB)

2

The results in the text files look fine to me and do not exhibit the
behavior you complain about in your message or in the Word file. For
example, the second file reports the forward and reverse energy
barriers to be:

0.6349601 0.63496006

Perhaps you just need to spend a little more time familiarizing
yourself with the meaning of the LAMMPS output. I suggest you read
this doc page:

http://lammps.sandia.gov/doc/neb.html

Aidan

3

Hi Aidan,

I was not complaining at all, simply seeking an explanation for apparently inconsistent results. I understand entirely which values refer to the EBF and EBR.

It is unusual that the final converged energy barrier is dependent on all of the parameters - the number of replicas, the spring constant, the timestep, the energy tolerance, the number of NEB atoms (*i.e. atoms which experience a spring force)… The initial interpolation IS consistent, and actually matches the experimental value. This is the reason for my comment, not because I do not understand the output!! Finally, the referenced EBF should be closer to 1.27 eV, which is closely approximated by the initial state.

Since you refer to the text files, perhaps you would like to refer to the text files for the Word examples?

Another compelling example, is my Case 3. They consistent of identical initial and final states. However, depending on the number of atoms that is used for the fix neb command, the outcome is enormously different. Based on only the NEB atoms (the defect atoms and their neighbours), I get a positive “plateau” converged energy barrier. If I use all atoms, it converges with a NEGATIVE energy barrier, and plateaus below the initial energy state.

Perhaps you can “interpret” these results for me?

Nathaniel

CASE3_WITH_NEBATOMSONLY.txt (4.81 KB)

CASE2_LAMMPS_COMMUNITY.txt (4.8 KB)

CASE3_WITH_ALL_ATOMS.txt (4.81 KB)

CASE1_LAMMPS_COMMUNITY.txt (4.81 KB)

4

Hi again,

I should state that I re-checked the literature and found that the activation energy for vacancy migration is roughly 0.60eV, so my result is reasonable. I had made an error previously, because I was referring to an old literature review I had performed (the 1.27eV also includes the energy of vacancy formation).

However I still do not understand why the final result is so substantially influenced by the number of replicas. I did not include the case where I used 4 replicas…

In other words:

4 partitions: 0.461eV
16 partitions: 0.5959 eV
40 partitions: 0.6350 eV
1000 partitions: ???

Infinite partitions: asymptote, or infinity??? ( I would hope asymptote is the case… but is infinite partitions practical - obviously not)?

Other parameters are also marginally important. However, even the energy tolerance only marginally influences the result (whereas the number of replicas, or number of “nebatoms” dramatically influences the results.

Can anyone refer to a good way to validate the results obtained? Or is this just a tool to compare different cases, maintaining all NEB factors as constant in each case?

thanks,
Kind regards,
Nathaniel

5

Hi Nathaniel,

So, the txt and doc files are from different runs? Glad to hear it.
You should try to only ask one question at a time. The shotgun
approach only works when your resources are cheap, definitely not the
case here. I looked at all your txt files and observed a few things:

1. The two 1atom runs both worked well, so that's good. They exhibit
several nice properties that you should expect to see in any
well-converged NEB run:
        a. The climbing replica is near the middle of the band
        b. MaxReplicaForce is small
        c. MaxReplicaForce is the same as GradVc i.e. all replicas are
better converged than the climbing replica
        d. GradV0 and GradV1 are very small

2. Your complaint about needing 1000 replicas is not supported. The 16
and 40 replica runs are not highly converged. If you run them longer,
the differences will get smaller.

3. Apparently your result with 1 atom is surprisingly close to the
literature value (DFT?).

4. Your Case 3 runs are both flawed. The initial gradient in the last
replica GradV1 is large. I can't say for sure, but it is possible
that the last replica found a way to relax towards the first replica,
thereby eliminating the barrier. You should try to relax the end state
before starting NEB.

5. Another problem with Case 3 All Atoms is that the last replica is
lower in energy than the first replica. This suggests that your
presumed stable initial state is really not very stable at all, when
all the atoms are allowed to move.

6. One more comment. The space of possible minimum energy pathways
between two minima in a system with 3N degrees of freedom is huge.
NEB is trying to find one of those MEPs. There is no guarantee it will
find "the right one." In other words, from a purely mathematical
point of view you are trying to do something that is impossible. A
lot of humility is required, but it is still possible to make
progress. Start with the 1atom result, and figure out how you can add
more degrees of freedom and still converge to the right answer. Hint,
consider using something other than file-style final.

There is probably lots more to be learned about this system, but the
mailing list is not a good place to do that. You need to think harder
about what you are seeing, or else just stick with the nice simple 1
atom result.

Aidan