Hydrogen molecule

In this tutorial, how to execute pw.x and how to converge computational parameters, such as kinetic energy cutoff and box size are described. How to optimize the molecular geometry is also described.

Preliminary

Let us begin with H2, the simplest molecule. In a standard plane-wave pseudopotential code, three dimensional periodic boundary condition is employed, and an isolated molecule is modeled by placing it in a large box in such a way that the electrostatic interactions with neighboring images can be negligible. If the box is large enough, the Brillouin zone becomes small (as it scales as (2p/a)**3), thus the Brillouin zone is sampled only at the Gamma point [(0,0,0)]. In this case, the accuracy of the calculation is controlled by increasing the kinetic energy cutoff, which corresponds to increasing the real-space mesh (resolution of the wave function in real space) and by increasing the box size. The latter is to confirm that the interaction between image molecules is negligibly small and that the model mimics the molecule in itsan isolated state.

SCF calculation

Let us begin with an SCF calculation of H2 having the experimental bond length of 0.74 Angstrom in a cubic box with the edges of 10 Angstrom. We use the SG15 optimized norm-conserving pseudopotentials with the cutoff energy of 30 Ry. Here the input file for this calculation is called "h2_scf.in" and looks as follows:

&control
   calculation   = 'scf'
   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       = 10.d0
   nat     = 2
   ntyp    = 1
   nbnd    = 4
   ecutwfc = 30.0
   ecutrho = 120.0
/
&electrons
   diagonalization = 'david'
   conv_thr        = 1.0e-12
   mixing_beta     = 0.1
/
ATOMIC_SPECIES
 H 0.0000  H_ONCV_PBE-1.0.UPF
ATOMIC_POSITIONS (angstrom)
H      0.000000000000      0.000000000000      0.00000000000
H      0.000000000000      0.000000000000      0.74000000000
K_POINTS (gamma)

Note "pseudo_dir" and "outdir" should be modified according to the environment. To perform the SCF calculation, we use pw.x. Here let us assume we are going to use pw.x in $HOME/QE/src/qe-6.2.1 and execute the following

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

When the SCF convergence is archived, the following can be found in the output file (h2_scf.out).

!    total energy              =      -2.31431279 Ry
     Harris-Foulkes estimate   =      -2.31431279 Ry
     estimated scf accuracy    <          3.7E-13 Ry

     The total energy is the sum of the following terms:

     one-electron contribution =      -3.79278238 Ry
     hartree contribution      =       2.01958218 Ry
     xc contribution           =      -1.37197167 Ry
     ewald contribution        =       0.83085908 Ry

     convergence has been achieved in  15 iterations

If "tprnfor = .true." the following lines can also be found

     Forces acting on atoms (cartesian axes, Ry/au):

     atom    1 type  1   force =     0.00000000    0.00000000   -0.03012765
     atom    2 type  1   force =     0.00000000    0.00000000    0.03012765

Convergence study

Convergence wrt cutoff

To study the convergence of the total energy, we calculate the total energy as a function of kinetic energy cutoff ("ecutwfc") by fixing the molecular structure. Note that in the present case normconserving pseudopotential is used and ecutrho is fixed to 4 times larger than ecutwfc (ecutrho = 4 * ecutwfc). Results are as follows

# ecutwfc (Ry) total energy (Ry)
  30.0         -2.31431279
  40.0         -2.32593152
  50.0         -2.33045981
  60.0         -2.33201780
  70.0         -2.33257560
  80.0         -2.33275476
  90.0         -2.33280832
  100.0        -2.33281717

Convergence wrt box size

Then let us calculate the total energy as a function of the cell edge. In this study, ecutwfc=60.0 Ry is used.

# a (Angstrom) total energy (Ry)
5.0            -2.33788501
10.0           -2.33201780
15.0           -2.33202479
20.0           -2.33202183
25.0           -2.33202063
30.0           -2.33202053

Optimizing the bond length (1)

Optimizing the bond length (2)

トップ   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS