UT theoretical chemistry code forum

Hello,

I am attemtping to run some NEB/Dimer calculations with a large slab (181 atom total).
In order to make these calculations practical, I have been forced to add;
ALGO=FAST
LREAL=A
NSIM=4
to get any results at a reasonable time.
However I wonder what these settings will be doing to the code produced by Professor Henkelman and his group, in terms of convergence or reaching the correct transition state.

As of this moment one NEB calculation and two dimer calculations for the same path have also produced a qualitatively similiar transition state stucture but none have converged.
I had hoped (naievly?) that I could produce the correct TS with these settings and then corrected the energy with a subsequent calculation using default VASP settings.

Is it possible that the use of the 'fast settings' is stopping the calculation from proceeding correctly? Certainly it is making the calculation proceed more quickly. Or do I just need to wait for a longer time period for the TS to be found?

Many Thanks for any reply,



---Results of DIMCAR file for one calculation (resumed by replacing POSCAR with CONTCAR, and MODECAR by NEWMODECAR)---
Run: No1
Step Force Torque Energy Curvature Angle
1 0.93029 4.23658 -841.66558 -1.14820 6.88745
1 0.93029 4.04594 -841.66558 -1.33522 7.54077
1 0.93029 3.93909 -841.66558 -1.55116 3.99649
1 0.93029 3.33027 -841.66558 -1.58875 5.64526
2 0.76996 5.82522 -841.66937 -1.62759 7.69137
2 0.76996 6.19712 -841.66937 -1.78638 3.20428
2 0.76996 5.92385 -841.66937 -1.81422 1.61379
2 0.76996 5.83550 -841.66937 -1.82668 2.21946
3 0.87115 6.84286 -841.68563 -1.64484 -1.89902
3 0.87115 6.62387 -841.68563 -1.66639 2.38831
3 0.87115 6.82852 -841.68563 -1.66687 -1.04269
3 0.87115 6.78812 -841.68563 -1.66393 -0.12660
4 0.56593 7.51658 -841.69168 -0.17688 -5.26238
4 0.56593 7.51609 -841.69168 -0.15830 -1.58383
4 0.56593 7.57281 -841.69168 -0.13400 -0.12663
4 0.56593 7.57886 -841.69168 -0.13131 -0.04012
5 0.94247 7.64980 -841.69745 -0.71144 0.93606
5 0.94247 7.53578 -841.69745 -0.75951 -0.12281
5 0.94247 7.61227 -841.69745 -0.74188 1.08147
5 0.94247 7.62293 -841.69745 -0.80783 -0.37453
6 0.59822 7.73434 -841.70186 0.10583 -2.13714
6 0.59822 7.66185 -841.70186 0.11951 -0.44254
6 0.59822 7.68960 -841.70186 0.13066 -0.02770
6 0.59822 7.69122 -841.70186 0.12861 -0.05730
7 1.42741 6.65331 -841.69463 1.61149 -3.67174
7 1.42741 6.90839 -841.69463 1.64698 -0.48946
7 1.42741 6.78707 -841.69463 1.66903 0.20359
7 1.42741 6.91323 -841.69463 1.66389 -0.05350

Run: No2
Step Force Torque Energy Curvature Angle
1 2.61504 8.31185 -841.65323 3.46095 7.17148
1 2.61504 8.70025 -841.65323 2.47308 8.93185
1 2.61504 7.11221 -841.65323 1.69398 10.05366
1 2.61504 6.41413 -841.65323 1.25154 4.12855
2 2.96452 5.93151 -841.58078 0.88042 6.69672
2 2.96452 8.39693 -841.58078 0.77810 8.35420
2 2.96452 6.38603 -841.58078 0.34009 6.42709
2 2.96452 4.94667 -841.58078 -0.08053 4.34487
3 2.46190 4.58448 -841.64906 0.77501 9.79908
3 2.46190 5.80668 -841.64906 0.48672 5.39046
3 2.46190 4.75331 -841.64906 0.02254 2.77639
3 2.46190 3.52371 -841.64906 -0.17823 3.84788
4 1.76334 6.41776 -841.67902 1.23938 6.89719
4 1.76334 8.48847 -841.67902 1.26740 5.00503
4 1.76334 3.40856 -841.67902 0.99245 6.13151
4 1.76334 4.73352 -841.67902 0.94250 3.44075
5 1.80469 7.22108 -841.61309 -0.32111 12.30048
5 1.80469 8.14264 -841.61309 -0.61186 6.11802
5 1.80469 4.92011 -841.61309 -0.99815 6.81515
5 1.80469 5.46823 -841.61309 -1.11281 3.98679
6 1.97160 5.78832 -841.64270 -1.06581 14.36118
6 1.97160 7.07648 -841.64270 -1.36489 3.33882
6 1.97160 3.55540 -841.64270 -1.76720 4.23498
6 1.97160 3.35267 -841.64270 -1.91090 2.67489
7 1.62263 5.98730 -841.66725 -0.67812 13.81606
7 1.62263 8.25427 -841.66725 -0.70925 4.14955
7 1.62263 5.89056 -841.66725 -0.75888 4.81049
7 1.62263 5.30416 -841.66725 -1.36576 3.60098
8 0.78107 7.76665 -841.68870 -0.14496 13.34619
8 0.78107 7.38683 -841.68870 -0.35105 5.89227
8 0.78107 4.71578 -841.68870 -0.55597 5.69780
8 0.78107 4.57943 -841.68870 -0.69759 4.46563
9 0.66283 7.49425 -841.69482 -0.85882 15.97305
9 0.66283 9.51917 -841.69482 -0.66620 4.15247
9 0.66283 6.21579 -841.69482 -0.83140 3.78268
9 0.66283 6.59462 -841.69482 -0.99738 6.89302

Run: No3
Step Force Torque Energy Curvature Angle
1 0.72993 9.94363 -841.69947 0.71889 10.94416
1 0.72993 12.92972 -841.69947 -0.04371 11.57521
1 0.72993 7.98453 -841.69947 -0.36732 2.19454
1 0.72993 6.73489 -841.69947 -0.47136 2.18831
2 0.40371 9.89070 -841.70129 0.40813 3.10282
2 0.40371 10.17441 -841.70129 0.32200 4.73522
2 0.40371 11.87498 -841.70129 0.33368 5.83648
2 0.40371 12.41256 -841.70129 0.44366 1.06109
3 1.61062 12.60386 -841.68798 2.75647 2.70587
3 1.61062 11.79408 -841.68798 2.73943 3.05205
3 1.61062 10.48933 -841.68798 2.61704 2.83768
3 1.61062 9.44270 -841.68798 2.45680 2.43934
4 2.84345 8.35024 -841.61147 3.90909 2.66987
4 2.84345 8.06435 -841.61147 3.85312 4.63246
4 2.84345 8.97426 -841.61147 3.48812 7.91011
4 2.84345 10.32518 -841.61147 3.36929 12.80998
5 2.92773 11.05383 -841.46629 3.08306 10.97800
5 2.92773 10.13407 -841.46629 2.92012 11.14227
5 2.92773 10.42716 -841.46629 2.88057 18.43009
5 2.92773 11.96853 -841.46629 2.94386 43.11650
6 3.93549 14.27979 -841.28190 4.56058 57.13172
6 3.93549 11.93295 -841.28190 3.97326 61.48398
6 3.93549 9.21658 -841.28190 3.26823 60.92482
6 3.93549 7.33502 -841.28190 2.59981 60.13447
7 4.23884 7.82969 -841.03882 2.91808 58.72140
7 4.23884 5.20232 -841.03882 1.25080 28.80764
7 4.23884 8.82547 -841.03882 0.57005 6.04428
7 4.23884 4.81322 -841.03882 0.03926 3.60428

Still proceeding...

These dimer calculations do not look good.

First, you can always set ALGO, LREAL, NSIM, NPAR ... to improve efficiency.

The dimer calculations are not near a saddle, and they are not obviously converging towards one. Something doesn't look quite right here. Do you have ediff set to a low value (around 1e-8?). If you post your INCAR file, we might see a problem or suggest different settings.

As a check, try starting with your initial structure and increase the value of DRotMax to 10, for example. Make sure that the Torque (3rd column) drops to a smallish value and the Curvature (5th column) also drops. If you are near a saddle, the curvature should become negative and stay negative. Then, as it translates, the force (2nd column) should drop as you converge to a saddle.

If you have already done a rough NEB, you can use the neb2dim.pl script to start a dimer run fairly near to the saddle.

Overall, I would suggest using a conservative optimizer (e.g. IBRION=3, POTIM=0.2) and just make sure that you can systematically move towards a saddle (assuming you are starting near one). If your starting point is not near a saddle, you will see the curvature become positive at times and the configuration will climb up the potential in search of a saddle. This can take a while and it not recommended for an expensive system.

Thank you for the reply,

I will try your suggestion of increasing value of DRotMax to 10 and see if that makes any difference.

---Here is the INCAR File that I am running the Dimer calculations with---
System=Dimer

ISMEAR=0
SIGMA=0.05
PREC=Medium
ISPIN=2
MAXMIX=80
NFREE=10
ENCUT=275.0
NSW=30
EDIFF=1E-7
EDIFFG=-0.01

# Disabling standard VASP optimisers
IBRION=3
POTIM=0

# Dimer Method
IOPT=2
ICHAIN=2
Ddr=5E-3
DRotMax=4
DFNMax=1
DFNMin=0.01

LREAL=A
NSIM=4
ALGO=FAST

Some more information that may help:

If a calculation reaches the end of the NSW cycles I have specified, I add in ICHARG=1 and ISTART=1 in the INCAR file and overtake the older POSCAR and MODECAR by CONTCAR and NEWMODECAR respectively.
I am using the USPP_GGA standard VASP potentials and 1-KPoint (located in the centre of my reciprocal space, ideally more once the structure converges).

Method employed:
1) Built a 4-image NEB, ran with the standard NEB setup described on the VASP TST Tools website (LGLOBAL=.TRUE. and LCLIMB=.False.). Until the geometries were sufficiently corrected (to fix some warping of the structures due to a required rotation, coming from the linear extrapolation method to build the images)
2) Extended the 4-image run to 8-images, extrapolating between the steps already optermised
3) Ran until a (~60-80 ionic steps) crude spline had been mapped out that showed a definite saddle.
4) Produced a dimer run, using the neb2dim.pl for this NEB result and ran the calculations (This was what was shown in the original post results)
5) While the dimer method was proceeding, enabled the climbing method (LCLIMB=.TRUE.) for the NEB and allowed it to continue to run
6) NEB method pushed the high energy structure to the top of the spline, but seemed to 'dance' around the top of the spline, not converging
7) Created another dimer run, using a snapshot of the climbing method NEB result at the top of the spline as input. Produced a MODECAR file between the second highest and highest energy structure.
This seemingly also had troubles (though it has run for less time)

As you can see I have tried a few different approaches to finding the TS. Perhaps some of this additional explanation into what I have done will bring the problem to light.

I you send us the dimer and/or NEB calculations, we'll take a look and try to see what's going on.

Hello Professor Henkelman,

I have modified slightly what I am doing with these calculations and will report if this makes a difference.

I have reduced my unit cell from the original 181 atoms to 110. This has reduced the isolation of the molecule on silicon I am looking at but hopefully still within acceptable limits. I was probably too generous with the slab size before and perhaps the size was to large to ever give the dimer method any chance of converging. Indeed after increasing the DRotMax to 10, there was a definite improvement in the apparent trend to convergence but the additional CPU time required was significant.

Thank you for you offer to look over my files but at this stage I am happy to continue to play with them, if I have no luck with this second run then perhaps I will take you up on your offer.