[PYTHON] Perform DFT calculation with ASE and GPAW

Let's perform DFT calculation using ASE and GPAW. As a simple example, let's find the atomization energy of a hydrogen molecule.

--ASE is a Python module for atomic simulation. Many features are implemented. Various codes can be called using the Calculator class as an interface. --GPAW is a DFT calculation code based on the PAW method and ASE. Written in C and Python. Widely used.

The following example was executed using Jupyter Lab.

Installation

Install ASE and GPAW with pip in the Python environment built with Anaconda.

$ pip install ase --user
$ pip install gpaw --user

Also download the GPAW PAW dataset.

$ gpaw install-data $HOME/gpaw-data

This is the end of installation. At this stage, my environment is as follows.

$ python --version
Python 3.7.6
$ pip list installed | grep -e ase -e gpaw
ase                                3.19.1
gpaw                               20.1.0

GPAW calculation conditions

Instantiate the GPAW class under the following conditions.

--Cutoff energy: 300eV --Function: PBE

from gpaw import GPAW, PW

gpaw = GPAW(mode=PW(300), xc='PBE')

Create a model of hydrogen molecule

The ʻase.build module contains small molecule presets. You can create a ʻAtoms object, which is a three-dimensional molecular model, simply by passing the molecular formula to the molecule function. Since I want to calculate the isolated system this time, I call the center method of the created ʻAtoms` object and add a 3.0A vacuum layer to the model.

from ase.build import molecule

atoms_h2 = molecule('H2')
atoms_h2.center(vacuum=3.)

It can be visualized with the view function, so let's check the structure of the molecule.

from ase.visualize import view

view(atoms_h2)

Perform the calculation

If you associate the GPAW Calculator with the ʻAtoms object and call the get_potential_energy` function, the GPAW calculation will be executed.

atoms_h2.set_calculator(gpaw)
e_h2 = atoms_h2.get_potential_energy() # takes several seconds

Calculation of hydrogen atom

The same calculation is performed for the hydrogen atom.

atoms_h = molecule('H')
atoms_h.center(vacuum=3.)
atoms_h.set_calculator(GPAW(mode=PW(300), xc='PBE', hund=True))
e_h = atoms_h.get_potential_energy()

View results

Finally, let's summarize the results. The atomization energy $ \ Delta E $ of H _{2 </ sub> can be obtained by the following equation.}

\Delta E = 2 E_{\rm H} - E_{\rm H_2}

At the same time, I will display the molecular structure on the Notebook. Use the ʻase.visualize.plot` module, which is a wrapper for matplotlib.

import matplotlib.pyplot as plt
from ase.visualize.plot import plot_atoms

delta_e = 2 * e_h - e_h2
print(f'Atomization Energy: {delta_e: 5.2f} eV')

fig, ax = plt.subplots(1, 2, figsize=(8,6))

title = f'Calculation Result for {atoms_h2.get_chemical_formula()}\n' + \
        f'Total Energy : {atoms_h2.get_potential_energy(): 5.2f} eV'
ax[0].set_title(title)
ax[0].axis('off')
plot_atoms(atoms_h2, ax=ax[0], rotation='90x')

title = f'Calculation Result for {atoms_h.get_chemical_formula()}\n' + \
        f'Total Energy : {atoms_h.get_potential_energy(): 5.2f} eV'
ax[1].set_title(title)
ax[1].axis('off')
plot_atoms(atoms_h, ax=ax[1], rotation='90x')

plt.show()

reference

Atomic Simulation Environment https://wiki.fysik.dtu.dk/ase/
GPAW https://wiki.fysik.dtu.dk/gpaw/