UT theoretical chemistry code forum

Hello everyone,
I am doing a simple vacancy-mediated diffusion calculation for a bulk system. I would like to calculate the prefactor ab-initio by generating the dynamical matrix and calculating the frequency modes. I have calculated the transition state using the NEB method. My question is.
If I use the script dymseldsp.pl to calculate the atoms to be included I use the initial state POSCAR as the first one and the saddle point POSCAR as the second POSCAR. So while calculating the DYNAMICAL matrix, the DISPLACECAR file will be the same for both the runs.. Is this right?
Also, how to determine the convergence ? Is there a set criteria?

Regards,
PG

Hello everyone,
I managed to set up the dynamical matrix calculation using the DISPLACECAR file obtained by the dymseldsph.pl script and I included 40 atoms. I ran two calculations one for the initial state and one for the saddle point using a displacement of 0.001. However for the initial state I see quite a few imaginaory frequencies. The output of freq.dat is as follows.
225.555301 cm^{-1} ... 1
167.349604 cm^{-1} ... 1
147.995844 cm^{-1} ... 1
138.353004 cm^{-1} ... 1
125.659931 cm^{-1} ... 1
64.713378 cm^{-1} ... 1
43.670636 cm^{-1} ... 1
34.330809 cm^{-1} ... 1
30.558527 cm^{-1} ... 1
24.621742 cm^{-1} ... 1
21.344564 cm^{-1} ... 1
17.955455 cm^{-1} ... 1
16.946047 cm^{-1} ... 1
15.748004 cm^{-1} ... 1
14.037267 cm^{-1} ... 1
12.352689 cm^{-1} ... 1
11.227785 cm^{-1} ... 1
9.277767 cm^{-1} ... 1
7.816974 cm^{-1} ... 1
0.000003 cm^{-1} ... 0
0.000003 cm^{-1} ... 0
0.000003 cm^{-1} ... 0
0.000003 cm^{-1} ... 0
0.000003 cm^{-1} ... 0
0.000004 cm^{-1} ... 0
0.000004 cm^{-1} ... 0
1.192523 cm^{-1} ... 0
64.137814 cm^{-1} ... 0
77.793091 cm^{-1} ... 0
88.960582 cm^{-1} ... 0
104.603070 cm^{-1} ... 0
111.006227 cm^{-1} ... 0
112.237804 cm^{-1} ... 0
126.094405 cm^{-1} ... 0
130.791377 cm^{-1} ... 0
139.511754 cm^{-1} ... 0
143.663884 cm^{-1} ... 0
155.172254 cm^{-1} ... 0
159.671237 cm^{-1} ... 0
177.198437 cm^{-1} ... 0
188.103646 cm^{-1} ... 0
202.655438 cm^{-1} ... 0
224.448462 cm^{-1} ... 0
231.515938 cm^{-1} ... 0

I am not sure what this means? Is there a problem with my convergence? The NEB calcualtions were well-converged with the residual forces on the saddle point less than 0.005 eV/Ang.

-Regards,
PG

i didnt try the dynamical matrix calculation; but i tried to use ibrion=5 as a quick test;

the extremely small negative frequencies are generally representing lateral displacment, or translation which are not real; only the one roughly along rxn coordinates, strech of an old bond is real

Dear Kai,
Thanks for your reply. Could you please elaborate? I have never done these calculations before. How do we determine which of the frequencies are real? There seem to be a lot of imaginaory frequencies. Also, I am interested in calculating the prefactor for which I obtained the dynamical matrix of the initial and the saddle points and then tried using dymprefactor.pl. But i get an error message " the initial matrix contains imaginary frequencies". I am using the VTST tools for all my calculations. Do you think it is better to use the IBRION= 5 within VASP for the calculations?

Thanks,
PG

Hello,
I am trying to find the problem with my dynamical matrix calculation. In the previous calculations I did not set ISYM = 0 and I think that was the problem for getting so many imaginary frequencies in the initial state. I repeated the calculations with ISYM = 0 and calculated the Dynamical matrix. I tried for two different displacements (0.001) and (0.005) keeping the same number of atoms (60). For displacement by 0.005 Ang, the freq.dat looks fine as below.
92.611177 cm^{-1} ... 0
95.798917 cm^{-1} ... 0
102.822516 cm^{-1} ... 0
111.750235 cm^{-1} ... 0
112.737276 cm^{-1} ... 0
114.813537 cm^{-1} ... 0
116.406092 cm^{-1} ... 0
117.806745 cm^{-1} ... 0
121.055968 cm^{-1} ... 0
122.040592 cm^{-1} ... 0
124.561748 cm^{-1} ... 0
128.657307 cm^{-1} ... 0
129.346438 cm^{-1} ... 0
130.510711 cm^{-1} ... 0
133.237403 cm^{-1} ... 0
134.093800 cm^{-1} ... 0
134.988669 cm^{-1} ... 0
136.311983 cm^{-1} ... 0
139.024017 cm^{-1} ... 0
139.754638 cm^{-1} ... 0
141.558529 cm^{-1} ... 0
142.534256 cm^{-1} ... 0
142.744540 cm^{-1} ... 0

However, changing the displacement to 0.001 gives me a completely different set with one imaginary frequency.
52.497053 cm^{-1} ... 1
58.405475 cm^{-1} ... 0
61.294414 cm^{-1} ... 0
69.356062 cm^{-1} ... 0
75.085752 cm^{-1} ... 0
83.061264 cm^{-1} ... 0
86.016231 cm^{-1} ... 0
89.176865 cm^{-1} ... 0
90.362848 cm^{-1} ... 0
94.483989 cm^{-1} ... 0
98.410340 cm^{-1} ... 0
99.654795 cm^{-1} ... 0
103.837027 cm^{-1} ... 0
108.079008 cm^{-1} ... 0
108.774139 cm^{-1} ... 0
108.830177 cm^{-1} ... 0
111.771196 cm^{-1} ... 0
113.860633 cm^{-1} ... 0
116.875441 cm^{-1} ... 0
117.970821 cm^{-1} ... 0
120.934304 cm^{-1} ... 0
121.132538 cm^{-1} ... 0
122.785509 cm^{-1} ... 0

I am not sure what the problem is. How does one determine the convergence in such calculations? These calculations take too much time and so I would like to know a better way to test convergence.. Could someone help me with this?

Regards,
PG

As you go to smaller displacements, your value of ediff will also have to be small. For a very small displacement of 0.001, you will need ediff=1e-8.

You should find converged frequencies between 0.05 and 0.005 Ang, and also between 0.005 and 0.001 if your forces are accurate enough.

Dear Graeme,
Thanks for the reply. I will try to restart my calculation with the new settings. The saddle point in my calculation are converged well until the force on each atom <0.005 eV/Ang. I guess this should be a good force convergence.

Regards,
PG

Regarding computing the vibrational frequencies at the saddle point: Let's say that we are done with the NEB part and we obtained a saddle point configuration. Does one need first to minimize the energy of the saddle point configuration under periodic Boundary conditions before the setup of the dynamical matrix calculations? If yes, why would be that necessary?
Thanks in advance

Assuming that you are calculating the frequencies at the saddle point within vasp, there is no need to apply boundary conditions.

That said, the output configuration from an NEB calculation should be consistent with the boundary conditions.

So as in the original post of this thread, I'm modeling oxygen vacancy diffusion in bulk oxide and tried to compute the pre-exponential and the energy barrier. The NEB part converged with forces less than 0.01 eV/A.
However, when I computed the vibrational frequencies at the initial state, I obtained three "small" imaginary frequencies and the rest are real. At the saddle point I obtained 3 small imaginary frequencies + one large imaginary frequencies and the rest are real. So for the initial state I got the following imaginary modes :
283 f/i= 0.019799 THz 0.124403 2PiTHz 0.660434 cm-1 0.081883 meV
284 f/i= 0.019799 THz 0.124403 2PiTHz 0.660434 cm-1 0.081883 meV
285 f/i= 0.019799 THz 0.124403 2PiTHz 0.660434 cm-1 0.081883 meV

And for the saddle point , I obtained:
282 f/i= 0.013194 THz 0.082899 2PiTHz 0.440100 cm-1 0.054565 meV
283 f/i= 0.013235 THz 0.083161 2PiTHz 0.441487 cm-1 0.054737 meV
284 f/i= 0.247360 THz 1.554207 2PiTHz 8.251034 cm-1 1.022998 meV
285 f/i= 8.323981 THz 52.301112 2PiTHz 277.658096 cm-1 34.425229 meV


I checked the convergence with respect to POTIM (0.001 up to 0.01), I used EDIFF =1e-8 and 1e-7 , I used PREC=High and accurate, I tried ISYM=0 and ISYM=2. Regardless of what I do, I always obtain these three small imaginary frequencies. Notice also, that I displace all the atoms in a 3D periodic solid. My first Guess, this is just a numerical error and those three imaginary frequencies can be neglected. Any input regarding that will be very appreciated.

Another point, in VASP forums somebody recommended performing geometry optimization for the Saddle point configuration (AFTER NEB convergence) before actually trying to compute the saddle point vibrational frequencies. Is that a recommended practice and why?

Thank you very much!

If you have no frozen atoms, you will find three zero modes corresponding to translation of the system. These should be removed from both the initial and saddle in your prefactor calculation. Freezing a single atom in your calculation will also serve to remove these modes.

If the climbing image is well converged, there is no need to do anything before a dynamical matrix calculation. If, however, you are running an expensive NEB with many images and you are impatient to calculate an accurate saddle and prefactor, then you can converge just the saddle point in a separate calculation. There are a few ways to do this. You can run a single-image NEB using the climbing image as your moving image. You can use the neb2dim script and run a dimer calculation. FInally, you can try using the quasi-newton or (l)bfgs optimizers to converge to the saddle. This last option carries the risk that you will move away from a saddle if you are not close enough.

Regardless of how you do it, you want the force on the saddle geometry to be low (<0.01 eV/Ang) before calculating vibrational frequencies.

Thank you very much for your answer and explanation! That helps a lot!