Frequencies in a dynamical matrix calculation

[paspanin]

Hello everybody!

I have searched the forum for replies to my question but none of the proposed suggestions seems to work. I am simulating a vacancy-mediated diffusion in a 4x4x4 bcc supercell. My aim is to find the prefactor after obtaining the dynamical matrix and finding the frequency modes using your VTST tools. First of all, I have used the NEB method to obtain the migration barrier and transition state. The initial equilibrium state and the transition state have been converged using EDIFF= 1E-8.

Firstly, I created a DISPLACECAR using the script the script dymseldsph.pl. I started by displacing only one of the nearest-neighbour atoms to the vacancy by 0.001 Ang, using 1.0 Ang as a radius from that atom. This would be the atom that moves to fill the vacancy during the migration process. This means only three DOF for each of the equilibrium and transition states. Thus, NSW=3+1=4 in INCAR. I set ISYM=0 as well. After using the dymmatrix.pl script, I obtain three real frequencies for the equilibrium initial state and one imaginary and two real frequencies for the transition state, as expected. The dymprefactor.pl script allows me to obtain a value for the prefactor.

Since the comparison to experimental values does not seem quite right, I was wondering if more atoms should be included to create the DISPLACECAR file, since it is certain that more than one atom moves. I thought of expanding the radius from the central atom to be used with the dymseldsph.pl script to 3.0 Ang and thus eight atoms are displaced by 0.001 Ang, meaning NSW=24+1=25. However, in my first attempt to obtain the frequency modes for the equilibrium state, all frequencies are imaginary and identically zero.

I have also tried to create the DISPLACECAR file using the dymseldsp.pl script and the equilibrium and transition state POSCARs as POSCAR1 and POSCAR2 files, to no avail.

I would like to ask if I am forgetting an obvious step that should be taken when more than one atom is displaced. Maybe it is enough to displace just one of the closest atoms in this kind of vacancy-mediated processes and I should be happy with the value found for the prefactor in the one-atom calculation. And so find better experiments to compare to!

Of course, I'd be glad to post my INCAR and DISPLACECAR files if required.

Thanks a lot in advance!

Cheers
Roberto

[graeme]

Hi Pavel,

You should always increase the size of your dynamical matrix to check for convergence. As you add more atoms away from the diffusing atom, the extra modes at the minimum and saddle will be similar and so cancel in the prefactor.

Everything in your description sounded fine until you say that the 24 DOF frequencies are all imaginary or zero. This does not make any sense. Check for some simple problems. First, make sure that the calculation finished. I'm not sure if our script warns you if there is missing data or not. It may generate a matrix with zeros for modes that are not calculated. Second, check the Hessian matrix itself and see if there are any zeros, missing data, or crazy values.

I expect there is some technical problem, rather than anything to do with the physics.

Also, for solid state diffusion, your prefactor should be close to 10^12 - 10^13 /s; it would be surprising get values much outside this range. Therefor, the frequency calculation is unlikely to change any agreement or disagreement with experiment. Nevertheless, it would be good to have a converged prefactor.

One more point: if you have no frozen atoms, remember that there will be 3 zero modes associated with translation. I think that our prefactor script removes these for you, but if not, they should not be included in the prefactor calculation. Technically they cancel, but noise can make them small and positive or imaginary, which will mess up the prefactor ratio.

[paspanin]

Dear Graeme,

Thanks a lot for your reply! You were absolutely right. I made a stupid mistake in the calculation and used the wrong executable, one not compiled using your VTSTTools package. The prefactor obtained for the single atom displacement was also of the order of 10^13 Hz as you mentioned.

Now I've redone both equilibrium and saddle point simulations with the DISPLACECAR created with eight atoms around a vacancy. It seems that I still have a problem, because I find two imaginary modes at the equilibrium point and more than three imaginary modes at the saddle point, respectively:

Equilibrium point:
25.492611 cm^{-1} ... 1
14.183293 cm^{-1} ... 1
0.000003 cm^{-1} ... 1
0.000001 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000001 cm^{-1} ... 0
0.000002 cm^{-1} ... 0
112.677596 cm^{-1} ... 0
132.453746 cm^{-1} ... 0

Saddle point:

36.268326 cm^{-1} ... 1
25.871894 cm^{-1} ... 1
6.038746 cm^{-1} ... 1
3.720160 cm^{-1} ... 1
3.548007 cm^{-1} ... 1
3.376108 cm^{-1} ... 1
11.675914 cm^{-1} ... 0
45.710378 cm^{-1} ... 0
103.910225 cm^{-1} ... 0
112.213908 cm^{-1} ... 0
112.764183 cm^{-1} ... 0
134.185418 cm^{-1} ... 0
148.284670 cm^{-1} ... 0
151.293108 cm^{-1} ... 0
153.182875 cm^{-1} ... 0
154.019781 cm^{-1} ... 0
154.485880 cm^{-1} ... 0
156.140260 cm^{-1} ... 0
156.657454 cm^{-1} ... 0
159.282982 cm^{-1} ... 0
161.684293 cm^{-1} ... 0
163.330305 cm^{-1} ... 0
168.740199 cm^{-1} ... 0
178.074114 cm^{-1} ... 0

With these two imaginary frequencies in the initial state, the message "Initial state freq.dat has imaginary frequencies." appears

I would be grateful for any suggestions. And please don't hesitate to ask if some other particular information is needed.

Thanks in advance!

Roberto

[graeme]

Hi Roberto,

Things are still a little crazy. Those zero frequencies indicate that the finite difference calculation of the minimum did not finish or had some technical problem. If you look at the Hessian matrix, you should see a number of displacements for which no force was calculated.

The saddle may be ok, but only if you have structure surrounded by vacuum to give 3 rotation and 3 translation zero modes. If it is a bulk or surface structure with constrained rotational or translational modes, there should not be so many near-zero modes.

[paspanin]

Hi Graeme!

Thanks a lot again for your reply! I've checked the Hessian matrix as shown in freq.mat and it's full of zeros! So there are many displacements for which no frequency was calculated. However, I don’t know where to look to check if there was a problem with the calculation at the minimum, as you suggest. Should I redo the previous relaxation at step 1 to find if there are errors or warnings of any kind? I guess the problem comes from the structure at the minimum, and not from the actual TST Tools DynMat calculation, doesn’t it?
Again, the saddle should be wrong as well, since it is not a structure in a vacuum, but a vacancy created in a bulk system.

Thank you very much for your suggestions!

Roberto

[graeme]

Check to see if vasp completed the total number of force calls corresponding to the number of displacements that you specified. If not, see if it quit due to a time limit, or something else. If it claimed to converge, set EDIFF very low (e.g. 1e-10) so that the calculation will not terminate prematurely.

[paspanin]

Hi again and thanks for your answer, Graeme!
As you suggested, I re-checked the number of force calls in my previous try, and instead of 25 force calls, corresponding to NSW=25 (8 atoms x 3DOF)+1, I found only 3. There is no walltime limit set in the machine, a local cluster, and I found no signs of premature termination of the job, the messages "reached required accuracy - stopping structural energy minimisation - writing wavefunctions" are there at the end of the OUTCAR file. Thus, I decided to try setting EDIFF=1E-10. Please find below the diff between the freq.mat files before/after:

< 25.255906 cm^{-1} ... 1
< 13.779987 cm^{-1} ... 1
---
> 25.492611 cm^{-1} ... 1
> 14.183293 cm^{-1} ... 1
4c4
< 0.000000 cm^{-1} ... 1
---
> 0.000001 cm^{-1} ... 1
22,24c22,24
< 0.000001 cm^{-1} ... 0
< 112.835029 cm^{-1} ... 0
< 131.367465 cm^{-1} ... 0
---
> 0.000002 cm^{-1} ... 0
> 112.677596 cm^{-1} ... 0
> 132.453746 cm^{-1} ... 0

Indeed almost none! The two imaginary frequencies at the equilibrium point are still there. The number of force calls is still 3.

Remember that I created the DISPLACECAR using the dymseldsp.pl script and as POSCARs the equilibrium and saddle point arrangements, respectively? I also tried setting as the central atom one of the 8 neighbouring atoms to the vacancy in this bcc cell and using the dymseldsph.pl script. Maybe this assymetry causes some confusion that ends up the job earlier than expected. Should I try to place a fictitious atom in the center of the vacancy and rebuild my DISPLACECAR?

Thank you very much for your patience and most valuable input!

Greets!
Roberto

[graeme]

Ah, my bad, I mean to set EIDFFG=1e-10 (not ediff, which should be about 1e-8).
Sorry about that (and thanks for your patience)!

[paspanin]

Dear Graeme,

Thank you very much for your excellent piece of advice!! Setting EDIFFG=-10 made a change and now there are 24 force calls, exactly the number of displacements in my DISPLACECAR file. As a result, all frequencies are real at the saddle point; however, the first one at the equilibrium point is imaginary. Actually, I would expect the opposite behaviour. Is this normal, related to a translation, for instance? Although I guess there should be three imaginary frequencies, and not one, in that case. Should I thus simply delete it before calculating the prefactor? (The prefactor.pl scripts complains, as "Initial state freq.dat has imaginary frequencies."). Or is there still something wrong?

Thank you very much again!

Cheers

Roberto

[graeme]

There should be one negative mode at the saddle and all positive modes at the minimum. Yes, you can remove 3 near-zero modes in both calculations.