Water molecule with PHONON
Tutorials
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 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
&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 using the input file above 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 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.