Gas formation energy calculation using DFT#
Gas formation Energy#
Formation energies describe the energy required for (or released by) the formation of a chemical species from some reference species at a reference state.
You may notice that “raw” DFT energies are typically quite large for molecules and surfaces. This is because the DFT energy is calculated relative to a much higher energy reference state. Calculating formation energies allows us to remove the often inconvenient DFT reference and specify one of our own.
The gas formation energy, \({\Delta}E_{f}\), for water can be obtained by calculating the difference between the energy of the formed gas molecule and the total energies of its components. For example, for \( \text{H}_{2}\text{O} \):
\({\Delta}E_{f} (H_{2}O) = E_{H_{2}O} - E_{H_{2}} - 1/2 \cdot E_{O_{2}}\)Here, we have specified two references (H2,O2) which is equal to the number of unique elements in H2O
Calculating gas formation Energy in SPARC#
H2O formation energy
Building Atoms
In the following script, we will make Atoms objects for H2O, H2, and O2 molecules.
import numpy as np
from ase.build import molecule, bulk
#building O2 molecule
O2 = molecule('O2')
O2.set_cell(np.array([10, 10, 10]))
O2.center()
#building H2 molecule
H2 = molecule('H2')
H2.set_cell(np.array([10, 10, 10]))
H2.center()
#building H2O molecule
H2O = molecule('H2O')
H2O.set_cell(np.array([10, 10, 10]))
H2O.center()
DFT calculation
In the following script, we define a SPARC calculator object and attach it to an ase.Atoms object so that ase knows how to get the data. After defining the calculator, a DFT calculation can be run. The potential energy can be obtained by the method: get_potential_energy(). The potential energy result can be stored in a text file e.g. ‘energy.txt’.
Details on how to submit these calculations can be seen in Lecture ‘Running_SPARC_on_PACE’.
O2 calculation
from ase.optimize import BFGSLineSearch
from ase.units import Bohr,Hartree,mol,kcal,kJ,eV
from sparc import SPARC
import numpy as np
image = O2
parameters = dict(
EXCHANGE_CORRELATION = 'GGA_PBE',
D3_FLAG=1, #Grimme D3 dispersion correction
SPIN_TYP=1, #spin-polarized calculation
KPOINT_GRID=[1,1,1], #molecule needs single kpoint !
ECUT=800/Hartree, #set ECUT (Hartree) or h (Angstrom)
#h = 0.15,
TOL_SCF=1e-5,
RELAX_FLAG=1, #Do structural relaxation (only atomic positions)
TOL_RELAX = 5.00E-04, #convergence criteria (maximum force) (Ha/Bohr)
PRINT_FORCES=1,
PRINT_RELAXOUT=1)
parameters['directory'] = 'O2'
calc = SPARC(atoms = image, **parameters)
image.set_calculator(calc)
eng = image.get_potential_energy()
with open('O2_energy.txt', 'w') as f:
f.write(str(eng))
Successfully registered sparc format with ase.io!
/home/hice1/nyu49/.local/lib/python3.9/site-packages/sparc/calculator.py:471: UserWarning: You have specified one of FD_GRID, ECUT or MESH_SPACING, conversion of h to mesh grid is ignored.
warn("You have specified one of FD_GRID, ECUT or MESH_SPACING, "
Step 0
Step 1
Step 2
Step 3
Step 4
{}
[0, 1] [0, 1]
/home/hice1/nyu49/.local/lib/python3.9/site-packages/sparc/sparc_parsers/utils.py:56: UserWarning: Key POISSON_SOLVER not in validator's parameter list, ignore value conversion!
warn(
/home/hice1/nyu49/.local/lib/python3.9/site-packages/sparc/sparc_parsers/utils.py:56: UserWarning: Key VERBOSITY not in validator's parameter list, ignore value conversion!
warn(
H2 calculation
from ase.optimize import BFGSLineSearch
from ase.units import Bohr,Hartree,mol,kcal,kJ,eV
from sparc import SPARC
import numpy as np
image = H2
parameters = dict(
EXCHANGE_CORRELATION = 'GGA_PBE',
D3_FLAG=1, #Grimme D3 dispersion correction
SPIN_TYP=1, #spin-polarized calculation
KPOINT_GRID=[1,1,1], #molecule needs single kpoint !
ECUT=800/Hartree, #set ECUT (Hartree) or h (Angstrom)
#h = 0.15,
TOL_SCF=1e-5,
RELAX_FLAG=1, #Do structural relaxation (only atomic positions)
TOL_RELAX = 5.00E-04, #convergence criteria (maximum force) (Ha/Bohr)
PRINT_FORCES=1,
PRINT_RELAXOUT=1)
parameters['directory'] = 'H2'
calc = SPARC(atoms = image, **parameters)
image.set_calculator(calc)
eng = image.get_potential_energy()
with open('H2_energy.txt', 'w') as f:
f.write(str(eng))
Step 0
Step 1
Step 2
Step 3
{}
[0, 1] [0, 1]
H2O calculation
from ase.optimize import BFGSLineSearch
from ase.units import Bohr,Hartree,mol,kcal,kJ,eV
from sparc import SPARC
import numpy as np
image = H2O
parameters = dict(
EXCHANGE_CORRELATION = 'GGA_PBE',
D3_FLAG=1, #Grimme D3 dispersion correction
SPIN_TYP=1, #spin-polarized calculation
KPOINT_GRID=[1,1,1], #molecule needs single kpoint !
ECUT=800/Hartree, #set ECUT (Hartree) or h (Angstrom)
#h = 0.15,
TOL_SCF=1e-5,
RELAX_FLAG=1, #Do structural relaxation (only atomic positions)
TOL_RELAX = 5.00E-04, #convergence criteria (maximum force) (Ha/Bohr)
PRINT_FORCES=1,
PRINT_RELAXOUT=1)
parameters['directory'] = 'H2O'
calc = SPARC(atoms = image, **parameters)
image.set_calculator(calc)
eng = image.get_potential_energy()
with open('H2O_energy.txt', 'w') as f:
f.write(str(eng))
Step 0
Step 1
Step 2
{}
[2, 0, 1] [1, 2, 0]
/home/hice1/nyu49/.local/lib/python3.9/site-packages/sparc/sparc_parsers/out.py:108: UserWarning: Key atomic mass from run information appears multiple times in your outputfile!
warn(
Finally, obtain H2O formation energy
import os
from ase.io import read
def readFile(path):
f = open(path, 'r')
content = f.readlines()
return content
def read_energy(path):
energy = float(readFile(path)[0])
return energy
E_H2O = read_energy('H2O_energy.txt')
E_H2 = read_energy('H2_energy.txt')
E_O2 = read_energy('O2_energy.txt')
E_f = E_H2O - E_H2 - 1/2*E_O2
print('The formation energy of H2O is {0:1.3f} eV.'.format(E_f))
print("The obtained value (-2.462 eV) is quite close to the experimental value of -2.712 eV.")
The formation energy of H2O is -2.462 eV.
The obtained value (-2.462 eV) is quite close to the experimental value of -2.712 eV.