Replace MeshPotential by PMEPotential (#15)

Also added calculator for Ewald summation `EwaldPotential`  and direct interaction without a cell `DirectPotential`

---------

Co-authored-by: Philip Loche <ploche@physik.fu-berlin.de>

.readthedocs.yaml

@@ -25,4 +25,5 @@ python:
     - method: pip
       path: .
       extra_requirements:
+        - examples
         - metatensor

README.rst

@@ -22,9 +22,8 @@ You can install *MeshLode* using pip with
 
 You can then ``import meshlode`` and use it in your projects!
 
-We also provide bindings to `metatensor
-<https://lab-cosmo.github.io/metatensor/latest/>`_ which can optionally be installed
-together and used as ``meshlode.metatensor`` via
+We also provide bindings to `metatensor <https://docs.metatensor.org/latest/>`_ which
+can optionally be installed together and used as ``meshlode.metatensor`` via
 
 .. code-block:: bash
 

docs/src/conf.py

@@ -55,10 +55,11 @@
 autodoc_typehints_format = "short"
 
 intersphinx_mapping = {
+    "ase": ("https://wiki.fysik.dtu.dk/ase/", None),
     "python": ("https://docs.python.org/3", None),
     "numpy": ("https://numpy.org/doc/stable/", None),
     "torch": ("https://pytorch.org/docs/stable/", None),
-    "metatensor": ("https://lab-cosmo.github.io/metatensor/latest/", None),
+    "metatensor": ("https://docs.metatensor.org/latest/", None),
 }
 
 # -- Options for HTML output -------------------------------------------------

docs/src/references/calculators/directpotential.rst

@@ -0,0 +1,6 @@
+DirectPotential
+###############
+
+.. autoclass:: meshlode.DirectPotential
+    :members:
+    :undoc-members:

docs/src/references/calculators/ewaldpotential.rst

@@ -0,0 +1,6 @@
+EwaldPotential
+##############
+
+.. autoclass:: meshlode.EwaldPotential
+    :members:
+    :undoc-members:

docs/src/references/calculators/index.rst

@@ -14,9 +14,8 @@ Calculators return the representations as a :py:obj:`List` of :py:class:`torch.T
 We also provide a return values as a :py:class:`metatensor.TensorMap` in
 :ref:`metatensor`.
 
-.. automodule:: meshlode.calculators
-
 .. toctree::
    :maxdepth: 1
+   :glob:
 
-   meshpotential
+   ./*

docs/src/references/calculators/pmepotential.rst

@@ -0,0 +1,6 @@
+PMEPotential
+############
+
+.. autoclass:: meshlode.PMEPotential
+    :members:
+    :undoc-members:

docs/src/references/lib/index.rst

@@ -6,6 +6,6 @@ are used for the meshLODE calculators.
 
 .. toctree::
    :maxdepth: 1
+   :glob:
 
-   fourier_convolution
-   mesh_interpolator
+   ./*

docs/src/references/lib/kvectors.rst

@@ -0,0 +1,6 @@
+Kvectors
+========
+
+.. automodule:: meshlode.lib.kvectors
+    :members:
+    :undoc-members:

docs/src/references/metatensor/directpotential.rst

@@ -0,0 +1,6 @@
+DirectPotential
+###############
+
+.. autoclass:: meshlode.metatensor.DirectPotential
+    :members:
+    :undoc-members:

docs/src/references/metatensor/ewaldpotential.rst

@@ -0,0 +1,6 @@
+EwaldPotential
+##############
+
+.. autoclass:: meshlode.metatensor.EwaldPotential
+    :members:
+    :undoc-members:

docs/src/references/metatensor/index.rst

@@ -15,5 +15,6 @@ For a plain :py:class:`torch.Tensor` refer to :ref:`calculators`.
 
 .. toctree::
    :maxdepth: 1
+   :glob:
 
-   meshpotential
+   ./*

docs/src/references/metatensor/pmepotential.rst

@@ -0,0 +1,6 @@
+PMEPotential
+############
+
+.. autoclass:: meshlode.metatensor.PMEPotential
+    :members:
+    :undoc-members:

examples/neighborlist_example.py

@@ -0,0 +1,123 @@
+"""
+Computations with explicit Neighbor Lists
+=========================================
+
+This example will explain how to use the metatensor branch of Meshlode with an attached
+neighborlist to an :py:class:`metatensor.torch.atomistic.System` object.
+"""
+
+# %%
+
+import math
+
+import numpy as np
+import torch
+from ase import Atoms
+from ase.neighborlist import neighbor_list
+from metatensor.torch import Labels, TensorBlock
+from metatensor.torch.atomistic import NeighborListOptions, System
+
+import meshlode
+
+
+# %%
+# Define simple example structure having the CsCl structure and compute the reference
+# values. PMEPotential by default outputs the types sorted according to the atomic
+# number. Thus, we input the compound "CsCl" and "ClCs" since Cl and Cs have atomic
+# numbers 17 and 55, respectively.
+
+types = torch.tensor([17, 55])  # Cl and Cs
+positions = torch.tensor([[0, 0, 0], [0.5, 0.5, 0.5]])
+charges = torch.tensor([-1.0, 1.0]).reshape(-1, 1)
+cell = torch.eye(3)
+
+# %%
+# Define the expected values of the energy
+
+n_atoms = len(types)
+madelung = 2 * 1.7626 / math.sqrt(3)
+energies_ref = -madelung * torch.ones((n_atoms, 1))
+
+# %%
+# We first define general parameters for our calculation MeshLODE.
+
+atomic_smearing = 0.1
+cell = torch.eye(3)
+mesh_spacing = atomic_smearing / 4
+interpolation_order = 2
+
+
+# %%
+# Generate neighbor list using ASE's :py:func:`neighbor_list()
+# <ase.neighborlist.neighbor_list>` function.
+
+sr_cutoff = np.sqrt(3) * 0.8
+struc = Atoms(positions=positions, cell=cell, pbc=True)
+nl_i, nl_j, nl_S, nl_D = neighbor_list("ijSD", struc, sr_cutoff)
+
+
+# %%
+# Convert ASE neighbor list into suitable format for a Metatensor system.
+
+neighbors = TensorBlock(
+    values=torch.from_numpy(nl_D.astype(np.float32).reshape(-1, 3, 1)),
+    samples=Labels(
+        names=[
+            "first_atom",
+            "second_atom",
+            "cell_shift_a",
+            "cell_shift_b",
+            "cell_shift_c",
+        ],
+        values=torch.from_numpy(np.vstack([nl_i, nl_j, nl_S.T]).T),
+    ),
+    components=[Labels.range("xyz", 3)],
+    properties=Labels.range("distance", 1),
+)
+
+
+# %%
+# Define the system.
+
+system = System(types=types, positions=positions, cell=cell)
+
+# %%
+# Attach charges to the system.
+
+data = TensorBlock(
+    values=charges,
+    samples=Labels.range("atom", charges.shape[0]),
+    components=[],
+    properties=Labels.range("charge", charges.shape[1]),
+)
+system.add_data(name="charges", data=data)
+
+# %%
+# Attach ``neighbors`` to ``system`` object.
+
+nl_options = NeighborListOptions(cutoff=sr_cutoff, full_list=True)
+system.add_neighbor_list(options=nl_options, neighbors=neighbors)
+
+pme = meshlode.metatensor.PMEPotential(
+    atomic_smearing=atomic_smearing,
+    mesh_spacing=mesh_spacing,
+    interpolation_order=interpolation_order,
+    subtract_self=True,
+    sr_cutoff=sr_cutoff,
+)
+potential = pme.compute(system)
+
+# %%
+# The total energy is just the sum of all atomic energies
+
+print(potential)
+
+# total_energy_metatensor = torch.sum(potential[0].values)
+
+# # %%
+# # Compare against reference Madelung constant and reference energy:
+
+# print("Using the metatensor version")
+# print(f"Computed energies on each atom = {potential[0].values.tolist()}")
+# print(f"Reference Madelung constant = {madelung:.3f}")
+# print(f"Total energy = {total_energy_metatensor[0].values}")