Tune PME

The tuning is based on the following error formulas:

\[\Delta F_\mathrm{real} \approx \frac{Q^2}{\sqrt{N}} \frac{2}{\sqrt{r_{\text{cutoff}} V}} e^{-r_{\text{cutoff}}^2 / 2 \sigma^2}\]
\[\Delta F_\mathrm{Fourier}^\mathrm{PME} \approx 2\pi^{1/4}\sqrt{\frac{3\sqrt{2} / \sigma}{N(2p+3)}} \frac{Q^2}{L^2}\frac{(\sqrt{2}H/\sigma)^{p+1}}{(p+1)!} \times \exp{\frac{(p+1)[\log{(p+1)} - \log 2 - 1]}{2}} \left< \phi_p^2 \right> ^{1/2}\]

where \(N\) is the number of charges, \(Q^2 = \sum_{i = 1}^N q_i^2\), is the sum of squared charges, \(r_{\text{cutoff}}\) is the short-range cutoff, \(V\) is the volume of the simulation box, \(p\) is the order of the interpolation scheme, \(H\) is the spacing of mesh points, and \(\phi_p^2 = H^{-(p+1)}\prod_{s\in S_H^{(p)}}(x - s)\), in which \(S_H^{(p)}\) is the \(p+1\) mesh points closest to the point \(x\).

torchpme.tuning.tune_pme(charges: Tensor, cell: Tensor, positions: Tensor, cutoff: float, neighbor_indices: Tensor, neighbor_distances: Tensor, full_neighbor_list: bool = False, prefactor: float = 1.0, exponent: int = 1, nodes_lo: int = 3, nodes_hi: int = 7, mesh_lo: int = 2, mesh_hi: int = 7, accuracy: float = 0.001) → tuple[float, dict[str, Any], float][source]

Find the optimal parameters for torchpme.PMECalculator.

For the error formulas are given elsewhere. Note the difference notation between the parameters in the reference and ours:

\[\alpha = \left(\sqrt{2}\,\mathrm{smearing} \right)^{-1}\]
Parameters:

  • charges (Tensor) – torch.tensor of shape (len(positions, 1)) containing the atomic (pseudo-)charges

  • cell (Tensor) – torch.tensor of shape (3, 3), where cell[i] is the i-th basis vector of the unit cell

  • positions (Tensor) – torch.tensor of shape (N, 3) containing the Cartesian coordinates of the N particles within the supercell.

  • cutoff (float) – float, cutoff distance for the neighborlist

  • neighbor_indices (Tensor) – torch.tensor with the i,j indices of neighbors for which the potential should be computed in real space.

  • neighbor_distances (Tensor) – torch.tensor with the pair distances of the neighbors for which the potential should be computed in real space.

  • full_neighbor_list (bool) – If set to True, a “full” neighbor list is expected as input. This means that each atom pair appears twice. If set to False, a “half” neighbor list is expected.

  • prefactor (float) – electrostatics prefactor; see Prefactors for details and common values.

  • exponent (int) – \(p\) in \(1/r^p\) potentials, currently only \(p=1\) is supported

  • nodes_lo (int) – Minimum number of interpolation nodes

  • nodes_hi (int) – Maximum number of interpolation nodes

  • mesh_lo (int) – Controls the minimum number of mesh points along the shortest axis, \(2^{mesh_lo}\)

  • mesh_hi (int) – Controls the maximum number of mesh points along the shortest axis, \(2^{mesh_hi}\)

  • accuracy (float) – Recomended values for a balance between the accuracy and speed is \(10^{-3}\). For more accurate results, use \(10^{-6}\).

Returns:

Tuple containing a float of the optimal smearing for the :class: CoulombPotential, a dictionary with the parameters for PMECalculator and a float of the optimal cutoff value for the neighborlist computation, and the timing of this set of parameters.

Return type:

tuple[float, dict[str, Any], float]

Example

>>> import torch

To allow reproducibility, we set the seed to a fixed value

>>> _ = torch.manual_seed(0)
>>> positions = torch.tensor([[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]])
>>> charges = torch.tensor([[1.0], [-1.0]])
>>> cell = torch.eye(3)
>>> neighbor_distances = torch.tensor(
...     [0.9381, 0.9381, 0.8246, 0.9381, 0.8246, 0.8246, 0.6928],
... )
>>> neighbor_indices = torch.tensor(
...     [[0, 1], [0, 1], [0, 1], [0, 1], [0, 1], [0, 1], [0, 1]]
... )
>>> smearing, parameter, timing = tune_pme(
...     charges,
...     cell,
...     positions,
...     cutoff=1.0,
...     neighbor_distances=neighbor_distances,
...     neighbor_indices=neighbor_indices,
...     accuracy=1e-1,
... )
class torchpme.tuning.pme.PMEErrorBounds(charges: Tensor, cell: Tensor, positions: Tensor)[source]

Error bounds for torchpme.PMECalculator.

Note

The torchpme.tuning.pme.PMEErrorBounds.forward() method takes floats as the input, in order to be in consistency with the rest of the package – these parameters are always float but not torch.Tensor. This design, however, prevents the utilization of torch.autograd and other torch features. To take advantage of these features, one can use the torchpme.tuning.pme.PMEErrorBounds.err_rspace() and torchpme.tuning.pme.PMEErrorBounds.err_kspace(), which takes torch.Tensor as parameters.

Parameters:

  • charges (Tensor) – torch.tensor of shape (len(positions, 1)) containing the atomic (pseudo-)charges

  • cell (Tensor) – torch.tensor of shape (3, 3), where cell[i] is the i-th basis vector of the unit cell

  • positions (Tensor) – torch.tensor of shape (N, 3) containing the Cartesian coordinates of the N particles within the supercell.

Example

>>> import torch
>>> positions = torch.tensor(
...     [[0.0, 0.0, 0.0], [0.4, 0.4, 0.4]], dtype=torch.float64
... )
>>> charges = torch.tensor([[1.0], [-1.0]], dtype=torch.float64)
>>> cell = torch.eye(3, dtype=torch.float64)
>>> error_bounds = PMEErrorBounds(charges, cell, positions)
>>> print(
...     error_bounds(
...         smearing=1.0, mesh_spacing=0.5, cutoff=4.4, interpolation_nodes=3
...     )
... )
tensor(0.0011, dtype=torch.float64)
err_kspace(smearing: Tensor, mesh_spacing: Tensor, interpolation_nodes: Tensor) → Tensor[source]

The Fourier space error of PME.

Parameters:
Return type:

Tensor

err_rspace(smearing: Tensor, cutoff: Tensor) → Tensor[source]

The real space error of PME.

Parameters:
Return type:

Tensor

error(cutoff: float, smearing: float, mesh_spacing: float, interpolation_nodes: float) → Tensor[source]

Calculate the error bound of PME.

\[\text{Error}_{\text{total}} = \sqrt{\text{Error}_{\text{real space}}^2 + \text{Error}_{\text{Fourier space}}^2\]
Parameters:

  • smearing (float) – if its value is given, it will not be tuned, see torchpme.PMECalculator for details

  • mesh_spacing (float) – if its value is given, it will not be tuned, see torchpme.PMECalculator for details

  • cutoff (float) – if its value is given, it will not be tuned, see torchpme.PMECalculator for details

  • interpolation_nodes (float) – The number n of nodes used in the interpolation per coordinate axis. The total number of interpolation nodes in 3D will be n^3. In general, for n nodes, the interpolation will be performed by piecewise polynomials of degree n - 1 (e.g. n = 4 for cubic interpolation). Only the values 3, 4, 5, 6, 7 are supported.

Return type:

Tensor