Frequency calculations

[baki]

Hi everybody
Recently I try to calculated the frequency of normal CH4
I used the following InPUT file
general:
ISTART = 1
ICHARG = 0
SYSTEM = clean (111) surface
ENMAX = 400
ALGO=V
NPLANE=4
NPAR=4
LPLANE=.TRUE.
EDIFFG=-0.025
EDIFF=0.000001
PREC=ACCURATE
ISYM=0
dynamic:
NSW=1
IBRION=5
POTIM=0.02
NFREE=2
~
I got several very big negative frequency , but I don't get the mode around 1500 and 1300 cm-1.
My box size was about 12 A.
best

[admin]

-- all vibration frequencies should be positive.
-- PREC should be set to High
-- check whether the box size was large enough
-- check whether your converged equilibrium geometry is reasonable
-- mind that the used type of the Pseudopotentials (LDA/GGA) has significant
influence on the bond strengths and hence on the vibrational spectra.

[baki]

Hi
What is the main difference between the PREC=HIGn and Prec=ACCurate, if I don't use the LREAL=.TRUE. or LREAL=.AUTO. option?
1. I used PAW/PBE pseudopotentials.
2. My geometry is OK. Input was some equilibrated CH4.
3. I try to use the ISYM=0 or ISYM=2 -no effect
My POSCAR file was the following
12.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 12.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 12.0000000000000000
4 1
Selective dynamics
Direct
0.1747173311070082 0.0833346212016710 0.0833340940928887 T T T
0.0528745470065775 0.1694916606418204 0.0833333889773981 T T T
0.0528720796287767 0.0402545266053759 0.0087195881664012 T T T
0.0528725788596435 0.0402536844633011 0.1579469233004847 T T T
0.0833309633979967 0.0833321737545028 0.0833326721294988 T T T
best

[admin]

the PREC tag determines the default settings of the variables
ENCUT,
NGX(F), NGY(F), NGZ(F),
ROPT
please have a look at the VASP online manual for further details
http://cms.mpi.univie.ac.at/vasp/vasp/node102.html

concerning the vibrational frequencies:
The errors you see are most likely DFT related
and / or related to the neglect of zero point vibrations (quantum effect):
i.e. you treat the hydrogen as a classical atom, whereas it is so light,
that one should treat it as a quantum particle. When you treat it
as quantum particle the average bond length increases somewhat, since
the potential energy surface is not exactly quadratic (anharmonic effect).
VASP does not allow to treat hydrogen as a quantum particle (actually
there practically no code that can do this...

[graeme]

I have also struggled with this problem, and have a few things to add. The first is that increasing the accuracy of the calculations do not seem to help. Going to a very high cutoff, large box, fine fft grid, very accurate forces (ediff=1e-8), and small finite difference displacement steps all give similarly converged frequencies, with positive modes (there can be negative modes as well) around:

4251, 3554, 2961, 2945, 2944, 2286, 1754, 894, 890

where the experimental values are

3019, 2917, 1534, 1306

One interesting point is that the vasp modes are fairly sensitive to the orientation of gas phase molecules in the cell. This is particularly apparent with a linear molecule such as H2, in which negative modes show up if the molecule is not aligned along the Cartesian axes, even for a very large box of 20 Ang in each direction.

Another point is that a 10 second calculation with spartan gives fairly nice modes:

3114, 3110, 3108, 2966, 1512, 1512, 1309, 1292, 1287

so this rules out the quantum effects for the hydrogen nucleus.

Now, a very interesting thing happens if you look at CH4 above a metal surface. The following modes were calculated with medium precision, default energy cutoff, small box, and regular fft grid:

3102, 3076, 3073, 2963, 1501, 1499, 1399, 1298, 1297

They are vastly improved! There is something interesting going on here, and it would be great to have a vasp guru help figure it out.

BTW, these frequencies were calculated using the finite difference dynamical matrix module of Blas Uberuaga: http://theory.cm.utexas.edu/vtsttools/ . However, I think this is a vasp issue, and not related to the details of how the modes are calculated.
<span class='smallblacktext'>[ Edited Sat Jul 02 2005, 08:50AM ]</span>

[scerwin]

Update (27 July): since the tables in this post apparently can't be formatted properly, I suggest you instead read the more legible version at
http://cst-www.nrl.navy.mil/~erwin/vasp ... ibrational.

-scerwin

------------------------------------------------------------------------------------------------------

Note: I can't seem to format the tables in this message properly. They look fine under "Preview" but not as an actual post! Forum designers, can you help? -scerwin

------------------------------------------------------------------------------------------------------

I am new to vibrational calculations, although not to vasp, and have
a few comments to add.

First, a few words about the "negative" (imaginary, actually)
frequencies that people are reporting. These are expected, or at
least are not surprising. For a molecule in a box, there should in
principle be six zero frequencies, corresponding to three rigid
translational modes and three rigid rotational modes. (For a
collinear molecule like H2, there are only two rigid rotational
modes.)

In numerical calculations, "zero" will be manifested as a small number
-- possibly real, possibly imaginary -- whose magnitude is expected to
scale with POTIM. This is exactly what happens in vasp. Here are the
rigid translational (t) and rotational (r) frequencies, for different
POTIM, calculated for the CH4 molecule:
<pre>
POTIM imaginary frequencies real frequencies
----- --------------------- ----------------
0.20 7.1 t(3), 532.3 r(3)
0.10 1.5 t(3), 266.6 r(3)
0.05 0.4 t(3), 123.1 r(3)
0.02 0.5 t(2), 6.0 t(1), 12.0 r(2) 6.7 r(1)
</pre>
I'll get to the vasp parameters I used later. For now, just notice
that the three rotational-mode frequencies scale with POTIM, as
expected. This scaling (and the degeneracies) begins to fall apart
for POTIM=0.02 -- a warning that this value might be too small for
calculating reliable numerical derivatives with PREC=ACCURATE
(everything recovers nicely if you manually increase the cutoff
energy, though).

Second, for the purpose of evaluating vasp results for the physically
meaningful modes, we should compare to the converged DFT all-electron
(AE) frequencies for CH4. These are given in a 1997 paper by Mark
Pederson (http://link.aps.org/abstract/PRB/v55/p7454, see Table III):
<pre>
irrep Exp AE/PBE (deg)
----- --- ------------
T2 3158 3081 (3)
A1 3037 2968 (1)
E 1567 1510 (2)
T2 1357 1268 (3)
</pre>
For my vasp calculations I am using PAW/PBE "standard" potentials ("C"
and "H"), a 6-Angstrom box, k=Gamma, PREC=ACCURATE, EDIFF=1E-6,
NFREE=2, ISMEAR=-1, SIGMA=0.1, and various values for POTIM. Here are
the results compared to Mark's:
<pre>
PAW/PBE
----------------------------
-- different POTIM values --
irrep Exp AE/PBE (deg) 0.10 0.05 0.02
----- --- ------------ -------- -------- --------
T2 3158 3081 (3) 3093 (3) 3104 (3) 3106 (3)
A1 3037 2968 (1) 2983 (1) 2993 (1) 2994 (1)
E 1567 1510 (2) 1495 (2) 1514 (2) 1519 (2)
T2 1357 1268 (3) 1279 (3) 1279 (3) 1287 (3)
</pre>
As you can see, the PAW/PBE frequencies are within about 1% of the
AE/PBE values. Interestingly, the results for POTIM=0.02 seem fine,
so maybe my warning above wasn't necessary.

So, how can we understand the several reports on this discussion
thread of CH4 frequencies quite different from these? The answer is
the occupation number scheme. Using the default Methfessel-Paxton
scheme (ISMEAR=1,SIGMA=0.2; which it appears "baki" is using) one gets
the following horrible results:
<pre>
ISMEAR=1,SIGMA=0.2
irrep Exp AE/PBE (deg) ------------------
----- --- ------------ --- POTIM=0.1 ----
T2 3158 3081 (3) 3121 (3)
A1 3037 2968 (1) 2979 (1)
E 1567 1510 (2) 932 (2)
T2 1357 1268 (3) 1337 (3)
</pre>
This is very strange behavior, because the HOMO-LUMO gap for CH4 is
over 8 eV in vasp, and so any scheme should give the same occupation
numbers (2 or 0). But Methfessel-Paxton has a problem. For the
unperturbed geometry, it indeed correctly gives exactly 2 or exactly
0. Same for perturbations of the carbon atom. But for the hydrogen
perturbations, the occupation numbers are as follows:
<pre>
band No. band energies occupation
1 -16.4853 2.00000
2 -9.0946 2.00015
3 -8.6978 2.06013
4 -8.6324 1.93972
5 -0.7906 0.00000
6 2.6483 0.00000
7 2.6622 0.00000
8 3.2954 0.00000
</pre>
This set of occupation numbers artifically decreases the total energy
below its correct value (by having occupation > 2 for the lower-energy
state #3). The consequence of this is to severely soften the twofold
"E" vibrational mode, as shown in the earlier table above.

I suspect this is also the reason that having a metal surface in the
calculation somehow "fixes" the bad frequencies. But the underlying
problem is with the occupation numbers.

Bottom line: the vasp vibrational frequencies for CH4 are fine, but
don't use Methfessel-Paxton occupation numbers.



<span class='smallblacktext'>[ Edited Wed Jul 27 2005, 05:01PM ]</span>

[baki]

Hi scervin
Than you your quit good advice. I try to use the Fermi smearing and it gives very good results for alkane,(metan to hexane). Applying this method give frequencies which are very close to experimental or AE one. I have only a small problem with this. I want to use these data as a reference calculation comparing the adsorbed and free state of alkane. The vasp manual suggest to use for adsorption purpose the MP schema. Now I try to repeat my calculation wit fermi schema and investigated the difference. (I also will recalculated the normal mode of methanol, formaldehyde, furane and formic acis where I don't have any problem with the MP schema
Best
baki

[scerwin]

Dear baki,

I would not recommend using MP occupation numbers for any vibrational problem, whether a free molecule or an adsorbed one. The important thing is that the occupation numbers should not change AT ALL when the small displacements are made to calculate the dynamical matrix. If they do, your force constants will be wrong and so will your frequencies.

-scerwin