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)):

&CONTROL

calculation='scf'

title='H2O'

prefix='H2O'

verbosity='high'

restart_mode='from_scratch'

nstep=1000

iprint=1

tstress=.true.

tprnfor=.true.

disk_io='default'

pseudo_dir='/afs/enea.it/software/qu_esp/pseudo/'

outdir='./your/scratch_dir/',

tstress=.true.

tprnfor=.true.

forc_conv_thr=1.0d-4

etot_conv_thr=1.0d-5

&SYSTEM

ibrav = 12,

a = 20.000,

b = 20.000,

c= 20.000,

cosab=-0.500000,

nat = 3,

ntyp = 2,

ecutwfc = 60.0 ,

ecutrho = 480.0 ,

smearing = 'marzari-vanderbilt' ,

nspin=1,

occupations = 'smearing' ,

degauss = 0.001,

&ELECTRONS

electron_maxstep = 100,

conv_thr = 1.0d-12 ,

mixing_mode = 'plain' ,

mixing_beta = 0.3d0 ,

&IONS

ion_dynamics='bfgs'

upscale=20.0

ATOMIC_SPECIES

O 15.999 O.pz-rrkjus.UPF

H 1.0079 H.pz-rrkjus.UPF

ATOMIC_POSITIONS angstrom

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

Notes:

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 with EXE= pw.

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

phonons of H2O in super cell at Gamma ! title line

&inputph

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 with EXE= ph .

6) Create the H2O.dynmat.in input file for dynmat.x:

&input

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

asr='simple',

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.

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.

12) The zero-point energy associated to vibrations can be easily evaluated, as: