Water molecule with PHONON - Francesco Buonocore CMAST Website 3.0

Francesco Buonocore
ENEA - TERIN-ICT Division, HPC Lab
F. Buonocore
Go to content

Water molecule with PHONON


Water vibrational frequencies calculation using phonon (ph.x in the QUANTUM ESPRESSO package)

In the following tutorial it is shown how to calculate frequencies of the vibrational modes of H2O with Pwscf & Phonon codes included in the QUANTUM ESPRESSO package installed on CRESCO.
The details on keywords are found in the QUANTUM ESPRESSO web-page and ph.x input file description.
1) Build up the H2O model with your preferred atomistic model editor.
2) Make a relax of H2O following the example of the relax calculation in the graphene tutorial. Here we provide the relaxed atomic coordinates and unit cell parameters at the following step 3.
To run a job on CRESCO6 the command line is:
bsub -n NP -q QUEUE -o lsf_file.out -e file.err qe 6.3.0 par EXE input.inp output.out
where NP is the number of processors (NP=48 is one node on CRESCO6), QUEUE=cresco6_48h24 (24h queue), cresco6_48h48 (48h queue, NP >= 48), cresco6_h144 (144h queue, NP can be less than 48), EXE = executable (e.g., pw, pp, bands, projwfc, ph).

3) Create your pwscf input file H2O.scf.in as follows (XCrySDen visualization is shown on the right)):

  ibrav = 12,
  a = 20.000,
  b = 20.000,
  c= 20.000,
  nat = 3,
  ntyp = 2,
  ecutwfc = 60.0 ,
  ecutrho = 480.0 ,
  smearing = 'marzari-vanderbilt' ,
  occupations = 'smearing' ,
  degauss = 0.001,
  electron_maxstep = 100,
  conv_thr = 1.0d-12 ,
  mixing_mode = 'plain' ,
  mixing_beta = 0.3d0 ,
       O    15.999 O.pz-rrkjus.UPF
       H     1.0079 H.pz-rrkjus.UPF
O        2.525856153   1.482855778   8.928868428
H        2.241743628   2.141988100   9.587993065
H        3.499400219   1.513156121   8.959138507
K_POINTS automatic
1 1 1 0 0 0

i) The convergence threshold of SCF calculation conv_thr must be very low (1.0D-12), smaller than for total energy calculation.
ii) The k-point cannot be set to gamma, because ph.x cannot use the gamma point algorithms (tricks). You chose to artificially set a mesh with one only k-point and no shift (1 1 1 0 0 0).
3) Launch the scf calculation using the input file above with EXE= pw.

4) Create your phonon input file H2O.ph.in as follows:

phonons of H2O in super cell at Gamma ! title line
tr2_ph=1.0d-14, ! the threshold
prefix='H2O', ! must be the same than scf calculation
amass(1)=15.999, ! If not specified, masses are read from data file
amass(2)=1.0079, ! If not specified, masses are read from data file
outdir='./your/scratch_dir/', ! must be the same than scf calculation
fildyn='H2O.dyn', ! name of output file
trans= .true. ! the phonons are computed

5) Launch the phonon calculation using the input file above with EXE= ph .
6) Create the H2O.dynmat.in input file for dynmat.x:

fildyn='H2O.dyn', ! dynamical matrix produced by ph.x

7) Impose that the acoustic sum rules (asr) on the elements of the dynamical matrix are satisfied and diagonalize the matrix with the code dynmat.x using the input file above (EXE=dynmat).

8) The diagonalized dynamical matrix is found in the output file dynmat.out, with the vibrational frequencies expressed in THz and cm-1 (that are also found in H2O.dynmat.out). The phonon eigendisplacements are found in dynmat.mold (Molden Format) and dynmat.axsf (XCrySDen format). We extract the following table from H2O.dynmat.out:

# mode   [cm-1]    [THz]      IR
   1   -221.53   -6.6412    0.0000
   2    -52.08   -1.5614    0.0000
   3     -8.84   -0.2650    0.0000
   4      9.60    0.2879    0.0000
   5     92.89    2.7849    0.0000
   6    173.99    5.2161    0.0000
   7   1547.07   46.3800    0.0000
   8   3795.94  113.7995    0.0000
   9   3914.63  117.3577    0.0000

Because we have N=3 atoms in the system, we have 3N-6= 3 vibrational degrees of freedom. The first six frequencies are related to translations and rotations. Compare the theoretical vibrational frequencies to the experiments:

# mode                             [cm-1]          exp [cm-1]            err
H-O-H bend                      1547.07        1711                     -9.5%
H-O symmetric strech        3795.94           3730                       1.7%
H-O symmetric strech         3914.63          3851                       1.6%

We observe a good agreement with experiments, but the H-O-H bending mode suffers a little underestimation.

Casaccia Research Center, Rome, Italy
TERIN-ICT Division, HPC Lab
Back to content
Cookies Policy