Ctrl+K

System#

class system.System#

Class that represents a system with periodic boundary conditions to facilitate orbital-free density functional theory calculations using Pytorch functions. This enables the use Pytorch’s autograd to compute gradient dependent terms via auto-differentiations. The calculations are also GPU compatible due to the use of Pytorch functions.

__init__(box_vecs, shape, ions, terms, units='b', coord_type='cartesian', Rc=None, pme_order=None, device=device(type='cpu'))#
Parameters:

  • box_vecs (torch.Tensor) –

    Lattice vectors \(\mathbf{a},~\mathbf{b},~\mathbf{c}\) with input format

    [[\(a_x\), \(a_y\), \(a_z\)], [\(b_x\), \(b_y\), \(b_z\)], [\(c_x\), \(c_y\), \(c_z\)]]

  • shape (torch.Size or iterable) – Real-space grid shape

  • ions (list) –

    Ion information, with each sublist containing the following

    [name (string), path_to_pseudopotential_file (string), ionic_coordinates (torch.Tensor)]

    For example, ['Al', 'pseudopots/Al.recpot', torch.tensor([[0,0,0]])]

  • units (string) – Units of a for angstrom or b for bohr.

  • Rc (None or float) – Cutoff radius for ion-ion interaction summation (in bohr). Default behaviour (Rc=None) uses the heuristic \(R_c = 12 h_\text{max}\) where \(h_\text{max}\) is the largest interplanar spacing for the cell.

  • coord_type (string) – Whether cartesian or fractional coordinates were used to represent the ionic coordinates in “ions”

  • pme_order (None or even int) – The spline order of the particle-mesh Ewald scheme for a quasi-linear scaling computation of the structure factor. The default value of ‘None’ indicates that the exact quadratic scaling method is used for computing the structure factor.

  • device (torch.device) – Device to store System tensors in. Default behaviour is to use CPU.

classmethod ecut2shape(energy_cutoff, box_vecs)#

Computes the shape of grid for a lattice given an energy cutoff.

Parameters:

  • ecut (float) – Energy cutoff given in eVs

  • box_vecs (torch.Tensor) – Lattice vectors given in Angstroms

Returns:

Real space grid shape

Return type:

tuple

set_device(device=device(type='cpu'))#

Moves all System tensors to the specified device. By default, the device is set to a GPU if is avaliable, otherwise the device is set as a CPU.

Parameters:

device (torch.device) – Device

place_ions(ion_coords, coord_type='cartesian', units='a', initialization=False)#

Places ions at given positions.

Parameters:

  • ion_coords (torch.Tensor) – Ionic coordinates

  • coord_type (string) – Whether cartesian or fractional coordinates were used to represent the ionic coordinates

  • units (string) – Units of a (Angstrom) or b (Bohr)

  • initialization (bool) – Whether this function is called during the initialization of a new system object or not

set_lattice(box_vecs, units='a', initialization=False)#

Sets lattice vectors.

Parameters:

  • box_vecs (torch.Tensor) – Lattice vectors

  • units (string) – Units of a (Angstrom) or b (Bohr)

  • initialization (bool) – Whether this function is called during the initialization of a new system object or not

set_potential(pot)#

Sets the external potential to a given one.

Parameters:

pot (torch.Tensor) – Given external potential (must match the system’s shape attribute)

initialize_density()#

Initializes the density to that of a uniform density profile.

set_density(den)#

Sets the electron density to a given one.

Parameters:

den (torch.Tensor) – Given electron density (must match the system’s shape attribute)

set_electron_number(N)#

Sets the number of electrons in the system.

Parameters:

N (float) – Number of electrons

detach()#

Detaches all core variables from any computational graphs they may be in. This is used primarily as a clean-up tool to avoid autograd-related bugs.

device()#

Returns the device that the system’s parameters are stored in.

Returns:

Device

Return type:

torch.device

ion_count()#

Returns the number of ions in the system.

Returns:

Number of ions

Return type:

int

electron_count()#

Returns the number of electrons in the system.

Returns:

Number of electrons

Return type:

int

lattice_vectors(units='a')#

Returns the lattice vectors of the system.

Parameters:

units (string) – Units of a (Angstrom) or b (Bohr)

Returns:

Lattice vectors

Return type:

torch.Tensor

ions()#

Returns list of ions in the system.

Returns:

List of ions

Return type:

list

cartesian_ionic_coordinates(units='a')#

Returns the Cartesian ionic coordinates of the system.

Parameters:

units (string) – Units of a (Angstrom) or b (Bohr)

Returns:

Cartesian ionic coordinates

Return type:

torch.Tensor

fractional_ionic_coordinates()#

Returns the fractional ionic coordinates of the system.

Returns:

Fractional ionic coordinates

Return type:

torch.Tensor

ionic_potential(units='Ha')#

Returns the ionic potential / external potential of the system.

Parameters:

units (string) – Units of energy (‘Ha’ or ‘eV’)

Returns:

Ionic potential

Return type:

torch.Tensor

density(requires_grad=False)#

Returns the electron density of the system.

Parameters:

requires_grad (bool) – Whether the returned density can be differentiated (used for training kinetic functionals for example)

Returns:

Electron density in \(\text{bohr}^{-3}\)

Return type:

torch.Tensor

check_density_convergence(method='dEdchi')#

Utility function that computes the measures of density convergence in atomic units. I.e. how well the Euler equation,

\[\frac{\delta E}{\delta n(\mathbf{r})} = \mu,\]

is obeyed. \(\mu\) is the chemical potential.

If the method = 'dEdchi', this function computes the quantity

\[\text{Max}\left(\left|\frac{\delta E}{\delta \chi(\mathbf{r})}\right|\right) ~~\text{where} ~~ n(\mathbf{r}) = \frac{N_e \chi^2(\mathbf{r})}{\int~d^3\mathbf{r}'\chi^2(\mathbf{r}')}\]

If the method = 'euler', this function computes the quantity

\[\text{Max}\left(\left|\mu - \frac{\delta E}{\delta n(\mathbf{r})}\right|\right) ~~\text{where} ~~ \mu = \frac{1}{N_e} \int~d^3\mathbf{r} \frac{\delta E}{\delta n(\mathbf{r})} n(\mathbf{r})\]

\(N_e\) is the number of electrons.

Parameters:

method (string) – dEdchi or euler

Returns:

A measure of density convergence

Return type:

float

functional_derivative(type='density', requires_grad=False)#

Returns, in atomic units, the functional derivative of the system,

\[\frac{\delta E}{\delta n(\mathbf{r})} ~~\text{or} ~~ \frac{\delta E}{\delta \chi(\mathbf{r})} ~~\text{where} ~~ n(\mathbf{r}) = \frac{N_e \chi^2(\mathbf{r})}{\int~d^3\mathbf{r}'\chi^2(\mathbf{r}')}.\]
Parameters:

  • type (str) – Functional derivative computed with respect to density or chi

  • requires_grad (bool) – Whether the returned energy has requires_grad = True or not (used for training kinetic functionals for example)

Returns:

Functional derivative

Return type:

torch.Tensor

chemical_potential()#

Returns the chemical potential of the system.

Returns:

Chemical potential

Return type:

torch.Tensor

energy(units='Ha', requires_grad=False)#

Returns the energy of the system.

Parameters:

  • units (string) – Units of energy (Ha or eV)

  • requires_grad (bool) – Whether the returned energy has requires_grad = True or not (used for training kinetic functionals for example)

Returns:

Energy

Return type:

float or torch.Tensor (depending on requires_grad)

volume(units='b3')#

Computes the cell volume \(\Omega\).

Parameters:

units (string) – Units of the pressure returned (b3 or a3)

Returns:

Volume

Return type:

float

pressure(units='Ha/b3', requires_grad=False)#

Computes, via autograd and/or Xitorch, the pressure

\[P = \frac{dE[n]}{d\Omega}\]

where the cell volume is \(\Omega\).

Parameters:

  • units (string) – Units of the pressure returned (Ha/b3, eV/a3 or GPa)

  • requires_grad (bool) – Whether the returned pressure has requires_grad = True or not (used for training kinetic functionals for example)

Returns:

Pressure

Return type:

float or torch.Tensor (depending on requires_grad)

enthalpy(units='Ha')#

Returns the enthalpy, \(H=U+PV\), of the system.

Parameters:

units (string) – Units of energy (Ha or eV)

Returns:

Enthalpy

Return type:

float

bulk_modulus(units='Ha/b3', requires_grad=False)#

Computes, via autograd and Xitorch, the bulk modulus

\[K = \Omega \frac{d^2 E[n]}{d \Omega^2}\]

where the cell volume is \(\Omega\).

Parameters:

  • units (string) – Units of the bulk modulus returned (Ha/b3, eV/a3 or GPa)

  • requires_grad (bool) – Whether the returned bulk modulus has requires_grad = True or not (used for training kinetic functionals for example)

Returns:

Bulk modulus

Return type:

float or torch.Tensor (depending on requires_grad)

eos_fit(f=0.05, N=9, eos='bm', verbose=False, plot=False, **den_opt_kwargs)#

Perform a Birch-Murnaghan equation of state fit for a given system. The volume of the system is taken to be an approximation of the true equilibrium volume.

The parameters are returned in the following order

  1. Equilibrium bulk modulus, \(K_0\) [GPa]

  2. Equilibrium bulk modulus derivative wrt pressure, \(K'_0\) [dimensionless]

  3. Equilibrium energy, \(E_0\) [eV]

  4. Equilibrium volume, \(V_0\) [ų]

Parameters:

  • f (float) – Fraction by which the volume is compressed and stretched by (default: 0.05)

  • N (int) – Number of energy-volume data points used for the fit (default: 9)

  • eos (string) – bm (default) for Birch-Murnaghan or m for Murnaghan

  • verbose (boolean) – Whether the energy-volume values are printed

  • plot (boolean) – Whether an energy-volume plot is made

  • den_opt_kwargs – Arguments for density optimization. The default values are: {'ntol': 1e-10, 'n_conv_cond_count': 3, 'n_method': 'LBFGS', 'n_step_size': 0.1, 'n_maxiter': 1000, 'conv_target': 'dE', 'n_verbose': False, 'from_uniform': False}

Returns:

Fitted parameters, Fitting errors

Return type:

ndarray

forces(units='Ha/b')#

Computes, via autograd, the forces

\[F_{\alpha, i} = - \frac{dE[n]}{dR_{\alpha, i}}\]

where \(\alpha\) is an ion index, \(i \in \{x,y,z\}\) and \(R_{\alpha, i}\) are the Cartesian ionic coordinates.

Parameters:

units (string) – Units of the forces returned (Ha/b or eV/a)

Returns:

Forces

Return type:

torch.Tensor

stress(units='Ha/b3')#

Computes, via autograd, the stress tensor

\[\sigma_{ij} = \frac{1}{\Omega} \sum_k \frac{\partial E[n]}{\partial h_{ik}} h_{jk}\]

where \(h_{ij}\) are matrix elements of a matrix whose columns are lattice vectors and the cell volume is \(\Omega\).

Parameters:

units (string) – Units of the stress tensor returned (Ha/b3, eV/a3 or GPa)

Returns:

Stress tensor

Return type:

torch.Tensor

elastic_constants(units='Ha/b3')#

Computes, via autograd and Xitorch, the elastic constants (Birch coefficients)

\[C_{ijk\ell} = \frac{\partial \sigma_{ij}}{\partial \epsilon_{k\ell}} = \sum_m \frac{\partial \sigma_{ij}}{\partial h_{km}} h_{\ell m}\]

where \(h_{ij}\) are matrix elements of a matrix whose columns are lattice vectors.

Parameters:

units (string) – Units of the elastic constants returned (Ha/b3, eV/a3 or GPa)

Returns:

Elastic constants

Return type:

torch.Tensor

force_constants(primitive_ion_indices, units='eV/a2')#

Computes, via autograd, the force constants

\[\Phi_{\alpha, i, \beta, j} = - \frac{dF_{\alpha, i} }{dR_{\beta, j}}\]

where \(\alpha, \beta\) are ion indices, \(i,j \in \{x,y,z\}\) and \(R_{\beta, j}\) are the Cartesian ionic coordinates.

Parameters:

  • units (string) – Units of the force constants returned (Ha/b2 or eV/a2)

  • primitive_ion_indices (list) – List of indices corresponding to ions in the primitive cell. primitive_ion_indices begins counting from 0

Returns:

Force constants with shape [len(primitive_ion_indices), N_ions, 3, 3]

Return type:

torch.Tensor

set_Rc(Rc=None)#

Sets the cutoff radius for the ion-ion interaction energy. Default behaviour (Rc=None) uses the heuristic \(R_c = 12 h_\text{max}\) where \(h_\text{max}\) is the largest interplanar spacing for the cell.

Parameters:

Rc (None or float) – Cutoff radius for ion-ion interaction summation (in bohr).

optimize_density(ntol=1e-07, n_conv_cond_count=3, n_method='LBFGS', n_step_size=0.1, n_maxiter=1000, conv_target='dE', n_verbose=False, from_uniform=False, potentials=None)#

Performs an orbital-free density optimization procedure to minimize the energy via a direct optimization using a Pytorch-based modified LBFGS [10.1109/eScience.2018.00112] optimizer, or a Two-Point Gradient Descent (TPGD) [IMA Journal of Numerical Analysis, Volume 8, Issue 1, January 1988, Pages 141–148] algorithm.

Parameters:

  • ntol (float) – Convergence tolerance for density optimization. The optimization procedure stops when the convergence target (conv_target) variable is lower than ntol for a number (n_conv_cond_count) of consecutive iterations (default value is 1e-7)

  • n_conv_cond_count (int) – ‘Convergence condition count’, which is the number of times the convergence condition has to be met in consecutive iterations before the optimization terminates

  • n_method (string) – LBFGS or TPGD

  • n_step_size (float) – Step size for the optimizers

  • n_maxiter (int) – Maximum number of density optimization iterations

  • conv_target (string) – dE (energy difference), dEdchi (Max \(|\delta E/\delta \chi|\)), or euler (Max \(|\mu-\delta E/\delta n|\))

  • n_verbose (bool) – Whether the density optimization progress is printed out or not

  • from_uniform (bool) – Whether the density optimization begins from a uniform density or not. When ‘False’, the behaviour is to compare the energy computed with the current density and the energy computed with a uniform density and initialize the density to whichever that yields a lower energy.

  • potentials (function) – Explicitly implemented functional derivatives (can be used to check that analytically derived and explicitly implemented functional derivatives are correct)

optimize_geometry(ftol=0.02, stol=0.002, g_conv_cond_count=3, g_method='LBFGSlinesearch', g_step_size=0.1, g_maxiter=1000, g_verbose=False, **den_opt_kwargs)#

Performs a geometry optimization to minimize the energy by varying the ionic positions and/or lattice vectors, minimizing the ionic forces and stress in the process.

Set ftol = None to only vary lattice vectors. Set stol = None to only vary ionic positions.

Parameters:

  • ftol (float) – Force tolerance in eV/Å - optimization terminates when the largest force component goes below this value for a number (g_conv_cond_count) of consecutive iterations

  • stol (float) – Stress tolerance in eV/ų - optimization terminates when the largest stress component goes below this value for a number (g_conv_cond_count) of consecutive iterations

  • g_conv_cond_count (int) – ‘Convergence condition count’, which is the number of times the convergence condition has to be met in consecutive iterations before the optimization terminates

  • g_method (string) – LBFGSlinesearch, LBFGS, RPROP or TPGD

  • g_step_size (float) – Step size for the optimizer

  • g_maxiter (int) – Maximum number of geometry optimization iterations

  • g_verbose (bool) – Whether the geometry optimization progress is printed out or not

  • den_opt_kwargs – Arguments for density optimization. The default values are: {'ntol': 1e-10, 'n_conv_cond_count': 3, 'n_method': 'LBFGS', 'n_step_size': 0.1, 'n_maxiter': 1000, 'conv_target': 'dE', 'n_verbose': False, 'from_uniform': False}

Returns:

Whether the optimization was successful or not

Return type:

bool

optimize_parameterized_geometry(params, parameterized_geometry, ftol=0.02, stol=0.002, g_conv_cond_count=3, g_method='LBFGSlinesearch', g_step_size=0.1, g_maxiter=1000, g_verbose=False, param_string=None, **den_opt_kwargs)#

Performs a parameterized geometry optimization to minimize the energy by varying the parameters that describe the ionic positions and/or lattice vectors.

Set ftol = None to only vary lattice vectors. Set stol = None to only vary ionic positions.

Parameters:

  • params (torch.Tensor) – The parameters that describe the geometry

  • parameterized_geometry (function) – Function that accepts params as input and returns box_vecs, frac_ion_coords, that is the lattice vectors and fractional ionic coordinates based on the parameters (box_vecs must be in units of Bohr)

  • ftol (float) – Force tolerance in eV/Å - optimization terminates when the largest force component goes below this value for a number (g_conv_cond_count) of consecutive iterations

  • stol (float) – Stress tolerance in eV/ų - optimization terminates when the largest stress component goes below this value for a number (g_conv_cond_count) of consecutive iterations

  • g_conv_cond_count (int) – ‘Convergence condition count’, which is the number of times the convergence condition has to be met in consecutive iterations before the optimization terminates

  • g_method (string) – LBFGSlinesearch, LBFGS, RPROP or TPGD

  • g_step_size (float) – Step size for the optimizer

  • g_maxiter (int) – Maximum number of geometry optimization iterations

  • g_verbose (bool) – Whether the geometry optimization progress is printed out or not

  • param_string (function) – For printing out the parameter values during the optimization

  • den_opt_kwargs – Arguments for density optimization. The default values are: {'ntol': 1e-10, 'n_conv_cond_count': 3, 'n_method': 'LBFGS', 'n_step_size': 0.1, 'n_maxiter': 1000, 'conv_target': 'dE', 'n_verbose': False, 'from_uniform': False}

Returns:

Whether the optimization was successful or not

Return type:

bool

On this page
Show Source