UT theoretical chemistry code forum

I am trying to use cNEB to find the MEP for a molecule bridging two adatoms to go from the flat, bridging geometry to an upright geometry. I have converged both the initial and final state geometries to the convergence level that I am using for the cNEB. The difference in energy between the two states is ~0.75 eV. It appears that the saddle point is very close to the initial state. This makes sense because a bond would have to be broken before the molecule can rotate to the upright geometry. I am not familiar enough with NEB calculations to know if this could be the cause of my issues (seems like this could be a common phenomenon when searching for a MEP). The calculation seems to have been going ok at first and then seems to become "unstable" after ~250 steps (see attached images). The INCAR is given below:

INCAR
------------------------------------------------------------------------------
SYSTEM = flat to upright NEB

PREC = Normal
ENCUT = 400
ALGO = Fast
NELMDL = -3
LREAL = Auto
LWAVE = .FALSE. # do not write WAVECAR file
ISMEAR = -5

IMAGES = 5 # number of images, forces NEB method
LCLIMB = .TRUE. # use climbing image method
IOPT = 1 # LBFGS = Limited-memory Broyden-Fletcher-Goldfarb-Shanno
IBRION = 3 # tell VASP to do MD
POTIM = 0 # tell VASP zero time step

NELMIN = 4 # do a minimum of four electronic steps
EDIFF = 1E-5 # high accuracy
EDIFFG = -0.03 # accuracy of ions low
NSW = 1000 # 100 ionic steps in ions

LPLANE = .TRUE.
LSCALU = .TRUE.
NSIM = 4
NPAR = 1
------------------------------------------------------------------------------

Is there anyone who has an idea as to what is causing this? Thanks in advance for any advice.

Mike Garvey
University of Wisconsin-Milwaukee

We would have to see the output to debug this further. For a start, however, try running with a conservative optimizer, such as IOPT=3. If this has the same problem, it would be good to look into this further.

But also, check the geometry when you see the forces start to rise. You may find that the band has found a lower energy pathway, and is falling off a shoulder towards it.

Professor Henkelman,

Thanks for your advice. I was going to let it keep running with the LBFGS optimizer since the forces and energies had started coming back down. However, I accidentally ran the nebresults.pl script in the directory where the job was running and crashed it (rookie mistake!). I have attached the MEP plot from the unconverged calculation for reference. As I said previously, it appears that the saddle point is very close to the initial state. I have also included the vaspout plot for image 1. If you have a chance, could you take a quick look at the plots and see if anything looks wrong?

I have decided to start the job over with the Quick-Min optimizer as you suggested. I have also changed EDIFF from 1E-5 to 1E-6 to make sure the forces are accurate enough. If that also has problems, I will definitely contact you with the output for help in determining the cause.

Thanks again for your help.

Professor Henkelman,

I started the calculation from the beginning using IOPT=3 as I said in my previous post. It has run almost 100 ionic steps and seems to be having trouble. The max force for all images started increasing from the beginning and stabilized at ~10 eV/Ang. The energy started climbing after ~12 steps for all images and started to level off at ~30 steps (it has increased by ~0.7-0.9 eV). Both the energy and max force for each image have not changed much for the last 70 steps. I have attached the vaspout.eps files so you can see the behavior more clearly. You mentioned that you might be willing to look at the output to figure this out. If so, how would you like me to get the output to you? I could post it here or email it directly to you. Please let me know what you would prefer.

Thanks,

Mike Garvey

Something doesn't look right. A .tar.gz file of the entire calculation (minus the CHG* WAV* files) would be the best for debugging this.

Professor Henkelman,

I have made tar.gz files of the two calculations (LBFGS and QM), but I am unable to upload them. They are ~100MB and ~33MB. Is that too large to post as an attachment?

Mike Garvey

My bad; I had a large upload file size limit, but my post size limit was too small. It should work now for files up to 256 MB.

Professor Henkelman,

That seems to have worked. Thank you.

Mike Garvey

These are cool calculations! I don't see any obvious problems with your set up; I'm going to try a few changes, and then I'll report back.

I think the problem you were seeing may just be due to a stiff C-N bond. Reducing the time step for quickmin, and/or using a lower INVCURV for LBFGS seem to give stable convergence of the NEB calculations. You can see in the attached calculations. I used cheaper settings (which I always recommend when doing NEB calculations, since it is easy to restart then with more accurate settings), using gamma point and a lower cutoff with softer potentials. The minr and minp are the reactant and product. You can see I did tests with fire and lbfgs and they converge as well. Then in the neb directory, I did one conservative run (in run1) using QM with a small time step of 0.001. I really didn't need to run for 1000 steps here, this was just the default that you had. Then I switched to lbfgs and the forces are dropping nicely.

I was a little worried that you might be missing the saddle at the start of the path since it is so long, so I did another neb calculation (in nebint) going between the initial state and the first image. Remarkably, there does not appear to be a barrier there, consistent with the full band. There is a pretty neat bond breaking, however, with some nice bending of the molecule coupled with motion of the Au adatom.

Professor Henkelman,

Thanks for all of the advice and help. Could you tell me what software you use to visualize the movie outputs? I currently use XCrySDen to make images of each step and then combine them in GIMP to make an animated gif. It works, but I doubt it is the most efficient method.

Mike Garvey

We have a local viewer, written by Rye Terrell, which can be downloaded as part of our tsase package;
see http://theory.cm.utexas.edu/henkelman/code/

One nice thing about this viewer is that you can run it across a network over x-windows.

Another option is VMD.

Hey Graeme,

How did you compile the tsase software onto your mac? Currently, I am unable to run the GUI through X-windows because of the absence of gtk....? I have attempted to brew install pygtk but with no luck. Any suggestions would be helpful! :-)

Hi Paul, I've found the macports is the best way to get these dependences working on the mac.
Graeme

Graeme,

Thank you for your last. However, upon searching gtk via macports 220 potential ports appear in my terminal. Moreover, there are a number of gtk options to choose from which are not descriptive enough to suggest selecting one to download over another. So what gtk have you ported onto your mac? I am using Mac OS X 10.8.