ASE Calculator

In the previous section, we learned the basics about ASE Atoms and how to generate structures.

In this section, we will learn how to handle Calculator, which is necessary for physical simulations.

Energy

The total energy \(E\) of the system is expressed in terms of the kinetic energy \(K\) and the potential energy \(V\).

\[E = K + V\]

Of these, the kinetic energy \(K\) of atoms can be calculated as follows, from classical mechanics.

\[K = \sum_{i=1}^{N} \frac{1}{2} m_i {\mathbf{v}}_i^2 = \sum_{i=1}^{N} \frac{{\mathbf{p}}_i^2}{2 m_i}\]

where \(m_i, \mathbf{v}_i, \mathbf{p}_i\) are the mass, velocity and momentum (\(\mathbf{p}=m\mathbf{v}\)) of each atom, respectively.

Those in bold represent vectors, where the velocity \(\mathbf{v}\) and momentum \(\mathbf{p}\) have three components in xyz coordinates.

On the other hand, the potential energy \(V\) is not easy to obtain exactly because it requires solving equations derived from quantum mechanics (e.g., Schrödinger equation, see column). There are various methods for determining potential energy, ranging from classical force fields, which is fast but limited in its applicability, to quantum chemical calculations such as the density functional theory (DFT), which takes longer but is more accurate (The differences between these methods are discussed below).

So why do we want to know about energy in the first place? In another word, what can we learn from energy?

First, we can know the stable structure of a substance, which is not self-evident in nature (as explained in chapter 2). By looking at the low and high energies, we can analyze whether the structure is realized in nature or not.

Second, we can learn how each atom moves. This will be explained in the next section.

Force

In classical mechanics, Newton’s equations of motion are expressed as follows

\[\mathbf{F} = m\mathbf{a}\]

This equation represents that given a force \(\mathbf{F}\), we know the acceleration \(\mathbf{a}\) applied to the system.

The position \(\mathbf{r}\), velocity \(\mathbf{v}\), and acceleration \(\mathbf{a}\) of each atom are related as follows.

\[\mathbf{v} = \frac{d\mathbf{r}}{dt}\]
\[\mathbf{a} = \frac{d\mathbf{v}}{dt} = \frac{d^2\mathbf{r}}{dt^2}\]

Therefore, when we know the force, we know the acceleration (= how the velocity changes), and as a result, we know how the position evolves over time. MD (Molecular Dynamics), which actually deals with time evolution in this way, is covered in chapter 6.

The forces can be expressed as the derivative of the potential energy with resprect to positions.

\[\mathbf{F} = -\frac{\partial V}{\partial \mathbf{r}}\]

In summary, the force \(\mathbf{F}\), required to know the time evolution \(\mathbf{r}(t)\) of the system, can be calculated if the potential energy \(V\mathbf(r)\) is known.

[Note] For the sake of simplicity, this tutorial has been explained using only knowledge of classical mechanics. For those who have studied analytical mechanics, the governing equations are defined by the Hamiltonian. Knowing the Hamiltonian (the function related to the energy mentioned above) allows us to describe the time evolution of the system.

Calculator class

We have explained that once the potential energy \(V(\mathbf{r})\) is determined, we know the governing equations of system, that is, how the atoms move in this world.

In ASE, the Calculator class is in charge of calculating the potential energy, and the calculation method can be switched by switching Calculator.

Let’s look at one example. The following is a calculation of each energy of H2 using the PFP calcultor provided by Matlantis.

The detailed usage is described later, but you can calculate each energy by setting Calculator for atoms.calc. In this example, we use below methods to calculate each energy.

  • Total energy \(E\): atoms.get_total_energy()

  • Potential energy \(V\): atoms.get_potential_energy()

  • Kinetic energy \(K\): atoms.get_kinetic_energy()

[1]:

from ase import Atoms

import pfp_api_client
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode


# print(f"pfp_api_client: {pfp_api_client.__version__}")

estimator = Estimator()
calculator = ASECalculator(estimator)

[2]:

atoms = Atoms("H2", [[0, 0, 0], [1.0, 0, 0]])
atoms.set_momenta([[0.1, 0, 0], [-0.1, 0, 0]])
atoms.calc = calculator

E_tot = atoms.get_total_energy()
E_pot = atoms.get_potential_energy()
E_kin = atoms.get_kinetic_energy()

print(f"Total Energy     : {E_tot:f} eV")
print(f"Kinetic Energy   : {E_kin:f} eV")
print(f"Potential Energy : {E_pot:f} eV")

Total Energy     : -3.808640 eV
Kinetic Energy   : 0.009921 eV
Potential Energy : -3.818561 eV

Atoms initially have an initial velocity of 0 and kinetic energy of 0, so an appropriate initial velocity is set using the set_momenta function.

You can see that $ E = K + V $.

The Calculator is necessary to calculate the total energy \(E\) and potential energy \(V\), and without the atoms.calc = calculator setting, the calculation cannot be performed and an error will be raised.

Calculator type

ASE supports a variety of calculators, from classical force fields to calculators using quantum chamical calculations. Some examples are listed below. Many other calculators are supported in addition to those listed here. Please check Supported calculators for details.

Category

Calculator

ASE embedded

Description

Classical force field

lj

✓

Lennard-Jones potential

morse

✓

Morse potential

emt

✓

Effective Medium Theory calculator

lammps

Classical molecular dynamics code

Quantum chemical calculations

gaussian

Gaussian based electronic structure code

vasp

Plane-wave PAW code

espresso

Plane-wave pseudopotential code

NNP (Neural Network Potential)

PFP

Potential provided by Matlantis

Calculators with ✓ in the “ASE embedded” column above can be used simply by installing ASE, other calculators require to install external program.

Calculators that perform quantum chemical calculations theoretically calculate potentials by solving the Schrödinger equation derived from quantum mechanics under a certain approximation (see also the column below). Special software is required and some calculations are time-consuming. It is used when you want to perform more accurate calculations.

https://en.wikipedia.org/wiki/Lennard-Jones_potential

The calculators belonging to classical force fields are methods that calculate potentials using manually created (empirical) function formulas, e.g., Lennard-Jones potential, Morse potential. The calculation speed is fast, but each method can handle different structures and physical phenomena. You need validation to check if the potential works for your interest.

The NNP type is intended to perform highly accurate potential prediction in a fast computation time by preparing many data computed in advance using quantum chemical calculations, and by supervised learning of the relationship between input structures and their potential energies. For more information on NNP, please refer to the following slides as well.

ASE embedded calculator

Here we use EMT calculator as an example.

[3]:

from ase.build import bulk
from ase.calculators.emt import EMT


calculator_emt = EMT()

atoms = bulk("Cu")
atoms.calc = calculator_emt

E_pot_emt = atoms.get_potential_energy()
print(f"Potential energy {E_pot_emt:.5f} eV")

Potential energy -0.00568 eV

The calculator can be switched to perform atomistic simulations in this way.

If the exact potential energy \(V(\mathbf{r})\) could be calculated quickly, there would be no need to switch calculators and the simulation could always be run using that, but in reality the exact potential energy cannot be calculated, and there is a tradeoff between accuracy and speed. Therefore, the choice must be made by the user of the atomistic simulation to suit the application.

PFP

The PFP provided by Matlantis is characterized by universality and high speed. It can be used to perform fast calculations for systems with a variety of structures in combinations of 55 elements.

PFP calculator is usually used in this tutorial throughout.

[4]:

import pfp_api_client
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode

estimator = Estimator(model_version="v2.0.0")
calculator = ASECalculator(estimator)

atoms.calc = calculator

E_pot_emt = atoms.get_potential_energy()
print(f"Potential energy {E_pot_emt:.5f} eV")

Potential energy -3.07547 eV

The absolute values of the energies are very different between EMT and PFP above, but their absolute values are arbitrary to shift when defining the energy of the system.

Note that the difference between the energies of two different systems consisting of the same number of atoms is meaningful, but the absolute values of each are not so meaningful.

Relation between Atoms and Calculator

In the ASE library, the atomic structure is represented by the Atoms class, and each basic physical property value (energy, force, stress, charge, etc.) for the Atoms can be calculated by setting Calculator to it. Calculator is set directly to atoms.calc.

Relation between Atoms and Calculator

c98abebd195d4ee0a8c8ba9e794a7972

The basic physical properties that can be calculated via Calculator and their calculation methods are as follows.

  • Potential energy: get_potential_energy

  • Force: get_forces

  • Stress: get_stress

  • Charge: get_charges

  • Magnetic moment: get_magnetic_moment

  • Dipole moment: get_dipole_moment

[5]:

calculator = ASECalculator(Estimator())

atoms = bulk("Pt") * (2, 2, 1)
atoms.calc = calculator

E_pot = atoms.get_potential_energy()
charges = atoms.get_charges()
forces = atoms.get_forces()
stress = atoms.get_stress()

print(f"E_pot {E_pot:.2} eV")
print(f"charges {charges} C")
print(f"forces {forces} eV/A")
print(f"stress {stress} eV/A^2")

E_pot -2.2e+01 eV
charges [ 1.56390826e-07  4.79747513e-08  1.49468789e-07 -3.53834395e-07] C
forces [[ 1.21045008e-06 -1.28531060e-07  1.16374716e-06]
 [-8.99028036e-07  9.10295516e-07 -7.65895014e-07]
 [ 8.04090586e-07 -1.76741016e-06  1.52298933e-06]
 [-1.11551263e-06  9.85645703e-07 -1.92084147e-06]] eV/A
stress [-5.98167504e-02 -5.98167416e-02 -5.98167977e-02  9.46340221e-09
 -2.35475377e-08 -5.36077254e-09] eV/A^2

PFP does not support the calculation of magnetic moments or dipole moments. (The magnetic moment is for electronic states and is not used in atomistic simulation, which is the scope of this tutorial.)

[6]:

# Some Calculator supports this, but PFP calculator raises error.

# magmom = atoms.get_magnetic_moment()
# dipole = atoms.get_dipole_moment()

The following [columns] contain advanced content, and can be skipped as you read through the tutorial.

[Column] Calcultor’s cache mechanism

Calculator stores the input atomic structure of the previous calculation in calculator.atoms and the result of the previous calculation in calculator.results. If you call the get_XXX method for the same atomic structure as before, it will skip the re-calculation.

[7]:

atoms.calc = calculator
calculator.reset()

[8]:

%time Epot = atoms.get_potential_energy()
%time Epot = atoms.get_potential_energy()

CPU times: user 3.06 ms, sys: 0 ns, total: 3.06 ms
Wall time: 94 ms
CPU times: user 375 µs, sys: 61 µs, total: 436 µs
Wall time: 383 µs

Comparing the wall time in the above example, we can see that the first run takes time on the order of milliseconds and calculations are performed, but the second run skips the calculations, so the energy is obtained in less than a millisecond.

calculator.reset() can be used to explicitly clear the cache of calculation results and recalculate them. In the following example, you can see that the second calculation still takes milliseconds of calculation time.

[9]:

calculator.reset()

print("----- 1st calc -----")
%time Epot = atoms.get_potential_energy()

print(f"After 1st calc   : {calculator.results}")
calculator.reset()
print(f"After reset      : {calculator.results}")

print("----- 2nd calc -----")
%time Epot = atoms.get_potential_energy()
print(f"After 2nd calc   : {calculator.results}")

----- 1st calc -----
CPU times: user 3.16 ms, sys: 285 µs, total: 3.44 ms
Wall time: 90.3 ms
After 1st calc   : {'energy': -21.712707825955075, 'forces': array([[-1.78798193e-07,  8.57089841e-08, -9.71775560e-07],
       [-3.81674838e-07,  2.40788247e-07,  4.98792271e-07],
       [ 9.00971539e-07, -2.72842101e-07,  1.12511234e-06],
       [-3.40498508e-07, -5.36551305e-08, -6.52129048e-07]]), 'charges': array([-2.17493621e-07, -2.45786026e-07, -9.83995818e-09,  4.73119684e-07]), 'calc_stats': {'elapsed_usec_infer': 36693, 'elapsed_usec': 83022, 'elapsed_usec_preprocess': 3930, 'n_neighbors': 400}, 'free_energy': -21.712707825955075, 'stress': array([-5.98166255e-02, -5.98166500e-02, -5.98166300e-02,  4.26570321e-08,
       -4.45518682e-08,  3.79886465e-08])}
After reset      : {}
----- 2nd calc -----
CPU times: user 1.66 ms, sys: 0 ns, total: 1.66 ms
Wall time: 57.7 ms
After 2nd calc   : {'energy': -21.712705946017945, 'forces': array([[ 9.52835494e-08,  9.18911048e-07, -3.93850900e-07],
       [-7.98279949e-07, -4.07897529e-07,  1.47059723e-06],
       [-3.51072639e-07, -3.76604209e-07,  3.40857921e-07],
       [ 1.05406904e-06, -1.34409311e-07, -1.41760425e-06]]), 'charges': array([-2.72528879e-08, -3.71969477e-08, -1.25965400e-07,  1.90415335e-07]), 'calc_stats': {'n_neighbors': 400, 'elapsed_usec_infer': 46355, 'elapsed_usec': 51052, 'elapsed_usec_preprocess': 3135}, 'free_energy': -21.712705946017945, 'stress': array([-5.98171226e-02, -5.98171468e-02, -5.98171011e-02,  1.97083025e-08,
       -1.43127693e-08,  1.31512884e-08])}

Note that if the atomic structure changes even slightly, the result of calculator.result is discarded and a new calculation is performed.

In the following example, atoms.get_potential_energy is called again after changing the interatomic distance of the hydrogen molecule, and in this case, Calculator detects the change in the input atomic structure and performs a new calculation. You can see that both calculations actually take a wall time on the order of milliseconds.

[10]:

atoms = Atoms(["H", "H"], positions=[[0, 0, 0], [0, 0, 0.8]])

atoms.calc = calculator
%time E_pot1 = atoms.get_potential_energy()
# --- calculator.atoms stores the previously calculated atoms
# print(calculator.atoms.positions)

# Change atomic distance to 2A
atoms.positions[1, 2] = 2.0
# --- The preveously calculated `calculator.atoms` and current calculate target `atoms` are
# different, so calculation is executed.
%time E_pot2 = atoms.get_potential_energy()
# print(calculator.atoms.positions)

print(f"E_pot1 {E_pot1:.2f} eV")
print(f"E_pot2 {E_pot2:.2f} eV")

CPU times: user 3.2 ms, sys: 0 ns, total: 3.2 ms
Wall time: 54 ms
CPU times: user 2.52 ms, sys: 0 ns, total: 2.52 ms
Wall time: 54.1 ms
E_pot1 -4.48 eV
E_pot2 -0.31 eV

[Column] PFP Calculator behavior

※This section is a description of Matlantis-specific behavior.

It describes the behavior of the ASECalculator provided by the pfp-api-client library.

In PFP, potential_energy and charge are calculated by NNP’s forward, and forces and stress are calculated by NNP’s backward.

5df7cb683dda49fb8c0ea6bbed1f50fd

[Column] Potential energy obtained by quantum chemical calculations

In the explanation of total energy at the top of this tutorial, the formula for kinetic energy was given, but there was no formula for potential energy, only explained that it is calculated by Calculator. How is this potential energy expressed theoretically?

If you want to simulate phenomena occurring at the microscopic atomic scale, you will need the theory of quantum mechanics in order to determine the energies, etc., taking into account the electronic states of the system. In the case of the steady state, the Schrödinger equation given below gives the wave function \(\Phi\) and total energy \(E\) for the entire system.

\[\mathcal{H} | \Phi \rangle = E | \Phi \rangle\]

Based on non-relativistic theory, the Hamiltonian \(\mathcal{H}\) for a system consisting of M nuclei and N electrons can be written as follows by atomic unit system.

\[\mathcal{H} = - \sum_{i=1}^N \frac{1}{2} \nabla^2_i - \sum_{A=1}^M \frac{1}{2 M_A} \nabla^2_A - \sum_{i=1}^N \sum_{A=1}^M \frac{Z_A}{r_{iA}} + \sum_{i=1}^N \sum_{j>i}^N \frac{1}{r_{ij}} + \sum_{A=1}^M \sum_{B>A}^M \frac{Z_A Z_B}{R_{AB}}\]

  • \(M_A\): Mass of nucleus \(A\).

  • \(R_{AB}\): Distance between nucleus \(A\) and nucleus \(B\).

  • \(r_{Ai}\): Distance between nucleus \(A\) and electron \(i\).

  • \(r_{ij}\): Distance between electron \(i\) and electron \(j\).

  • \(Z_A\): Atomic number of nucleus \(A\) = number of protons.

and all electrons have charge -1.

The first item corresponds to the kinetic energy of electrons, the second to the kinetic energy of nuclei, the third to the electrostatic interaction energy between electrons and nuclei, the fourth to the electrostatic interaction energy between electrons, and the fifth to the electrostatic interaction energy between nuclei.

Since it is difficult to solve this as is, and a Born-Oppenheimer approximation is often used. Electrons move faster than nuclei because they are much lighter than nuclei. Therefore, the idea is to find the steady state of the electron at that time, assuming that the atom is perfectly stationary.

Then the Schrödinger equation for N electrons is

\[\mathcal{H}_{\rm{elec}} = - \sum_{i=1}^N \frac{1}{2} \nabla^2_i - \sum_{i=1}^N \sum_{A=1}^M \frac{Z_A}{r_{iA}} + \sum_{i=1}^N \sum_{j>i}^N \frac{1}{r_{ij}}\]

and

\[\mathcal{H}_{\rm{elec}} | \Phi_{\rm{elec}} \rangle = E_{\rm{elec}} | \Phi_{\rm{elec}} \rangle \tag{1}\]

The total energy of the nucleus at stationary is written by

\[V = E_{\rm{elec}} + \sum_{A=1}^M \sum_{B>A}^M \frac{Z_A Z_B}{R_{AB}}\]

This corresponds to the potential energy.

If we can solve this, we should be able to obtain the potential energy from the fundamental laws of quantum mechanics alone, without introducing any parameters as in the classical force field. A calculation method based on quantum mechanics that does not depend on experimental values other than the fundamental physical constants is called first-principles calculation (ab initio calculation).

In reality, it is known that the Schrödinger equation (1) for this electron system cannot be solved analytically, and various approximations and assumptions are made by different methods. The Hartree-Fock method and DFT (density functional theory) are typical examples.

So, the potential energy \(V\) handled by ASE is a value that includes the electrostatic potential of electrons and nuclei plus the kinetic energy of electrons, and the kinetic energy \(K\) can be taken as the kinetic energy of nuclei only.

[Reference]

  • “MODERN QUANTUM CHEMISTRY Introduction to Advanced Electronic Structure Theory”, Attila Szabo and Neil S. Ostlund