UT theoretical chemistry code forum

I have been running NEB for C4H3(CH2OH)O* + * ---> C4H3(CH2)O*+ OH* step, the highest image has got forces lesser than -0.05 ev, my EDIFFG criterion, which i checked from outcar which reads as

stress matrix after NEB project (eV)
-52.47845 -0.42376 -0.07734
-0.42376 -51.58503 0.04778
-0.07734 0.04778 -5.94655
FORCES: max atom, RMS 0.006373 0.001843
FORCE total and by dimension 0.016169 0.004277
Stress total and by dimension 73.829130 52.478450
however, the NEB doesnot stop. Please suggest me what I am doing wrong in here.
I have attached the calculation here.

Thanks

This looks ok. The force on every image has not dropped below 0.05 eV/Ang, but if you are primarily interested in the saddle, you can see that it is well-converged in image 01 with a barrier of 0.3091 eV.

Note that your list of frozen atoms is messed up on image 00.

Yeah so I optimised the highest structure separately with very low POTIM. However, the optimised structure gave multiple imaginary frequencies like -100,-120 and -290 (-100 and -120 corresponds to vibration of the furan ring of C4H3(CH2)O*) which conveys that it is not a first order saddle point.

Is there any way that I can remove these negative frequencies so that I have a first order transition state here ?
I tried perturbing the TS state a liitle bit in structure in the direction of smaller frequencies(-100) but didnot work.

Thanks

Frequency calculations are a little tricky. Make sure that you use an ediff value around 1e-8. And the linear response routine in vasp is best, since it avoids issues related to the finite difference step size.

If you need to converge the TS structure even tighter, the neb2dim.pl script will set up a dimer calculation which will let you converge to the saddle with just a single image.

Thanks for the suggestion.
However, I didnot understand what u mean by using linear response routine in vasp. And do I need to use ediff 10-8 for frequency calculation ?
I am using this INCAR :

PREC = NORMAL
LREAL = AUTO
VOSKOWN = 1
EDIFF = 1E-04

ALGO = Very_Fast ##Updated 8/07 Very_Fast = IALGO 48, Normal = IALGO 38, Fast = some of each
NELM = 400
ISYM = 0
ISMEAR = 2; SIGMA = 0.2 ## broadening in eV. Should be about 0.2 for metals, 0.003 for molecules and semiconductors (oxides)

ISPIN = 2 #spin polarized calculation? 1= no, restricted, 2=yes, unrestricted
ENCUT = 400
EDIFFG = -.5E-01 stopping-criterion for IOM

NSW = 1000 number of steps for IOM
IBRION = 5 ionic relax: 0-MD 1-quasi-New 2-CG 5-frequency calc
ISIF = 2
POTIM = 0.01
LWAVE = .FALSE.
LCHARG = .FALSE.
NWRITE = 3
NFREE =2

Thanks

that's not going to work. you need accurate electronic structure for curvatures - set ediff to 1e-8
then use IBRION = 7 or 8

Thanks a lot !!
I will try that and see.

As for dimer calculation on the above mentioned state, I was getting following high forces and curvatures (so I aborted it):
(I used the default parameter values for Ichain = 2 with IOPT = 2)
Step Force Torque Energy Curvature Angle
1 0.03822 27.81351 -425.38946 21.07626 27.07616
1 0.03822 13.40899 -425.38946 6.64038 18.84903
1 0.03822 7.08688 -425.38946 2.80527 9.29710
1 0.03822 5.14913 -425.38946 1.81077 4.85174
2 0.84805 4.19966 -425.35825 1.57256 4.79895
2 0.84805 3.87299 -425.35825 1.23578 3.75358
2 0.84805 3.01986 -425.35825 1.02241 3.15127
2 0.84805 2.76094 -425.35825 0.82356 2.98896
3 1.34556 2.89367 -425.30298 0.86278 3.18985
3 1.34556 2.61263 -425.30298 0.69823 2.55405
3 1.34556 2.03671 -425.30298 0.53737 2.16842
3 1.34556 1.90963 -425.30298 0.49610 2.23910
4 1.68822 1.95925 -425.24266 0.39856 1.64464
4 1.68822 1.42717 -425.24266 0.31315 2.08230
4 1.68822 1.65212 -425.24266 0.28513 1.26501
4 1.68822 0.99739 -425.24266 0.22338 1.58642
5 1.91425 1.61310 -425.18156 0.49423 2.03362
5 1.91425 1.67383 -425.18156 0.43287 1.46765
5 1.91425 1.27533 -425.18156 0.41475 1.83991
5 1.91425 1.58606 -425.18156 0.37772 1.46969
6 2.09940 1.56038 -425.12133 0.27439 2.49506
6 2.09940 1.89886 -425.12133 0.20083 1.43574
6 2.09940 1.15224 -425.12133 0.15668 2.42761
6 2.09940 1.88362 -425.12133 0.11606 1.39666
7 2.23919 1.45373 -425.06607 0.03768 6.24354
7 2.23919 2.94214 -425.06607 -0.02154 2.32090
7 2.23919 1.65561 -425.06607 -0.16303 3.23083
7 2.23919 2.69940 -425.06607 -0.25987 2.32865
8 1.57655 2.57471 -425.04787 -0.24294 5.91015
8 1.57655 4.39905 -425.04787 -0.42970 3.16877
8 1.57655 2.54278 -425.04787 -0.69916 6.78064
8 1.57655 5.40921 -425.04787 -1.02489 3.97031

Let me hear your suggestions on this dimer calculation.
Thanks

Ok, increase DRotMax to about 20 and make sure that the curvature goes negative in the first step. In this calculation it did not find the negative mode quickly enough and so it wandered away from the saddle and then would have to come back. In your calculation the lowest modes drops from 21 to 1.8 in the 4 rotation steps, but it is critical that it goes negative to rapidly converge to the saddle.

As a side note, I recently fixed a bug in how the negative mode was being stored and read upon restart. It will not be an issue if you do the entire calculation in a single run, but if you are restarting a dimer calculation, it would make sense to do the update.

Thanks for the suggestion.
However, I ran vibration calculation with IBRION 8, but it ended in 1 NSW saying it cannnot be used with vanderwaals correction ON.
Please suggest.
Thanks

Ah, then yes, you may have to resort to the finite difference calculation. But the most important thing is that you need accurate forces by setting ediff=1e-8.

Thanks

As for the bug as u mentioned below:

"As a side note, I recently fixed a bug in how the negative mode was being stored and read upon restart. It will not be an issue if you do the entire calculation in a single run, but if you are restarting a dimer calculation, it would make sense to do the update."

do you mean I need to update my vasp_vtst with the instructions from this page http://theory.cm.utexas.edu/vtsttools/installation.html for dimer calculation ?

yes, that's right