第一原理分子動力学プログラム STATE Senri Wiki

Vibrational mode analysis of a hydrogen molecule †

In this tutorial, how to calculate vibrational frequencies of a molecule by density functional perturbation theory as implemented in PH (ph.x) by taking a hydrogen molecule (H2) in a large box as an example.

Preliminary †

in the parent directory build ph.x, if not, by typing

 $ make ph

Geometry optimization †

To calculate vibrational frequencies accurately, we should use fully optimized geometry with a tight convergence threshold. Following is an example input file for the relaxation of the hydrogen molecule. Here the molecule is put in a cubic box with the cell edges of 15 Angstroms.

 &control
  calculation  = 'relax'
  restart_mode = 'from_scratch',
  pseudo_dir   = '/home/ikutaro/QE/pseudo'
  outdir       = './tmp'
  prefix       = 'h2'
  tprnfor      = .true.
  forc_conv_thr = 1.d-4
 /
 &system
  ibrav   = 1
  A       = 15.d0
  nat     = 2
  ntyp    = 1
  ecutwfc = 50.0
  ecutrho = 400.0
  nbnd    = 8
 /
 &electrons
 diagonalization='david'
 conv_thr = 1.0e-12
 mixing_beta = 0.7
 /
 &ions
 /
ATOMIC_SPECIES
H     1.00794  h_pbe_v1.4.uspp.F.UPF
ATOMIC_POSITIONS (angstrom)
H  0.00000000  0.00000000  0.00000000
H  0.00000000  0.00000000  0.75000000
K_POINTS (gamma)

Then execute pw.x by executing the following (and/or alike):

$ mpirun -np 2 $HOME/QE/src/qe-6.2.1/bin/pw.x < h2_relax.in > h2_relax.out

SCF calculation †

We then perform a single-shot SCF calculation by using the optimized geometry, to obtain the unperturbed wave functions for the perturbation calculations. For this purpose, I recommend using a very tight convergence threshold for SCF such as 1.e-12 - 1.e-14. Following is an example of the input file.

 &control
  calculation  = 'scf'
  restart_mode = 'from_scratch',
  pseudo_dir   = '/home/k0227/k022700/QE/pseudo'
  outdir       = './tmp'
  prefix       = 'h2'
  tprnfor      = .true.
 /
 &system
  ibrav   = 1
  A       = 15.d0
  nat     = 2
  ntyp    = 1
  ecutwfc = 50.0
  ecutrho = 400.0
  nbnd    = 8
 /
 &electrons
 diagonalization='david'
 conv_thr = 1.0e-14
 mixing_beta = 0.7
 /
ATOMIC_SPECIES
H     1.00794  h_pbe_v1.4.uspp.F.UPF
ATOMIC_POSITIONS (angstrom)
H        0.000000000   0.000000000  -0.005787349
H        0.000000000   0.000000000   0.755787349
K_POINTS (automatic)
  1  1  1  0  0  0

Note that ph.x is NOT able to read real wave functions obtained with "K_POINTS (gamma)", and we are forced to create and use complex wave functions with

K_POINTS (automatic)
  1  1  1  0  0  0

"phcg.x" may be able to handle real wave functions, but the code is limited to spin-unpolarized LDA.

Execute pw.x by executing the following (and/or alike):

$ mpirun -np 2 $HOME/QE/src/qe-6.2.1/bin/pw.x < h2_scf.in > h2_scf.out

Vibrational mode analysis †

Having obtained the optimized geometry and the unperturbed (complex) wave function, we are at the stage to calculate the vibrational frequencies. Below is an (minimal) input file for the ph.x program.

H2 molecule
 &inputph
  outdir       = './tmp'
  prefix       = 'h2'
  fildyn       = 'h2.dyn'
  epsil        = .true.
  tr2_ph       = 1.d-18
 !recover      = .true.
 /
0.0 0.0 0.0

Here "outdir" and "prefix" should be the same as those in the SCF calculation. "fildyn" is the name of file, in which dynamical matrix is stored. "epsil" allows one to compute the dielectric matrix. Convergence threshold "tr2_ph" should be as tight as 1.e-18. With "recover = .true.", one is able to restart the vibrational calculation from an interrupted run. For a gas-phase molecule, only the gamma-point is needed to be considered. To do this, the q-point (0.0 0.0 0.0) is explicitly specified after the &inputph.../ section. To run ph.x the following command is executed

$ mpirun -np 2 $HOME/QE/src/qe-6.4.1/bin/ph.x < h2_ph.in > h2_ph.out

and one get the vibrational frequencies in the output file as:

     Diagonalizing the dynamical matrix

     q = (    0.000000000   0.000000000   0.000000000 ) 

 **************************************************************************
     freq (    1) =      -3.582291 [THz] =    -119.492378 [cm-1]
     freq (    2) =      -3.582290 [THz] =    -119.492341 [cm-1]
     freq (    3) =      -3.552544 [THz] =    -118.500120 [cm-1]
     freq (    4) =      -3.552544 [THz] =    -118.500117 [cm-1]
     freq (    5) =      -2.974216 [THz] =     -99.209174 [cm-1]
     freq (    6) =     131.121910 [THz] =    4373.756130 [cm-1]
 **************************************************************************

Negative frequency indicates the imaginary mode. Note that the hydrogen molecule is a linear diatomic molecule and only the vibrational frequency of the stretching mode is meaningful, and the calculated vibrational frequency is in good agreement with the experimental harmonic vibrational frequency (4401 cm^-1 from NIST). In addition, by using the "dynmat.x" program, vibrational modes can be visualized: The following is the minimalist input file (dynmat.in):

&input
 fildyn='h2.dyn'
 filout='dynmat.out'
/

By executing the following

$ $HOME/QE/src/qe-6.4.1/bin/dynmat.x < dynmat.in > dynmat.log

we obtain "dynmat.axsf","dynmat.mold", in which vibrational modes are stored, which can be visualized by using XCrySDen and/or Molden. When XCrySDen is not available, copy ~ikutaro/QE/tools/axsf2xsf.sh by doing

$ cp ~ikutaro/QE/tools/axsf2xsf.sh .

and execute

$ ./axsf2xsf.sh dynmat.axsf

then we get dynmat_#*.xsf for each mode. These file can be visualized with VESTA.