Two imaginary frequencies in the TS obtained from dimer

[spbatchu]

Hello to the admin team

I ran a dimer calculation to get the transition state(TS) of H-abstraction from beta-C of ring-opened intermediate of tetrahydrofuran on ZrO2. The forces in the dimer converge and the curvature is negative. I went ahead to do the vibrational frequency calculation of the structure(TS) obtained from the dimer and I find that there are two huge imaginary frequencies (1556.9i cm-1 and 567.0i cm-1) and I was initially expecting only 1 imaginary frequency.

I am wondering if there is a way to make changes to the calculation to get the right transition state from what I have. It will be of great help if you can provide suggestions. Thanks in advance

P.S: I have attached the dimer-related input and output files and the files for vibrations( I use ASE python to do the vibrational frequencies, hence I attach the relevant files).

[graeme]

The dimer calculation looks ok to me. The vibrational calculation, however, seems to be using displacements that are probably too large. I like to use accurate forces (ediff = 1e-7) and displacements on the order of 0.002 Ang. It looks like you are using something more like 0.03 Ang, which is given large forces that are almost certainly outside of the harmonic region of the potential. Do consider using vasp's built-in perturbation theory method, which does not rely on a finite difference step size.

[spbatchu]

Hello Dr. Graeme

Thank you so much for your reply. I will rerun the calculation using accurate forces (ediff=1e-7) and check if it helps. I use a displacement of 0.015 Angstrom for vibrational frequency calculation. I will try lowering this value to 0.002 Ang and get back to you. I will also try out the vasp built-in perturbation theory method.

I have another question, Dr. Graeme:
When I look at my DIMCAR file again, I see that the Torque doesn't drop consistently(DIMCAR text below). Do you think this could be causing the problem or is it okay to have such torque values?


Step Force Torque Energy Curvature Angle
1 0.08936 1.03986 -1802.42042 -14.43426 0.78777
1 0.08936 0.71384 -1802.42042 -14.40941 1.01516
2 0.06593 1.12903 -1802.42051 -14.39134 0.75013
2 0.06593 0.52979 -1802.42051 -14.46150 0.57040
3 0.06554 0.97011 -1802.42062 -14.43915 1.34459
4 0.07708 0.91831 -1802.42071 -14.48730 0.49409
5 0.13098 1.00762 -1802.42097 -14.57211 0.56570
5 0.13098 1.31426 -1802.42097 -14.57984 0.99045
6 0.12854 1.11725 -1802.42129 -14.65174 0.10854
6 0.12854 1.17772 -1802.42129 -14.78828 0.47180
7 0.08963 1.10372 -1802.42163 -14.60022 0.59287
7 0.08963 1.05889 -1802.42163 -14.68404 0.98635
8 0.08963 --- -1802.42163 --- ---

[graeme]

The key measure of convergence for the dimer calculation is the force, and the fact that the curvature is negative. Both are true for you. I don't see any problem with the dimer calculation. We'll see, but I'm pretty sure that a calculation of curvatures closer to the minimum, or using perturbation theory, will solve the problem of 2 negative modes.

[spbatchu]

Thank you, Dr. Graeme

I have tried using perturbation theory and VASP throws this warning saying that perturbation theory cannot be used for DFT+D calculations. I currently use VASP 5.4.1 and use the PBE+D3 setting for my calculations. I am not sure whether there is a way to make this perturbation theory calculation work and kindly suggest me in case there is one to your knowledge.

[graeme]

Ah, ok, then a finite difference calculation with more accurate forces and smaller displacements should be fine.

[spbatchu]

Hello Dr. Henkelman

I lowered the displacement to 0.002 Angstrom in the above-discussed example and I see that the 2nd imaginary frequency disappeared. Thanks for the suggestion.

I have tried the same idea of lowering the displacement to 0.002 Angstrom in other dimer calculations(which also have multiple imaginary freq.) and I am unable to remove the additional imaginary frequencies. These dimer calculations are attached below. A quick summary of these attachments are as follows:

1) There are 3 transition states -TS1, TS2, and TS3.
2) TS1 and TS2 have been obtained using NEB+dimer whereas TS3 has been obtained using NEB+climb.
3) For TS1, I use ediff of e-06 for dimer and potim of 0.002 Ang for vib. freq. For TS 2, I try both ediff of 10-6 and 10-8 with potim of 0.002 Ang
4) For TS3, I use ediff of e-06 and ediffg of -0.02 eV/Ang for the NEB+climb and potim of 0.002 Ang for vib freq.
5) To my understanding, all the dimer and climb calculations are converged.

And, I get the imaginary vibrational frequencies as follows:
TS1_dimer_e-06
92.496734i cm-1
102.668960i cm-1
296.614077i cm-1

TS2_dimer_e-06
76.546732i cm-1
89.649038i cm-1
286.467675i cm-1

TS2_dimer_e-08
192.227217i
349.603371i
490.980356i
1252.606139i
1516.792244i
3179.233453i

TS3_climb_e-06
34.638204i
70.599026i
297.00951i

TS2_dimer_e-08 looks weird since it gives many frequencies in spite of a lower ediff value. All the important files related to these calculations are attached below.

I am sorry if the number of transition states is many but it will definitely be of great help if you can take a look and provide suggestions to fix these transition states.

[graeme]

In all of these calculations, at least one of the ionic steps reaches the maximum electronic number of iterations (60). If this ever occurs in a vibrational calculation, you can't keep any of it because the force at that step is not converged. To fix this, you can increase NELM and/or change your electronic optimizer to make it more efficient.

By the way, there is never a downside to use EDIFF=1e-8 in frequency calculations. It may take a few more electronic iterations, but there is no reason not to have precise forces, especially when you have soft modes in your system, as you will with largish molecules dissociating on a surface. If you are struggling getting down to EDIFF=1e-8, you need to solve that problem before doing frequency calculations.

[spbatchu]

Thanks, Dr. Henkelman for the suggestion. I will change the optimizer and nelm and get back to you.

Just a follow-up question, in the files that I have sent, I see that the number of electronic iterations exceeds the maximum number only in the case of TS2_dimer_e-08, and in the other TS cases(TS1, TS3, and TS2_e-06), I did not see the electronic iterations exceeding the maximum(please correct me if I am wrong). What do you think could have led to the additional imaginary freq. in these cases where the number of electronic steps hasn't exceeded the maximum and the potim is 0.002?

[graeme]

Yes, that is right. I think there is also a problem with soft modes in the dissociation of a molecule from a surface. Check the frequencies in your product state - I would not be surprised if you see the same near-zero or negative modes in what should be a stable state. I hate to suggest this, but as well as always keeping ediff=1e-8 and making sure that you converge every ionic step, you might have to increase the displacement magnitude to account for the softest modes, maybe to 0.005 or 0.01 Ang. Ideally you could step along each mode with a distance that makes sense for the stiffness of that mode, but with Cartesian displacements, you have to consider both C-H bond frequencies and the near-zero motion of the dissociating molecule. You have a difficult task ahead - good luck. And remember, at some point, there are modes which will not be well-described by harmonic oscillators.

[spbatchu]

Thanks Dr. Henkelman. Yes, I do see imaginary frequencies in the final state in the case of TS1 and TS3, but not in TS2.

Just a follow-up question, by any chance is there a way (theory,/thumb-rule/heuristic) to estimate the stiffness of a mode and choose a corresponding/specific potim/displacement value (or) is trying out different potim/displacement values( such as 0.005,0.01 Ang) in a trial and error fashion the only way to go about this issue?

[graeme]

If you know the curvatures and the accuracy in the forces, then yes, you can determine the smallest displacement that is appropriate. Unfortunately, the displacements in vasp are along cartesian coordinates and so you can have a mix of both very soft and hard modes in a displacement. For example, the displacement of a H atom can be both part of a stiff C-H bond and a floppy mode of the molecule. So in practice, you basically want to look for convergence in the modes as you lower the displacement magnitude and before you get into the noise of the forces.

[spbatchu]

Hello Dr. Henkelman

I have redone the vib freq. calculation with a potim of 0.002 Ang and ediff of e-07 and I see that the additional imaginary frequencies are going away. Also, having a NELM value higher than 150 or 200 helps. Thanks a lot.

[graeme]

Very good to hear!