#author("2023-10-04T11:16:47+09:00","default:StatE","StatE")
#author("2023-10-04T11:18:40+09:00","default:StatE","StatE")
* Oxygen molecule [#k112b541]
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.
#Contents
** Preliminary [#bc970db3]
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 [#b9138eba]
*** Fixed moment calculation [#xe714aec]
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 [#ucd1be8e]
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 [#s024b02b]
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 [#f68e7f56]
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 [#dfacc88e]
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 (O.pbe-mt.UPF, for instance) and PAW (O.pbe-n-kjpaw_psl.1.0.0.UPF, for instance)) and see the differences. It is also important to converge the results with respect to the plane-wave cutoffs for each pseudopotentials used.
The atomization energy, in particular the oxygen atom, is sensitive to the computational setup. Try to use
- different functional (RPBE, BEEF-vdW, for instance)
- different pseudo)potentials
-- Normconserving (O.pbe-mt.UPF, for instance)
-- PAW (O.pbe-n-kjpaw_psl.1.0.0.UPF, for instance))
and see the differences.~
It is also important to converge the results with respect to the plane-wave cutoffs for each pseudopotentials used.
![[PukiWiki]](image/logo_state-senri3_80x80.png "[PukiWiki]"){kind=logo}
