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