Questions about DynMat calculation of gas-phase methanol

[jianmin]

Dear guys,

I tried to get the vibrational frequeies of methanol (CH3-OH, which has 6 atoms). The following is what I did.

Step 1: Fully relax one methanol molecule in a cubic box of 15 Angstroms side length. The key INCAR parameters are
NBANDS = 10
NSW = 150
IBRION = 2
EIDFF = 1.0E-7

In the OUTCAR file, the forces on ions are
POSITION TOTAL-FORCE (eV/Angst)
-----------------------------------------------------------------------------------
8.31684 8.01528 7.40000 0.000463 -0.000678 0.000000
8.37032 6.47450 6.50199 -0.000289 0.000105 -0.000434
8.37032 6.47450 8.29801 -0.000289 0.000105 0.000434
6.20195 6.11195 7.40000 0.000752 -0.001965 0.000000
7.97016 6.97500 7.40000 -0.001047 0.000969 0.000000
6.54041 7.02279 7.40000 0.000410 0.001464 0.000000
-----------------------------------------------------------------------------------
which, I believe, is accurate enough for the dynamical matrix calculation.


Step 2: perform a dynamical matrix calculation using the relaxed structure from step 1. The DISPLACECAR file looks like
0.001 0.001 0.001
0.001 0.001 0.001
0.001 0.001 0.001
0.001 0.001 0.001
0.001 0.001 0.001
0.001 0.001 0.001

And the key INCAR parameters are
IBRION = 3
POTIM = 0.0
ICHAIN = 1
NSW = 19 ! 3*6+1=19
EDIFF = 1.0E-6

Step 3: I issued the command
dymmatrix.pl DISPLACECAR OUTCAR

and I got the freq.dat which looks like
642.147364 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} ... 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.000013 cm^{-1} ... 0
1640.302054 cm^{-1} ... 0


Could this be right? because according to Chem. Eng. Sci. 59(2004)4679, the vibrational frequencies (cm^-1)of methanol are
3859, 3086, 3007, 2954, 1452, 1442, 1420, 1315, 1129, 1051, 1005, 285

The results I got are not even close. What do I do to solve this problem?

[graeme]

Try setting ediffg to a very small value, such as ediffg=-1e-8. This will force vasp to run all 19 iterations. My guess is that it only ran one, decided that it was converged, and then quit.

We'll work on a fix for this bug in our code.

[jianmin]

Prof. Graeme, thanks a lot for your reply.

I changed the EDIFFG to -1.0E-8 and then got all the frequencies as following:
352.208906 cm^{-1} ... 1
306.487546 cm^{-1} ... 1
87.364247 cm^{-1} ... 1
48.395025 cm^{-1} ... 1
35.040875 cm^{-1} ... 1
31.075402 cm^{-1} ... 1
9.693041 cm^{-1} ... 1
5.927690 cm^{-1} ... 1
2.150085 cm^{-1} ... 1
298.682645 cm^{-1} ... 0
997.634423 cm^{-1} ... 0
1076.182743 cm^{-1} ... 0
1384.823649 cm^{-1} ... 0
1468.123860 cm^{-1} ... 0
1534.949793 cm^{-1} ... 0
2759.019295 cm^{-1} ... 0
2967.691636 cm^{-1} ... 0
3730.150633 cm^{-1} ... 0

Are these results correct? Half of them are imaginary and half of them are real. It should have 3*6-6=12 regular vibrational modes, right?

[graeme]

Something is still wrong. Try two different tests

(a) increase the displacement to 0.005
(b) reduce ediff to 1e-8, using the original displacement of 0.001

If you get more reasonable modes with the more accurate forces (test b) then also try using a gaussian or fermi smearing with a low value of sigma. Make sure that your modes are converged with respect to small values of sigma.

Finally, just make sure that vasp always runs all 19 steps and that there are no zeros in the resulting dynamical matrix.

[jianmin]

Hi Prof. Graeme, I have done the following tests according to your suggestions.

(a) Increase the displacement to 0.005, and I got
177.080571 cm^{-1} ... 1
129.632840 cm^{-1} ... 1
48.307346 cm^{-1} ... 1
38.469833 cm^{-1} ... 1
26.444106 cm^{-1} ... 1
18.776305 cm^{-1} ... 1
12.635260 cm^{-1} ... 1
8.085865 cm^{-1} ... 1
3.005838 cm^{-1} ... 1
156.233757 cm^{-1} ... 0
1006.154487 cm^{-1} ... 0
1076.736518 cm^{-1} ... 0
1341.193465 cm^{-1} ... 0
1477.099550 cm^{-1} ... 0
1522.743368 cm^{-1} ... 0
2771.905633 cm^{-1} ... 0
3016.084249 cm^{-1} ... 0
3723.009372 cm^{-1} ... 0
I still got some imaginary modes.

(b) Reduce EDIFF to 1e-8, using the original displacement of 0.001
75.902602 cm^{-1} ... 1
6.196534 cm^{-1} ... 1
5.045446 cm^{-1} ... 1
1.996285 cm^{-1} ... 1
1.524069 cm^{-1} ... 1
0.657294 cm^{-1} ... 1
0.304232 cm^{-1} ... 1
4.305527 cm^{-1} ... 0
8.011507 cm^{-1} ... 0
62.020225 cm^{-1} ... 0
1007.813556 cm^{-1} ... 0
1054.311553 cm^{-1} ... 0
1326.451739 cm^{-1} ... 0
1431.061811 cm^{-1} ... 0
1462.841596 cm^{-1} ... 0
2933.008496 cm^{-1} ... 0
3054.688854 cm^{-1} ... 0
3747.087017 cm^{-1} ... 0
Similar to (a).

What should I do next?

[graeme]

I should have been a little more explicit: when you change sigma, use different orders of magnitude such as 0.1, 0.01, 0.001. See if there is any occupancy in the excited state.

I'll also give this a try and see what I get.

[jianmin]

Hi Prof. Graeme, I have run 3 jobs from comparison.

a 2x2x2 k mesh for the cubic box of 15 Angstrom side length, which has 2 k points.
EDIFFG = -1.0E-8
EDIFF = 1.0E-8
ISMEAR = 1
DISPLACEMENT = 0.001

_______________________________________________
Job 1: When SIGMA = 0.1
k-point 1 : 0.2500 0.2500 0.2500
band No. band energies occupation
1 -25.3962 2.00000
2 -16.8292 2.00000
3 -12.4722 2.00000
4 -10.6138 2.00000
5 -10.3314 2.00000
6 -7.9992 2.00000
7 -6.3114 1.99999
8 -0.7439 0.00000
9 -0.1284 0.00000
10 0.2627 0.00000
k-point 2 : -0.2500 0.2500 0.2500
band No. band energies occupation
1 -25.3962 2.00000
2 -16.8292 2.00000
3 -12.4722 2.00000
4 -10.6138 2.00000
5 -10.3314 2.00000
6 -7.9991 2.00000
7 -6.3114 2.00001
8 -0.7439 0.00000
9 -0.1261 0.00000
10 0.2404 0.00000

_______________________________________________
Job 2: When SIGMA = 0.01
k-point 1 : 0.2500 0.2500 0.2500
band No. band energies occupation
1 -25.3962 2.00000
2 -16.8292 2.00000
3 -12.4722 2.00000
4 -10.6138 2.00000
5 -10.3314 2.00000
6 -7.9992 2.00000
7 -6.3114 2.00000
8 -0.7439 0.00000
9 -0.1284 0.00000
10 0.2627 0.00000
k-point 2 : -0.2500 0.2500 0.2500
band No. band energies occupation
1 -25.3962 2.00000
2 -16.8292 2.00000
3 -12.4722 2.00000
4 -10.6138 2.00000
5 -10.3314 2.00000
6 -7.9991 2.00000
7 -6.3114 2.00000
8 -0.7439 0.00000
9 -0.1261 0.00000
10 0.2404 0.00000

_______________________________________________
Job 3: When SIGMA = 0.001
k-point 1 : 0.2500 0.2500 0.2500
band No. band energies occupation
1 -25.3962 2.00000
2 -16.8292 2.00000
3 -12.4722 2.00000
4 -10.6138 2.00000
5 -10.3314 2.00000
6 -7.9992 2.00000
7 -6.3114 2.00000
8 -0.7439 0.00000
9 -0.1284 0.00000
10 0.2627 0.00000
k-point 2 : -0.2500 0.2500 0.2500
band No. band energies occupation
1 -25.3962 2.00000
2 -16.8292 2.00000
3 -12.4722 2.00000
4 -10.6138 2.00000
5 -10.3314 2.00000
6 -7.9991 2.00000
7 -6.3114 2.00000
8 -0.7439 0.00000
9 -0.1261 0.00000

Conclusion: No excited state is occupied.

[graeme]

Another idea: try ISYM=0 in the INCAR

[graeme]

I hope that ISYM is the problem. I just ran the calculation and got the following modes:

1112.247957 cm^{-1} ... 1
393.774277 cm^{-1} ... 1
122.754647 cm^{-1} ... 1
73.095249 cm^{-1} ... 1
62.813129 cm^{-1} ... 1
179.409466 cm^{-1} ... 0
346.386662 cm^{-1} ... 0
496.782180 cm^{-1} ... 0
973.566811 cm^{-1} ... 0
1145.479796 cm^{-1} ... 0
1163.260906 cm^{-1} ... 0
1431.296239 cm^{-1} ... 0
1468.691291 cm^{-1} ... 0
1502.954331 cm^{-1} ... 0
1518.855128 cm^{-1} ... 0
2961.348592 cm^{-1} ... 0
3000.314115 cm^{-1} ... 0
3090.389060 cm^{-1} ... 0

I did not check for convergence with respect to displacement and sigma, but they already look a lot better than what you have reported. The symmetry flag can really mess up these calculations because if symmetry is detected for a displaced configuration, forces along the symmetric directions are artificially zero.

Here is my INCAR file. The displacement was 0.001 for each atom in each direction.

SYSTEM = methanol
IBRION=3
POTIM=0.0
ISYM=0
NSW=16
EDIFF=1E-06
EDIFFG=-0.0001
PREC=normal
LWAVE=.FALSE.
LCHARG=.FALSE.
ICHAIN=1

[jianmin]

Hi Prof. Graeme, thank you so much for all of your help. I set ISYM to 0 and it worked. The following is results from my job:

64.784000 cm^{-1} ... 1
50.187061 cm^{-1} ... 1
25.267127 cm^{-1} ... 1
4.768032 cm^{-1} ... 1
36.103207 cm^{-1} ... 0
43.456248 cm^{-1} ... 0
281.443955 cm^{-1} ... 0
1009.954328 cm^{-1} ... 0
1053.461913 cm^{-1} ... 0
1135.336745 cm^{-1} ... 0
1328.090165 cm^{-1} ... 0
1429.789579 cm^{-1} ... 0
1450.750545 cm^{-1} ... 0
1464.200341 cm^{-1} ... 0
2933.758479 cm^{-1} ... 0
2979.969844 cm^{-1} ... 0
3055.826147 cm^{-1} ... 0
3747.286098 cm^{-1} ... 0

These results are very close to the results in the literature.

[Ammi]

graeme wrote:
> I hope that ISYM is the problem. I just ran the calculation and got the
> following modes:
>
> 1112.247957 cm^{-1} ... 1
> 393.774277 cm^{-1} ... 1
> 122.754647 cm^{-1} ... 1
> 73.095249 cm^{-1} ... 1
> 62.813129 cm^{-1} ... 1
> 179.409466 cm^{-1} ... 0
> 346.386662 cm^{-1} ... 0
> 496.782180 cm^{-1} ... 0
> 973.566811 cm^{-1} ... 0
> 1145.479796 cm^{-1} ... 0
> 1163.260906 cm^{-1} ... 0
> 1431.296239 cm^{-1} ... 0
> 1468.691291 cm^{-1} ... 0
> 1502.954331 cm^{-1} ... 0
> 1518.855128 cm^{-1} ... 0
> 2961.348592 cm^{-1} ... 0
> 3000.314115 cm^{-1} ... 0
> 3090.389060 cm^{-1} ... 0
>
> I did not check for convergence with respect to displacement and sigma, but
> they already look a lot better than what you have reported. The symmetry
> flag can really mess up these calculations because if symmetry is detected
> for a displaced configuration, forces along the symmetric directions are
> artificially zero.
>
> Here is my INCAR file. The displacement was 0.001 for each atom in each
> direction.
>
> SYSTEM = methanol
> IBRION=3
> POTIM=0.0
> ISYM=0
> NSW=16
> EDIFF=1E-06
> EDIFFG=-0.0001
> PREC=normal
> LWAVE=.FALSE.
> LCHARG=.FALSE.
> ICHAIN=1

Hi,

It is very old thread. However, I found it when browsing the forum.

About the effect of sigma value (Gaussian smearing) on the calculated frequencies for my crystal structure:

Could we check directly the effect of different sigma values (0.05, 0.01, 0.001) on the calculated frequencies for the original optimized structure with for example sigma=0.1?

or should first optimize the structure with the relevant sigma value then do frequency calculation on that structure?

Thanks

[graeme]

I would suggest checking both convergence of the structure and the frequencies as a function of sigma. So what I mean by that is to reduce sigma, re-optimize and then recalculate frequencies and check for convergence. But this is not typically a difficult parameter to converge - as long as you see on the order of a meV/atom of electronic entropy, it should be fine.

[Ammi]

Thank you so much.

Do you mean that: continue with contcar (of old sigma) and the new sigma to re-optimize the structure and calculate the frequencies until convergence? With this I see only one ionic step and no further ionic relaxation for my structure.

With entropy less than about 1meV/atom for the equilibrium (minimum) structure (with the preliminary sigma value), I got three negative frequencies of about 0.01, 0.02, 0.80 THz, with vasp phonon calculation and I allow all atoms to be free. What about these frequencies, especially the last one (0.80 THz or about 26 cm-1)?

Thanks

[graeme]

Well, I have no information about the system and so I can't say much. Typically, if you relax all of the atoms then you should have positive frequencies, however, if you have all atoms free and you have a nice set of positive frequencies and those three near-zero frequencies, they could correspond to translation, which is not a problem. If you attach the calculation files, I could probably be more helpful.

[Ammi]

Thank you so much.

I checked carefully the frequencies that I got previously which were using Ibrion6 with isif3. In addition to the last three negative frequencies, there are some small positive frequencies.
I have also tried ibrion8 with isif2 and addgrid. There are also three negative frequencies and these small positive frequencies increased.

Could both results be reasonable without problems? I try to check if the structure is a minimum and a fairly optimized. But those frequencies at the end of each file, were confused me.

How can we decide that the negative frequencies are small and can be ignored? could those negative frequencies and especially that of about 0.80 THz be ignored when have these small positive frequencies?

Also, for my system the entropy is low, with reducing the smearing size, it seems no more change. I am checking it more.

Thanks