Ewald sum matrix

Ewald sum matrix [1] can be viewed as a logical extension of the Coulomb matrix for periodic systems, as it models the interaction between atoms in a periodic crystal through electrostatic interaction. However, in order to calculate the Coulomb interaction between the atomic cores in a periodic crystal, the Ewald summation technique has to be used together with a uniform neutralizing background charge.

Notice that the terms of the Ewald sum matrix in DScribe differ slightly from the original article. We correct an issue in the energy term related to the self-energy and the charged-cell correction. In the original article [1] (equation 20) this correction for the off-diagonal elements was defined as:

\[\phi^{\text{self}}_{ij} + \phi^{\text{bg}}_{ij} = -\frac{\alpha}{\sqrt{\pi}}(Z_i^2 + Z_j^2) -\frac{\pi}{2V\alpha^2}(Z_i + Z_j)^2~\forall~i \neq j\]

This term does however not correspond to the correct Ewald energy, as seen in the case of two interacting atoms, one of which has no charge: \(Z_i = 0\), \(Z_j \neq 0\). In this case the interaction should be zero, but the given term is non-zero and additionally depends on the screening parameter \(\alpha\). DScribe instead uses the corrected term:

\[\phi^{\text{self}}_{ij} + \phi^{\text{bg}}_{ij} = -\frac{\pi}{2 V \alpha^2} Z_i Z_j~\forall~i \neq j\]

Setup

Instantiating the object that is used to create Ewald sum matrices can be done as follows:

r_cut = 6.0
nmax = 8
lmax = 6

# Setting up the Ewald sum matrix descriptor
esm = EwaldSumMatrix(
    n_atoms_max=6,
)

# Creation
from ase.build import bulk

The constructor takes the following parameters:

Creation

After the Ewald sum matrix has been set up, it may be used on periodic atomic structures with the create()-method.

from ase.build import bulk

# NaCl crystal created as an ASE.Atoms
nacl = bulk("NaCl", "rocksalt", a=5.64)

# Create output for the system
nacl_ewald = esm.create(nacl)

# Create output for multiple system
al = bulk("Al", "fcc", a=4.046)
fe = bulk("Fe", "bcc", a=2.856)
samples = [nacl, al, fe]
ewald_matrices = esm.create(samples)            # Serial
ewald_matrices = esm.create(samples, n_jobs=2)  # Parallel

The call syntax for the create-function is as follows:

Note that if you specify in n_atoms_max a lower number than atoms in your structure it will cause an error. The output will in this case be a flattened matrix, specifically a numpy array with size #atoms * #atoms. The number of features may be requested beforehand with the get_number_of_features()-method.

In the case of multiple samples, the creation can also be parallellized by using the n_jobs-parameter. This splits the list of structures into equally sized parts and spaws a separate process to handle each part.

Examples

The following examples demonstrate usage of the descriptor. These examples are also available in dscribe/examples/ewaldsummatrix.py.

Accuracy

Easiest way to control the accuracy of the Ewald summation is to use the accuracy-parameter. Lower values of this parameter correspond to tighter convergence criteria and better accuracy.

ewald_1 = esm.create(nacl, accuracy=1e-3)
ewald_2 = esm.create(nacl, accuracy=1e-5)

Another option is to directly provide the real- and reciprocal space cutoffs:

ewald_3 = esm.create(nacl, r_cut=10, g_cut=10)

Total electrostatic energy

Let’s calculate the electrostatic energy of a crystal by using the information contained in the Ewald sum matrix.

# Ewald summation parameters
r_cut = 40
g_cut = 40
a = 3

# Calculate Ewald sum matrix with DScribe
ems = EwaldSumMatrix(
    n_atoms_max=3,
    permutation="none",
    flatten=False
)
ems_out = ems.create(al, a=a, r_cut=r_cut, g_cut=g_cut)

# Calculate the total electrostatic energy of the crystal
total_energy = ems_out[0, 0] + ems_out[1, 1] + ems_out[0, 1]

# Converts unit of q*q/r into eV
conversion = 1e10 * scipy.constants.e / (4 * math.pi * scipy.constants.epsilon_0)
total_energy *= conversion
print(total_energy)