My Community

Dear all,

I have a question related to a rather simple molecule, H2, for which I have tried to calculate its vibrational mode following geometry optimisation. The result was surprising, but I must also admit that I am relatively new to VASP. It turns out that there are two very large imaginary modes present, which are energetically nearly degenerate, plus the larges vibrational mode is around 9000 cm-1. Compared to experiment (~4000 cm-1), this is way off.

1 f = 8976.739914 cm-1
2 f = 0.012922 cm-1
3 f = 0.000123 cm-1
4 f/i= 0.004516 cm-1
5 f/i= 4711.668966 cm-1
6 f/i= 4711.706323 cm-1

OK, that doesn't look good. Now, I understand that there are a number of factors that would influence this result, and after skimming through the forum I found posts reporting similar issues. However, in these cases the imaginary frequencies were smaller. My result suggests that something is going pretty wrong. I understand that for linear species the molecule should be aligned with the axis, a requirement which is fulfilled in my case. The cutoff value should be ok too, I am using the same value for other, larger, organometallic systems at the moment. In fact I need the H2 molecule as a reference for calculating free energies, since H2 is involved in the reaction of these organometallic species, and hence I need an absolutely reliable result here.
Box size should be fine (please note that in the example below I have used a somewhat smaller cubic box for the phonon calculation [8 Ang], whilst the optimisation was done for 20 Ang.; there is not effect in the result of the phonon calculation when a larger box size is used. This was just a test). SO, presently I would be reluctant to change anything about the energy cutoff, since this would mean I'd have to repeat the other calculations I guess.

Is this an issue about the geometry convergence. I mean, it can't be so hard to optimise H2, and the obtain bond distance of 0.749 Ang is perfectly acceptable from a structural point of view. I don't know what I should do here to get a tighter geometry.

Any comments would be much appreciated. Which parameters do I need to play around with, and is there a way of getting the "right" results in the first shot? To my eye the parameters I have used don't look too bad after all, but the result is of course unacceptable. Are there some rules of thumb to set up these type of jobs?

Thanks in advance.


Tobias


Please find my input parameters for both the geometry relaxation and phonon calculations below.




INCAR (phonons)

SYSTEM = H2
ISMEAR = 0 ! Gaussian smearing
IBRION = 5 ! use the conjugate gradient algorithm
NFREE = 2 ! central differences
ENCUT 600
POTIM = 0.02 ! 0.02 A stepwidth
NSW = 1 ! ionic steps > 0
Van der Waals Interaction (vasp 5.3.3 patched verion):
IVDW = 12 ! switches between 0:off, 1: DFT-D3 and 2: TS-VDW (default=1)

POSCAR (phonons)

H
1.00000000000000
8.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 8.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 8.0000000000000000
H
2
Selective dynamics
Cartesian
0.0000000000000000 0.0000000000000000 0.0000000000000000 T T T
0.0000000000000000 0.0000000000000000 0.7498000000000000 T T T


INCAR (relaxation)

General:
SYSTEM = H2
! if enough K-points are present -5 should be used
ISMEAR = 0 ! 0: Gaussian -> semicond/insulator ; 1-N: MP; -5: Tetra+Blochl -> metal/very accurate energy/forces
SIGMA = 0.05
EDIFF = 1.0E-7
PREC = Normal ! to make ROPT=2.0E-4 for LREAL=Auto
LREAL = Auto ! cheaper for large cell's we are using
ENCUT = 600 !
ISYM = 2 # use symmetry as done for PAW PP

NELMIN = 8
LWAVE = .FALSE.
LCHARG = .FALSE.
! LVTOT = .TRUE. ! Write LOCPOT file for Toon
! LVHAR = .TRUE. ! but LOCPOT without exchange and correlation contribution for forcefields Toon
VOSKOWN = 1 ! important for GGA (PW91) for interpolation of XC...since we only use PBE=switch on (LDA: VOSKOWN=0)
LASPH = .TRUE. ! For VASP.5.X the aspherical contributions are properly accounted for in the Kohn-Sham potential
! as well. This is essential for accurate total energies and band structure calculations for f-elements
! (e.g. ceria), all 3d-elements (transition metal oxides), and magnetic atoms in the 2nd row (B-F atom),
! in particular if LDA+U or hybrid functionals or meta-GGAs are used, since these functionals often result
! in aspherical charge densities.

Van der Waals Interaction (vasp 5.3.3 patched verion):
IVDW = 12 ! switches between 0:off, 1: DFT-D3 and 2: TS-VDW (default=1) 11: zero-damping 12: Becke-Johnson

dynamic:
IBRION = 2 ! -1: Fix atoms; 0: MD; 2: ConjGrad relax; 44: improved dimer method
NSW = 500 ! Number ionic steps
ISIF = 0 # relax ions, no change cell shape, no change cell volume (volume changes require ENCUT*1.3)
POTIM = 0.2

parallel:
LPLANE = .TRUE.
NPAR = 1

POSCAR (relaxation)

H
1.00000000000000
20.000 0.000 0.000
0.000 20.000 0.000
0.000 0.000 20.000
H
2
Selective dynamics
Cartesian
0.00000 0.00000 0.00000 T T T
0.00000 0.00000 0.80000 T T T

What do the calculated gradients and Hessian look like?

Hi,


sorry for the late reply. This is what I get from the geometry relaxation, in terms of the forces

POSITION TOTAL-FORCE (eV/Angst)
-----------------------------------------------------------------------------------
0.00000 0.00000 0.02510 0.000000 0.000000 0.001442
0.00000 -0.00000 0.77490 0.000000 -0.000000 -0.001442
-----------------------------------------------------------------------------------
total drift: -0.000000 -0.000000 0.000000

Seems to be ok.

This is what it looks like in the phonon calculation, and this certainly doesn't look alright.

POSITION TOTAL-FORCE (eV/Angst)
-----------------------------------------------------------------------------------
0.00000 0.00000 0.02510 0.000000 0.000000 -3473.078382
0.00000 0.00000 0.75490 0.000000 -0.000000 3473.078382
-----------------------------------------------------------------------------------
total drift: -0.000003 -0.000055 0.025066

Likewise, the Hessian looks suspicious too:

SECOND DERIVATIVES (NOT SYMMETRIZED)
------------------------------------
1X 1Y 1Z 2X 2Y 2Z
1X 4199.843257 0.000000 -0.000368-4199.843257 0.000000 0.000368
1Y 0.000000 4199.775546 -0.004923 0.000000-4199.775546 0.004923
1Z 0.000000 -0.000000************ 0.000000 0.00000015249.984065
2X -4199.633797 0.000000 -0.000301 4199.633797 0.000000 0.000301
2Y 0.000000-4199.619188 -0.003569 0.000000 4199.619188 0.003569
2Z 0.000000 -0.00000015248.481182 0.000000 0.000000************

These numbers are far too large. Could you have a look at my inputs as well, to check that I am not making some random mistake. Btw, I have repeated these calcs in the meantime, with somewhat tighter convergence criteria for the geometry optimisation. Same result. This is puzzling.

Thanks

Tobias

Hi,

I just wanted to mention that the problem has been solved. A cluster inspection of the phonon output revealed that the energies were completely wrong due to the use of wrong pseudo-potentials. It was after all a mistake on my side since I didn't check the POTCAR file. I was essentially trying to describe H atoms with Rh potentials. Upon repeating the calculation with the correct potentials for H, I get a reasonable result.

Tobi

please consider how many degrees of freedom a 2-atomic molecule has
(basic physical chemistry)....
the imaginary nodes are nothing else but the drift of the molecule in the vacuum.