UT theoretical chemistry code forum

Hello all

I have been doing some simple NEB calculations for diffusion of H from an FCC to a HCP binding site on a Os(0001) closed packed surface. The INCAR I am using is as follows:
# Regular
SYSTEM = Os_NEB
ISYM = 0
NSW = 2000
EDIFF = 0.1E-6
EDIFFG = -0.01
ENCUT = 354
GGA = RP
LWAVE = .TRUE.
LCHARG = .FALSE.
PREC = medium
ISIF = 0
# ISPIN = 2

# Speedup
NPAR = 2
NSIM = 4
IALGO = 48
LPLANE = .TRUE.
LREAL = Auto

#NEB
ICHAIN = 0
IMAGES = 8
SPRING = -5.0
LCLIMB = .FALSE.
LTANGENTOLD = .FALSE.
LDNEB = .FALSE.

# Optimizer2
IBRION = 3
POTIM = 0.00
IOPT = 2
MAXMOVE = 0.10
FDSTEP = 5E-3

And the neb.dat file after using nebresults.pl script gives:
0 0.000000 0.000000 -0.001496 0
1 0.181737 0.017571 -0.002813 1
2 0.368447 0.037901 0.001344 2
3 0.551387 0.046485 0.000317 3
4 0.729682 0.049821 0.000008 4
5 0.906971 0.052864 -0.000198 5
6 1.086388 0.055378 -0.001679 6
7 1.272234 0.051754 -0.002708 7
8 1.463062 0.033695 0.008068 8
9 1.643894 0.018108 0.000795 9

What is bugging me is the look of the minimum energy path figure. As you can see from the neb.dat the force projection on to tangent for each image is very small but I thought that unless it is a minima or maxima it should be a little bit larger and I thought that for a smooth curve an increase in energy would have positive force and decrease would have a negative force but here it is just random with negative forces both before and after the transition state etc.

At least would anyone please look at my mep.eps figure and hopefully give me some advice of what might be wrong.

I have done similar calculations for other metals in Dacapo and there I always got a smooth curve, not with these tiny bumps all over the ME path.

Thanks,
sig

p.s. I attached the mep.eps image as well as the movie.xyz generated using nebresults.pl

For some reason, the forces do not appear to be correct -- the 3rd column of numbers are essentially all zero. These values determine the slope of the spline interpolation used to make the mep plot. This is why the spline is flat when it passes through each image point.

I have not seen this before and I'm not sure what is going wrong. Perhaps check that the vtstcode is linked into your vasp binary. I would also be happy to look at a tar.gz file of the run and see what went wrong. I imagine that there is a simple solution.

Thanks Graeme for a quick reply.

Attached is tar.gz of the results folder for the NEB calculation. I wasn't sure what you needed exactly from this all.

Any advice or suggestion of what might be wrong would be of great help. I have done a lot of NEB calculations in the past and never seen this either. I was wondering if we have not updated the scripts from you that generate the results (nebresults.pl) after we started using the latest version of vasp and perhaps read the wrong forces from the output or something.

Thank you very much for your help,

Sigga

Ah, that was confusing, but I see the problem.

If you look at your constraints, the first and last atoms are frozen in x and y, and free to move in z. The H atom is particularly problematic because it is moving in x and z along the band. In our implementation of the NEB, any forces which are in constrained degrees of freedom are zeroed. Our assumption is that frozen atoms are held frozen in the same position across the band. When this is violated, as it is in your band, there will be components of the NEB tangent in frozen degrees of freedom, but the force is set to zero so that the force projected onto the tangent will not be consistent with the change in energy.

In principle, we could treat this situation in a different way in our code, such as projecting the force onto the tangent before zeroing the force on frozen degrees of freedom. But I assume that this is just a mistake, and your intension was to have the first and last atoms free to move in all degrees of freedom.

A separate issue: check the distance between your two frozen layers. I see there are high forces repelling those layers. They are frozen so it is not a big problem in your calculation, but it probably indicates a poor lattice constant.

Ah thank you very much. I will look into this more closely.

Thank you very very much for your help.

Sigga