## 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.xsf","dynmat.mold", in which vibrational modes are stored, which can be visualized by using XCrySDen and/or Molden. |