mlipops.ChargeEquilibration#

class mlipops.ChargeEquilibration(coulomb: CoulombNC | CoulombRF | CoulombEwald)#

Compute atomic partial charges with charge equilibration.

This class implements the charge equilibration method as described in http://dx.doi.org/10.1103/PhysRevB.92.045131. It models atoms as Gaussian charge distributions. You provide three parameters for each atom: electronegativity, which describes its innate affinity for electrons; hardness, which describes its resistance to changes in charge; and radius, which is the width of the Gaussian charge distribution. It solves for the charges that minimize an energy function subject to constraints on the total charge.

In the simplest version, you provide a single number for the total charge of the system, and it constrains the atomic partial charges to add up to the correct value. This works well for single isolated molecules, but cannot accurately describe systems of multiple molecules. In particular, it unrealistically predicts fractional charges for molecules. Instead you can provide a list of the atoms that make up each molecule and the total charge of each one. It then solves the equations subject to a separate constraint for each molecule. This produces more realistic results, but requires that you know in advance how the atoms are divided into molecules and how the charge is divided among them.

You can optionally provide an external electric potential that should be used to polarize the atoms. This might be computed from a uniform external field, or by calling compute_potential() on a Coulomb calculation object to get the potential resulting from a set of external charges.

You can optionally specify a default charge for each atom. The energy function is then modified as described in https://doi.org/10.1021/jz3008485 to bias each atom towards its default charge. In some cases, using nonzero default charges (for example, the formal oxidation state of each atom) can lead to more accurate results.

This class offers a choice of method for solving the system of equations. By default it uses torch.linalg.solve(), which implements a direct algorithm. It is accurate and generally robust, but it can be slow, especially for large systems. Alternatively you can choose MINRES, an efficient iterative algorithm. It can be much faster in some cases, but this comes at the cost of somewhat lower accuracy.

__init__(coulomb: CoulombNC | CoulombRF | CoulombEwald)#

Create an object for performing charge equilibration.

Parameters:

coulomb (CoulombNC | CoulombRF | CoulombEwald) – the object used to compute Coulomb interactions between atoms. This determines what method is used for handling long range interactions and the neighbor list used for identifying interacting pairs.

forward(positions: Tensor, electronegativity: Tensor, hardness: Tensor, radius: Tensor, total_charge: float | Tensor | None = None, molecules: list | None = None, box_vectors: Tensor | None = None, potential: Tensor | None = None, default_charge: Tensor | None = None, batch: Tensor | None = None, solver: str = 'direct') → Tensor#

Perform charge equilibration to compute atomic partial charges.

Parameters:

  • positions (torch.Tensor) – a Tensor of shape (n_particles, 3) containing the Cartesian coordinates of each particle

  • electronegativity (torch.Tensor) – a Tensor of shape (n_particles,) containing the electronegativity ($chi$) of each particle

  • hardness (torch.Tensor,) – a Tensor of shape (n_particles,) containing the hardness ($J_{ii}$) of each particle

  • radius (torch.Tensor) – a Tensor of shape (n_particles,) containing the radius ($alpha$) of each particle

  • total_charge (float | torch.Tensor | None) – the total charge of the system. If batch is None, this should be a float or scalar Tensor. If batch is not None, this should be a Tensor of shape (n_systems,) containing the charge of each system. You must specify either total_charge or molecules, but not both.

  • molecules (list | None) – the list of molecules. Each element should be a tuple with two elements. The first element is a Tensor containing the indices of the particles that belong to the molecule. The second element is a float with the total charge of the molecule. You must specify either total_charge or molecules, but not both.

  • box_vectors (torch.Tensor | None) – if batch is None, a Tensor of shape (3, 3) containing box vectors defining the periodic box. If batch is not None, a Tensor of shape (n_systems, 3, 3) containing the box vectors for each system. If None, periodic boundary conditions are not used.

  • potential (torch.Tensor) – a Tensor of shape (n_particles,) containing the external electric potential at the location of each particle

  • default_charge (torch.Tensor | None) – a Tensor of shape (n_particles,) containing the default charge of every particle. If None, all default charges are 0.

  • solver (str) – the method to use for solving the system of equations. Options are ‘direct’ (use torch.linalg.solve() to directly compute the result) and ‘minres’ (an iterative solver that tends to be faster, especially for large systems, at the cost of slightly lower accuracy).

  • batch (torch.Tensor | None) – a Tensor of shape (n_particles,) containing the index of the system each particle belongs to. This must be sorted in ascending order, and every system must contain at least one particle. If None, the calculation is performed for a single system instead of a batch of systems.

Returns:

a Tensor of shape (n_particles,) containing the charge of each particle

Return type:

torch.Tensor