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:

from dscribe.descriptors import EwaldSumMatrix

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

The constructor takes the following parameters:

EwaldSumMatrix.__init__(n_atoms_max, permutation='sorted_l2', sigma=None, seed=None, sparse=False, dtype='float64')
Parameters:

  • n_atoms_max (int) – The maximum nuber of atoms that any of the samples can have. This controls how much zeros need to be padded to the final result.

  • permutation (string) –

    Defines the method for handling permutational invariance. Can be one of the following:

    • none: The matrix is returned in the order defined by the Atoms.

    • sorted_l2: The rows and columns are sorted by the L2 norm.

    • eigenspectrum: Only the eigenvalues are returned sorted by their absolute value in descending order.

    • random: The rows and columns are sorted by their L2 norm after applying Gaussian noise to the norms. The standard deviation of the noise is determined by the sigma-parameter.

  • sigma (float) – Provide only when using the random-permutation option. Standard deviation of the gaussian distributed noise determining how much the rows and columns of the randomly sorted matrix are scrambled.

  • seed (int) – Provide only when using the random-permutation option. A seed to use for drawing samples from a normal distribution.

  • sparse (bool) – Whether the output should be a sparse matrix or a dense numpy array.

Creation

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

# Creation
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:

EwaldSumMatrix.create(system, accuracy=1e-05, w=1, r_cut=None, g_cut=None, a=None, n_jobs=1, only_physical_cores=False, verbose=False)[source]

Return the Ewald sum matrix for the given systems.

Parameters:

  • system (ase.Atoms or list of ase.Atoms) – One or many atomic structures.

  • accuracy (float) – The accuracy to which the sum is converged to. Corresponds to the variable \(A\) in https://doi.org/10.1080/08927022.2013.840898. Used only if g_cut, r_cut and a have not been specified. Provide either one value or a list of values for each system.

  • w (float) – Weight parameter that represents the relative computational expense of calculating a term in real and reciprocal space. This has little effect on the total energy, but may influence speed of computation in large systems. Note that this parameter is used only when the cutoffs and a are set to None. Provide either one value or a list of values for each system.

  • r_cut (float) – Real space cutoff radius dictating how many terms are used in the real space sum. Provide either one value or a list of values for each system.

  • g_cut (float) – Reciprocal space cutoff radius. Provide either one value or a list of values for each system.

  • a (float) – The screening parameter that controls the width of the Gaussians. If not provided, a default value of \(\alpha = \sqrt{\pi}\left(\frac{N}{V^2}\right)^{1/6}\) is used. Corresponds to the standard deviation of the Gaussians. Provide either one value or a list of values for each system.

  • n_jobs (int) – Number of parallel jobs to instantiate. Parallellizes the calculation across samples. Defaults to serial calculation with n_jobs=1. If a negative number is given, the used cpus will be calculated with, n_cpus + n_jobs, where n_cpus is the amount of CPUs as reported by the OS. With only_physical_cores you can control which types of CPUs are counted in n_cpus.

  • only_physical_cores (bool) – If a negative n_jobs is given, determines which types of CPUs are used in calculating the number of jobs. If set to False (default), also virtual CPUs are counted. If set to True, only physical CPUs are counted.

  • verbose (bool) – Controls whether to print the progress of each job into to the console.

Returns:

Ewald sum matrix for the given systems. The return type depends on the ‘sparse’-attribute. The first dimension is determined by the amount of systems.

Return type:

np.ndarray | sparse.COO

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 n_atoms * n_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")
ems_out = ems.create(al, a=a, r_cut=r_cut, g_cut=g_cut)
ems_out = ems.unflatten(ems_out)

# 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)