UT theoretical chemistry code forum

Hi all,
I found a transition state with NEB and then optimized the structure with IBRION=1. I used EDIFFG=-0.03. I did not want to use even stricter force criteria, since the optimization needed more than 2000 geometry steps to converge.
I decided to run a Full Hessian calculation (120 atoms = 720 energy and force calculations with IBRION-5, EDIFF=1E-6, POTIM=0.015, NFREE=2).
Since the calculation took so much time, I decided to run a partial Hessian calculation in two ways:
i) By using the Dynamical matrix script. Displacement was chosen again to be 0.015 and 11 atoms were chosen (NSW = 11*3+1=34 )
ii) with IBRION=5 and running selective dynamics and choosing as active the same 11 atoms and freezing all other. Again the same settings as in DYNMAT were used: (same SIGMA, ISMEAR and EDIFF=1E-6, POTIM=0.015) and NFREE=2 .

The DYNMAT results were not agreeing very well with IBRION=5:
Full Hessian : 1 imaginary frequency at 241 cm-1 (plus of course low freqs between -10 and -1 cm-1 that correspond to translations)
Partial Hessian: 1 imaginary frequency at 216 cm-1
DYNMAT: 2 imaginary frequencies at 221 cm-1 and 61 cm-1.
The largest imaginary frequency corresponds to the breaking of a bond.
The other freq (-61 cm-1) corresponds to a rotation of a -OCH3 group.

My question is about the accuracy of the DYNMAT method.
If I have understood correctly, there is only one displacement in the DYNMAT method. In the IBRION=5 and NFREE=2, each atom is displaced in each direction by a positive and negative displacement.
However in DYNMAT only one displacement is chosen. So, the method is similar as setting NFREE = 1 . Am I too wrong ?
This setting is not suggested according the vasp manual.
So, I think that these discrepancies are caused by the fact that only one displacement is chosen in DYNMAT. Do you agree ?
Has anyone of you taken similar results?
Are you thinking of implementing a double displacement for the DYNMAT method ?

DYNMAT is really helpful for big systems, but because of these disagreements unfortunately I can not use it for my system.

Could you please make a comment or suggest anything ?

Thank you,
Andreas

I do not agree that central difference is necessary or optimal for normal mode calculations. In our tests, it is essential to find convergence with respect to accurate forces (small ediff) and a small displacement (POTIM, using IBRION=-5). Once you have found convergence, your modes should not be different if you have used forward or central difference. In my experience, central difference can not save you if your forces are not accurate enough or your displacements are too large. Looking at your parameters, both values appear too large to me. Try EDIFF between 1e-7 and 1e-8 and POTIM (or displacement) between 0.001 and 0.005. It is the large displacement that worries me the most.

You can run the DYNMAT code using central difference. This is done by doing a dynmat calculation with positive displacements in your DISPLACECAR file, and then again with negative displacements. The dymfit.pl script (see http://theory.cm.utexas.edu/vtsttools/scripts/ ) can be used to combine the resulting Hessian matrices into a single one. A linear fit using two matrices corresponds to the NFREE=2 central difference scheme, but you can use more displacements and a higher order fit if you want. But again, I do not recommend doing this -- a single forward difference calculation that is converged is more efficient than trying to correct for errors in the force or large displacements using more than one displacement.

Also, some sanity checks would be good. Look at your Hessian matrices to make sure that there are no zeros in them. If you didn't specify NSW correctly, or your job terminated before doing all the iterations, you can get zeros in the matrix and modes that do not make sense.

Dear Graeme,
thank you for your suggestions.

The calculation with displacement=0.001 and EDIFF=1e-7 still gives one additionally negative frequency.
I am still checking the other settings: EDIFF=1E-8 or negative displacement and also run the DYNMAT with a double displacement.

If you do not see the modes change between small displacement values and ediff values, I would trust them and aim for a tighter convergence to the saddle. You can use the neb2dim.pl script to run a less expensive dimer calculation from your NEB. Also, if the path is fairly linear, you might get away with 3 or just 1 image along the NEB. Also, if you are seeing very slow convergence, you might want to try some of our second order NEB optimizers, such as lbfgs.