Oxygen molecule †

In this tutorial, how to perform calculation(s) of a molecule with the spin degree of freedom (or spin unrestricted calculation of a molecule) is discussed by taking an oxygen (O2) molecule as an example.

↑

Preliminary †

In this tutorial, the O2 molecule is placed in a cubic box with the cell edges of 10 Angstroms. The GBRV ultrasoft pseudopotential, kinetic energy cutoffs of 40 Ry and 320 Ry for the wave functions and augmentation charge, respectively, and the Gamma-point are used.

↑

SCF calculation †

↑

Fixed moment calculation †

We start with an SCF calculation using the following input file. The bond length is set to 1.2 Angstrom.

  • o2_scf.in 
    O2
O2 molecule in a box
 &control
    calculation = 'scf'
    restart_mode = 'from_scratch'
    prefix = 'o2'
    pseudo_dir = '/home/ikutaro/QE/pseudo'
    outdir = './tmp/'
    verbosity = 'high'
 /
 &system
    ibrav = 1
    A = 10.d0
    nat = 2
    ntyp = 1
    nbnd = 12
    nosym = .true.
    nspin = 2
    tot_magnetization = 2.0
    starting_magnetization(1) = 0.5d0
    ecutwfc = 40.0
    ecutrho = 320.0
 /
 &electrons
    mixing_beta = 0.5
    conv_thr =  1.0d-10
 /
 &ions
 /
ATOMIC_SPECIES
O    15.9994   o_pbe_v1.2.uspp.F.UPF
ATOMIC_POSITIONS (angstrom)
 O      0.000000000000      0.000000000000      0.000000000000
 O      0.000000000000      0.000000000000      1.200000000000
K_POINTS (gamma)
    We also use a job script like following:
  • o2.sh 
    #$ -S /bin/bash
#$ -cwd
#$ -q xe2.q
#$ -pe x12 12
#$ -N O2
#$ -o O2.out
#$ -e O2.err
#
module load intel/2020.2.254 intelmpi/2020.2.254
#
INPUT_FILE='o2_scf.in'
OUTPUT_FILE='o2_scf.out'
#
PW_DIR=/home/ikutaro/QE/src/qe-7.2/bin
#
MPI_COMMAND=mpirun
PW_COMMAND=pw.x
PW=$PW_DIR/$PW_COMMAND
#
I_MPI_PIN=1
I_MPI_FABRICS=shm:dapl
OMP_NUM_THREADS=1
#
$MPI_COMMAND $PW < $INPUT_FILE > $OUTPUT_FILE

To submit the job execute

$ qsub o2.sh

In this example, we perform the fixed moment calculation to obtain the ground state triplet state of the O2 molecule (with the magnetic moment of 2.0 Bohr magneton) by setting

   tot_magnetization = 2.0

in the namelist &SYSTEM. Note that nonzero initial magnetization should be given to obtain spin polarized solution as

    starting_magnetization(1) = 0.5d0

When the convergence of the SCF calculation is achieved, we can confirm that the desired magnetization is obtained by assessing the output file (o2_scf.out) as:

     The total energy is the sum of the following terms:
     one-electron contribution =    -124.42909406 Ry
     hartree contribution      =      64.11909734 Ry
     xc contribution           =     -14.24282541 Ry
     ewald contribution        =      10.24736697 Ry

     total magnetization       =     2.00 Bohr mag/cell
     absolute magnetization    =     2.05 Bohr mag/cell

and electron configuration as (by setting verbosity = 'high'):

 ------ SPIN UP ------------


          k = 0.0000 0.0000 0.0000 ( 14336 PWs)   bands (ev):

   -32.8993 -20.2605 -13.4590 -13.4590 -13.4215  -6.5703  -6.5703  -0.4687
     0.8418   0.9710   0.9732   0.9749

     occupation numbers 
     1.0000   1.0000   1.0000   1.0000   1.0000   1.0000   1.0000   0.0000
     0.0000   0.0000   0.0000   0.0000

 ------ SPIN DOWN ----------


          k = 0.0000 0.0000 0.0000 ( 14336 PWs)   bands (ev):

   -31.7149 -18.4058 -12.4328 -11.5313 -11.5313  -4.1646  -4.1646  -0.3714
     1.0168   1.0686   1.0826   1.0826

     occupation numbers 
     1.0000   1.0000   1.0000   1.0000   1.0000   0.0000   0.0000   0.0000
     0.0000   0.0000   0.0000   0.0000

     highest occupied, lowest unoccupied level (ev):    -6.5703   -4.1646
↑

Fixed occupation calculation †

Alternatively, we can also perform calculations with fixed (desired) occupation number by using input file like:

  • o2_scf.in 
    O2
O2 molecule in a box
 &control
    calculation = 'scf'
    restart_mode = 'from_scratch'
    prefix = 'o2'
    pseudo_dir = '/home/ikutaro/QE/pseudo'
    outdir = './tmp/'
    verbosity = 'high'
 /
 &system
    ibrav = 1
    A = 10.d0
    nat  =  2
    ntyp = 1
    nbnd =12
    nosym = .true.
    nspin = 2
    occupations = 'from_input'
    starting_magnetization(1) = 0.5d0
    ecutwfc = 40.0
    ecutrho = 320.0
 /
 &electrons
    mixing_beta = 0.5
    conv_thr =  1.0d-10
 /
 &ions
 /
ATOMIC_SPECIES
O    15.9994   o_pbe_v1.2.uspp.F.UPF
ATOMIC_POSITIONS (angstrom)
 O      0.000000000000      0.000000000000      0.000000000000
 O      0.000000000000      0.000000000000      1.200000000000
K_POINTS (gamma)
OCCUPATIONS
1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0
1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
    The fixed moment calculation is enabled setting 
        occupations = 'from_input'
    and OCCUPATIONS 
    1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0
1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
    Note that the number of states (bands) should be consistent with nbnd in the namelist &SYSTEM.
↑

Structural relaxation †

Structural optimization can be performed by setting the calculation='relax' as

  • o2_relax.in 
    O2
O2 molecule in a box
 &control
    calculation = 'relax'
    restart_mode = 'from_scratch'
    prefix = 'o2'
    pseudo_dir = '/home/ikutaro/QE/pseudo'
    outdir = './tmp/'
    verbosity = 'high'
 /
 &system
    ibrav = 1
    A = 10.d0
    nat = 2
    ntyp = 1
    nbnd = 12
    nosym = .true.
    nspin = 2
    tot_magnetization = 2.0
    starting_magnetization(1) = 0.5d0
    ecutwfc = 40.0
    ecutrho = 320.0
 /
 &electrons
    mixing_beta = 0.5
    conv_thr =  1.0d-10
 /
 &ions
 /
ATOMIC_SPECIES
O    15.9994   o_pbe_v1.2.uspp.F.UPF
ATOMIC_POSITIONS (angstrom)
 O      0.000000000000      0.000000000000      0.000000000000
 O      0.000000000000      0.000000000000      1.200000000000
K_POINTS (gamma)
    To perform the structural relaxation 
     &ions
 /
    needs to be added in the input file. In this example, the default force threshold (forc_conv_thr = 1.e-3) and relaxation method (ion_dynamics = 'bfgs') are used.
↑

Calculation of the atomization energy †

Let us calculate the atomization energy, which is the energy gain to form the molecule from isolated atoms, defined as

\[ E_{\rm{a}} = E_{\rm{tot}}(\rm{O}_2) - 2 \times E_{\rm{tot}}(\rm{O}) \]

where \( E_{\rm{tot}}(\rm{O}_2) \) and \( E_{\rm{tot}}(\rm{O}) \) are total energies of (optimized) O2 molecule and O atom, respectively, and the total energy of the O atom should be that of the spin-polarized ground state. To perform the calculation of the spin polarized O atom in the ground state, an input file like following may be used:

  • o_scf.in 
    O
O atom in a box
 &control
    calculation = 'scf'
    restart_mode = 'from_scratch'
    prefix = 'o'
    pseudo_dir = '/home/ikutaro/QE/pseudo'
    outdir = './tmp/'
    verbosity = 'high'
 /
 &system
    ibrav = 1
    A = 10.d0
    nat  =  1
    ntyp = 1
    nbnd = 8
    nosym = .true.
    nspin = 2
    occupations = 'from_input'
    one_atom_occupations = .true.
    starting_magnetization(1) = 0.5d0
    ecutwfc = 40.0
    ecutrho = 320.0
 /
 &electrons
    mixing_beta = 0.1
    conv_thr =  1.0d-10
 /
 &ions
 /
ATOMIC_SPECIES
O    15.9994   o_pbe_v1.2.uspp.F.UPF
ATOMIC_POSITIONS (angstrom)
 O      0.000000000000      0.000000000000      0.000000000000
K_POINTS (gamma)
OCCUPATIONS
1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0
1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
    From my experience, it is better to use the fixed occupations with 
        one_atom_occupations = .true.
    Furthermore, the SCF convergence of the atomic calculation can be ridiculously difficult. In this case, try to use small mixing_beta and increase electron_maxstep (from default value (100) to 200-400).

After the SCF convrgence, we can calculate the binding energy, and I obtained -5.93 eV. Compare with the values in the literature, for example, Hammer et al., Phys. Rev. B 59, 7413 (1999).

↑

Exercise †

The atomization energy, in particular the oxygen atom, is sensitive to the computational setup. Try to use different functional (RPBE, for instance) and different (pseudo)potentials (normconserving and PAW) and see the differences. It is also important to converge the results with respect to the plane-wave cutoffs.

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