UT theoretical chemistry code forum

Dear All,

I have been attempting to use the dimer method implemented in VTST to identify transition states, as I found obtaining converged CI-NEB calculations to be rather more resource intensive, and not always reliable. Having initialised dimer calculations from partially converged CI-NEB calculations using the VTST script neb2dim.pl, I have run dimer calculations using a single gamma-centred k-point to obtain reasonably good TS geometries before restarting with a finer k-point mesh. Some of my gamma point dimer calculations have converged, others not, but I have a few questions that I would be most appreciative of any advice on:

1. In the cases where my gamma point dimer calculations converged, I noticed that the force reported in the DIMCAR isn't actually lower than the EDIFFG specified in the INCAR. For example, here's the last few lines of the DIMCAR for a dimer calculation that terminated correctly, when I had EDIFFG = -0.01 eVA-1:

7 0.19111 3.61373 -1354.29673 -9.32915 2.90197
7 0.19111 2.60226 -1354.29673 -9.48132 1.84870
7 0.19111 1.68391 -1354.29673 -9.57579 2.64930
7 0.19111 2.14196 -1354.29673 -9.64715 1.03717
8 0.13445 0.88171 -1354.29742 -9.69178 3.93975
9 0.08895 0.57143 -1354.29754 -9.72473 0.22764
10 0.05604 0.13000 -1354.29761 -9.69224 0.18518
11 0.05135 0.27221 -1354.29764 -9.68289 0.24851
12 0.04810 0.24867 -1354.29767 -9.68167 0.11577
13 0.04047 0.09904 -1354.29772 -9.64471 0.23648
14 0.05957 0.08409 -1354.29775 -9.63287 0.04620
15 0.04578 0.08291 -1354.29781 -9.63245 0.08409
16 0.04088 0.05698 -1354.29783 -9.63662 0.06725
17 0.04178 0.06975 -1354.29785 -9.64293 0.05002
18 0.03693 0.05630 -1354.29787 -9.65126 0.04628
19 0.03844 0.05533 -1354.29790 -9.67250 0.04773
20 0.03877 0.11386 -1354.29791 -9.65496 0.11680
21 0.04498 0.12808 -1354.29794 -9.68010 0.07782
22 0.03194 --- -1354.29797 --- ---

So, all of this looks promising - the curvature is consistently negative, the angle has decreased to a small number, and the forces are small, but I'm not sure why the calculation has finished here when the lowest force achieved, 0.03194 evA-1, is still bigger then my EDIFFG. I'm assuming I'm missing something here, so any advice would be most appreciated.

2. When restarting converged gamma point dimer calculations with a finer k-point mesh, I initialised the calculation by using the CENTCAR and NEWMODECAR from the pre-converged gamma point calculation as the POSCAR and MODECAR, respectively, for the calculation with a finer k-point sampling mesh, obviously I used the same POTCAR and mostly the same INCAR (just adjusting some parallelisation parameters). I had hoped that having obtained a pre-converged TS geometry and a (hopefully) accurate dimer direction that refining at the finer k-point sampling would be rather trivial, unfortunately so far the calculation doesn't appear to be going great, as the DIMCAR shows that the curvature has gone positive in the second dimer rotation loop:

Step Force Torque Energy Curvature Angle
1 1.05227 36.44992 -1353.72307 -9.29018 57.26713
1 1.05227 19.09210 -1353.72307 -26.58704 16.20399
1 1.05227 15.85489 -1353.72307 -29.15829 8.41732
1 1.05227 8.31876 -1353.72307 -30.06472 7.17612
2 0.72959 7.62066 -1354.07242 6.44210 21.89512
2 0.72959 8.78510 -1354.07242 3.34589 7.84344
2 0.72959 5.92504 -1354.07242 2.05094 21.94122

I imagine it's too early to tell (the calculation is still running, unsurprisingly much more slowly using the finer k-point sampling mesh), but I would be appreciative if any more experienced users can tell me if this is something to worry about and if the calculation is unlikely to converge. I can see the force has decreased from step 1 to step 2, so I intend to allow this calculation to continue to run to at least wait and see.

3. Some of my gamma point dimer calculations haven't been going as well - I presume simply because the initial dimer geometry and dimer mode from my partially converged CI-NEB wasn't very accurate. I've noticed the energy going up quite considerably in later steps compared to earlier ones, which is worrying as obviously this would give insanely high activation barriers for the process I am looking at:

Step Force Torque Energy Curvature Angle
1 2.06494 44.14745 -1347.04542 -8.01006 58.74842
1 2.06494 22.68473 -1347.04542 -32.59916 21.02689
1 2.06494 19.57650 -1347.04542 -36.12908 9.07447
1 2.06494 7.88115 -1347.04542 -37.96348 7.57726
2 1.78922 6.94262 -1347.04974 5.12735 22.43672
2 1.78922 7.75164 -1347.04974 2.27307 6.11048
2 1.78922 3.02646 -1347.04974 1.34965 10.65861
2 1.78922 3.58664 -1347.04974 0.78190 2.32461
.
.
.
.
37 39.61517 45.08142 -1338.33328 14.23220 -0.24989
37 39.61517 46.39824 -1338.33328 14.84377 -0.99662
37 39.61517 51.94720 -1338.33328 15.11193 -0.49949
37 39.61517 53.96301 -1338.33328 14.00385 1.88844
38 46.14136 99.54852 -1335.92757 14.54611 -42.28752
38 46.14136 125.53077 -1335.92757 -3.61029 -22.23385
38 46.14136 139.97238 -1335.92757 -9.81092 -6.03310
38 46.14136 157.78971 -1335.92757 -8.70032 2.29629

So that's a 10 eV difference, so clearly something is wrong. Fortunately I haven't wasted too much resources on this calculation since it has only been running overnight at the gamma-point, but I would be most appreciative of any advice on why this might be the case. Clearly, the energy increase is problematic and the forces are very large. The geometry for the dimer doesn't look crazy, though. I have set EDIFF = 1E-7 as recommended, so if this threshold is being met then the forces should be accurate, as suggested in the VTST documentation.

Thanks in advance for any comments or advice, and apologies if these are basic questions, I am relatively new to using the VTST dimer method.

These are all good questions - most of which I should be able to address. If you are willing to attach a .tar.gz of the three calculations (without CHG* / WAV*) then I can take a closer look.

Thanks Graeme, and apologies for the slow response - I was investigating some of my calculations. In some cases, the dimer calculation restarting from the finer k-point mesh did converge eventually. In some other cases, convergence failed because SCF convergence wasn't being achieved, hence inaccurate forces prevented the dimer from converging - in these cases, playing with the SCF algorithm, or adding empty bands, worked. However, I am still running into some problems where dimer convergence is failing, despite it having converged at the gamma point, and without there being any evident problems with SCF convergence.

Here's the DIMCAR output for the calculation that did converge at the gamma point:

Step Force Torque Energy Curvature Angle
1 0.25276 50.46858 -1357.13203 -9.56117 56.63579
1 0.25276 24.91699 -1357.13203 -31.54181 15.75328
1 0.25276 20.99383 -1357.13203 -34.12902 7.86696
1 0.25276 10.96020 -1357.13203 -35.15687 8.53903
2 0.14438 10.22297 -1357.88731 12.74356 12.46715
2 0.14438 13.63767 -1357.88731 10.48603 10.17775
2 0.14438 12.18330 -1357.88731 8.20194 10.31986
2 0.14438 12.76450 -1357.88731 5.73246 14.12355
3 0.49031 15.69399 -1357.88705 2.67437 10.92633
3 0.49031 11.38768 -1357.88705 -0.78915 26.07551
3 0.49031 17.77019 -1357.88705 -6.04654 6.96991
3 0.49031 5.94998 -1357.88705 -8.60411 12.29095
4 0.15429 11.86589 -1357.88769 -9.76489 4.65704
4 0.15429 4.53072 -1357.88769 -10.51881 1.86455
4 0.15429 2.67276 -1357.88769 -10.69167 3.84966
4 0.15429 4.20630 -1357.88769 -10.86025 1.14996
5 0.16332 1.42355 -1357.88740 -10.77613 2.94825
5 0.16332 2.40878 -1357.88740 -10.85177 0.93167
5 0.16332 0.96807 -1357.88740 -10.88279 0.38177
6 0.10724 0.60367 -1357.88742 -10.81793 0.89451
7 0.06226 1.02518 -1357.88747 -10.77543 0.27164
7 0.06226 0.34628 -1357.88747 -10.78017 0.70572
8 0.04553 0.57347 -1357.88749 -10.77134 0.23420
9 0.02086 --- -1357.88751 --- ---

And here's the output (so far) for the same system at a finer kpoint mesh (4x4x1):

Step Force Torque Energy Curvature Angle
1 7.67439 22.69155 -1354.90488 3.35225 19.50386
1 7.67439 51.56921 -1354.90488 4.11866 57.99481
1 7.67439 31.54808 -1354.90488 -2.90042 15.08302
1 7.67439 10.52776 -1354.90488 -5.27482 4.62245
2 9.62875 106.33004 -1356.89424 28.39508 22.96367
2 9.62875 27.96914 -1356.89424 14.90239 1.31204
2 9.62875 33.09297 -1356.89424 14.87617 0.41910
2 9.62875 34.78764 -1356.89424 14.87685 0.18581
3 22.63407 46.35925 -1356.24462 16.60335 1.83328
3 22.63407 21.49693 -1356.24462 15.10064 1.28660
3 22.63407 34.59352 -1356.24462 14.82734 5.52310
3 22.63407 60.90094 -1356.24462 13.84422 2.15100
4 12.39779 54.63559 -1355.19530 15.20024 0.11865
4 12.39779 54.86826 -1355.19530 15.19596 0.14343
4 12.39779 55.14321 -1355.19530 15.19924 0.01521
4 12.39779 55.17398 -1355.19530 15.20409 0.00742
5 12.87861 23.18660 -1353.49969 15.55081 14.63610
5 12.87861 30.36422 -1353.49969 14.39694 -0.84253
5 12.87861 29.90493 -1353.49969 14.41360 0.13028
5 12.87861 29.98038 -1353.49969 14.40053 -0.02814
6 16.34900 15.36078 -1351.44857 12.40958 25.33441
6 16.34900 10.87040 -1351.44857 9.50114 -0.00857
6 16.34900 10.87131 -1351.44857 9.50771 0.09618
6 16.34900 10.87518 -1351.44857 9.51445 0.08158

So, clearly this isn't going very well, the forces are large, the curvature is positive, and the total energy is getting pushed up by over 5 eV. Here are my INCAR settings, ignore the commented out lines as these pertain the the partially converged NEB calculation I did to initialise the dimer:

#ICHARG = 1
ENCUT = 650
LREAL = .FALSE.
ISMEAR = 1
SIGMA = 0.2
GGA = RP
NELMIN = 8
NELMDL = -5
#EDIFF = 1E-5
EDIFFG = -0.01
NSW = 400
ISIF = 2
#IBRION = 3
#POTIM = 0.01
LCHARG = .FALSE.
LWAVE = .FALSE.
ISPIN = 2
IVDW = 12
NCORE = 16
KPAR = 8
#NPAR = 32
ISYM = 0
ALGO = FAST
MAXMIX = 40
ADDGRID = .TRUE.
LORBIT = 11

#ICHAIN = 0
#IMAGES = 8
#SPRING = -5.0
#LCLIMB = .TRUE.


ICHAIN = 2
IOPT = 2
IBRION = 3
POTIM = 0.0
EDIFF = 1E-7

Let me know if you have any suggestions. I was also discussing with a colleague who suggested that the TS geometry obtained from the fully converged gamma-point calculations should be accurate, and therefore it might be unnecessary to restart the dimer with the finer k-point mesh, and that I could just run a single-point calculation with the finer k-point mesh to get an accurate activation barrier. Does this sound sensible? I can see the logic, but I'd like to be sure before pursuing that avenue, which could save me some work, but I wouldn't want to do this if I were not able to have a high level of confidence in the activation barriers obtained using this approach.

I could answer the first question. You should look at the value of FORCES: max atom, RMS in OUTCAR and when this value is below EDIFFG dimer run is stopped. This post may help http://146.6.145.114/forum/viewtopic.php?t=12193.
Hope some advice may be given by other experts about questions on convergence problems when using fine K-points.