UT theoretical chemistry code forum

Hello,

I am comparing the barrier energy for a simple system (Fe adatom migrating between hollow sites on a graphene sheet) calculated with two methods:

1. Running it normally with the CINEB method with 3x3x1 grid points (SUCCESSFUL)

2. Running a rough CINEB method with 1x1x1 grid points, then running only the saddle point image with the Dimer Method with 3x3x1 grid points (NOT SUCCESSFUL)

Unfortunately, the dimer method is taking the same, if not longer, to converge as the 3x3x1 CINEB run which defeats the point.

Questions:
1. I have been running the Dimer Method with 2 nodes, since I assumed that the two images run in parallel similarly to the CINEB method images. Is this accurate?

2. I've been initializing the Dimer Method with the exact CINEB saddle point image (no interpolation, not using neb2dim.pl) and then making a MODECAR pointing to the next image in the reaction pathway from the CINEB run with modemake.pl. Should I be setting this up differently?

3. Would it be better to run a CINEB with a few images with 3x3x1 kpoints with a conservative optimizer and then going to the Dimer Method?

I also ran the Dimer Method from the saddle point with 1x1x1 kpoints which converged very quickly (not surprising) and used the MODECAR generated from that to continue the run with 3x3x1 k points and still had a very long time.

I have attached both the initial CINEB 1x1x1 run and the Dimer Method attempts. Maybe I'm missing something else? Note: I killed the 2nd Dimer Method run after the DIMCAR forces started increasing consistently (round2 dir in Fe-hopping-dimer 3x3x1 run).

Thank you.

I took a look at your system and I think there are two complicating factors, neither of which really has to do with k-point sampling.

First, I was able to reproduce your NEB calculation (see neba.tar.gz). I used a somewhat reduced z-axis to save some time. It has a barrier of about 0.9 eV and shows a fairly straight diffusion path between the hollow sites over a C-C bridge.

As I was doing this calculation, I noticed some strange energetics both in things I tried and as compared to your calculation (as well as high forces in dimer calculations). I realized that my calculations would often get stuck in local electronic minima. You can see in your calculation that the spin state on the Fe atom changes along the reaction coordinate between 0 and 2 and 4. This is the first complicating factor.

To help simplify the problem, I did a closed shell calculation with ISPIN=1 (nebb.tar.gz). I let it run for a while and you can see that the energy barrier is lower, at about 0.8 eV. Also, if you look at the path, it breaks symmetry and the Fe atom diffuses to a C top site at the transition state, rather than a C-C bridge. This is the second complicating factor. It is difficult for a calculation to spontaneously break symmetry, but I believe that the bridge site is not a first-order saddle and upon relaxation, the transition state becomes the top site.

To make matters worse, there is a 3-fold symmetry at the top site and so this is a so-called monkey saddle. I'm not sure how the dimer method handles monkey saddles - my expectation is that it will get close to the saddle but may have a hard time actually converging - but I'm not sure.

So then I did some tests with the dimer method. In dima.tar.gz, which was started from the NEB calculation over the bridge site, the dimer eventually breaks symmetry and relaxes away from the bridge to the top, consistent with that being a lower energy true saddle as in nebb. Second, I tried turning on spin polarization in dimb,tar,gz and you can see that the forces and curvatures go crazy, and the spin state changes during the calculation.

So while I'm not providing completely converged calculations, I think I understand what's going on. Converging to specific spin states in vasp can be a pain - for that you might consider the occupation matrix control provided by Graeme Watson. Or, you can converge to the proper (monkey) saddle at the top site and then make sure that you have the lowest spin state, either with the dimer or neb methods.

Good luck - you managed to pick a challenging system here!