UT theoretical chemistry code forum

Hello all,

I am new to using VASP and your tools so please excuse my possibly rudimentary questions.

I am investigating the dissociation of a moderately sized organic molecule on a surface, and attempting to calculate a activation barrier using the DIMER method. I have not run a NEB on the system, just for clarification. To initialize the DIMER method I have used the nebmake.pl script to generate an interpolation between initial and final states, which were optimized. Using two images near what I think might be the transition state I use neb2dim.pl to get a MODECAR, and use the image I think is close to the transition state as image zero, i.e., the geometry used to start the DIMER calculation. This geometry consists of the molecule with one stretched bond. When running the calculation it appears to go in the right direction but the molecule appears to relax back to geometry similar to the initial state, but slightly perturbed and interacting with the surface more. So to make a short story long :) my questions:

1) For a complex molecule, say 10-20 atoms, how close should DIMER starting geometry be to the transition state for the method to work.

2) In generating the MODECAR file, is it better to use two geometries far away from each other, i.e., like say a transition state guess and the optimized molecule on the surface, or is it better to use two images close together on the reaction pathway, like two sequential images of a 7 image NEB run or interpolation.

3) Is it better to use a starting geometry for the DIMER that is on the down slope of the reaction pathway or the upslope? It seems when working with a large molecule it may be more beneficial to clearly indicate the bond being broken by using an image on the down slope of the reaction. I hope you know what I mean :)

My INCAR file: (could be wrong, I just started working with VASP)
NWRITE = 2
ISTART = 1
ISPIN = 1
GGA = 91
NPAR = 8
LMAXMIX = 6
ENCUT = 400.0
EDIFF = 1E-7
LREAL = Auto
NELM = 60
NELMIN = 6
EDIFFG = -0.001
IBRION = 3
NSW = 600
NBLOCK = 1
ALGO = FAST
ISMEAR = -5
SIGMA = 0.05
POTIM = 0.0
IDIPOL = 3
LDIPOL = .TRUE.
LDAU = .TRUE.
LDAUTYPE = 1
LDAUL = -1 3 -1 -1
LDAUU = 0 5 0 0
LDAUPRINT = 0
ICHAIN = 2
DdR = 0.005
DRotMax = 4 !// previous value = 1
DFNMin = 0.01
DFNMax = 1.0
IOPT = 2

The last few lines of my DIMCAR file:
30 0.51032 3.32724 -351.27833 -4.42478 1.20466
30 0.51032 2.96650 -351.27833 -4.29105 4.33231
30 0.51032 3.18411 -351.27833 -4.53838 1.13574
30 0.51032 3.43508 -351.27833 -4.54371 1.71452
31 0.44283 3.67865 -351.28294 -5.28193 -1.28728
31 0.44283 3.61000 -351.28294 -5.19690 0.75463
31 0.44283 3.61762 -351.28294 -5.19830 0.07899
31 0.44283 3.62484 -351.28294 -5.19830 0.01892
32 0.52486 4.75545 -351.28736 -2.94184 1.75873
32 0.52486 4.83132 -351.28736 -2.86920 0.90141
32 0.52486 4.89768 -351.28736 -2.86104 0.17169
32 0.52486 4.91327 -351.28736 -2.85899 0.03173
33 0.58855 4.30602 -351.29194 -4.26492 -2.84982
33 0.58855 3.93239 -351.29194 -4.22991 3.55128
33 0.58855 2.22569 -351.29194 -4.39157 1.34700
33 0.58855 2.80753 -351.29194 -4.37659 1.09530
34 0.44651 3.19381 -351.29711 -4.15048 -1.40087
34 0.44651 2.87161 -351.29711 -3.99770 0.57229
34 0.44651 2.98279 -351.29711 -3.99510 0.19831
34 0.44651 3.02503 -351.29711 -3.99032 0.10115
35 0.38314 3.29311 -351.29918 -3.05737 0.34689
35 0.38314 3.41736 -351.29918 -2.87860 0.46710
35 0.38314 3.50658 -351.29918 -2.87618 0.01487
35 0.38314 3.50953 -351.29918 -2.87424 0.00593
36 0.51575 3.62416 -351.30179 -1.98746 -2.88825
36 0.51575 3.33007 -351.30179 -1.88662 -0.20677
36 0.51575 3.30386 -351.30179 -1.89432 -0.16201
36 0.51575 3.30402 -351.30179 -1.89303 -0.03732
37 0.45975 3.01846 -351.30486 -3.83490 2.10664
37 0.45975 3.25606 -351.30486 -3.67007 0.52257
37 0.45975 3.35326 -351.30486 -3.66895 0.03006
37 0.45975 3.36026 -351.30486 -3.66753 0.00411

The calculation has been running for a while now.

Any help would be much appreciated.

Thanks
Siris

It sounds like you may not have a very accurate guess at the transition state to start your dimer calculation. This is not necessarily bad, but it will take longer to converge, and may not converge to the saddle point that you expect.

If you initiated the dimer using a linear NEB, the MODECAR will just be a vector between reactants and products. It will not matter which images you use to generate the MODECAR. It would be better to run a couple of iterations of the NEB (if you can afford it) until the forces are <1 eV/Ang. Then the neb2dim script can make a sensible guess at the transition state, based upon the energy of the NEB images. The MODECAR will then also be more accurate since it is interpolated to be at the point of maximum energy along the path.

1) is difficult to answer in general. There is a radius of convergence around each saddle point, but the distance will be system specific.

2) It would be best to use images as close to the estimated saddle point as possible.

3) It doesn't matter if you are higher or lower in energy than the saddle when converging, as long as there is a negative mode in the valley/ridge of the saddle. That said, there is only one negative mode going up the valley to the saddle and many other modes along the ridge from the saddle. So in practice, you almost always relax down in energy to the saddle point.

Fortunately, in your example, there is a stable negative mode. Also the mode does not change much during the rotation steps, and the dimer angle also does not change. This is bit strange given the high torque. You might want to make sure there are no problems caused by DFT+U. We don't have much experience running the dimer on these energy surfaces, but I imagine that they are somewhat rougher than standard DFT.

Follow up:

The problem calculation I was describing earlier has been worked on further. It turns out that there were two transition states very close to each other in geometry. One was not what we were expecting, thus was deemed an anomaly, however, was in fact a transition state. Interestingly, using images at the top of a finished NEB calculation to start DIMER runs, both of the transition states could be converged upon. However, the transition state with the lower barrier and possibly a more "flat" PES near the MEP maximum was less favored by the DIMER calculation.

All in all, the lesson learned is to use a transition state geometry as close to the real transition state geometry as possible when dealing with complex reaction pathways, and assume anything is possible :)

Also the DIMER method appears to work just fine in calculations using the +U method in VASP.