Skip to content

Commit

Permalink
Remove own System class
Browse files Browse the repository at this point in the history
  • Loading branch information
PicoCentauri committed Mar 4, 2024
1 parent 0bedb0a commit 31a3c38
Show file tree
Hide file tree
Showing 13 changed files with 274 additions and 158 deletions.
1 change: 0 additions & 1 deletion docs/src/references/lib/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -9,4 +9,3 @@ are used for the meshLODE calculators.

fourier_convolution
mesh_interpolator
system
6 changes: 0 additions & 6 deletions docs/src/references/lib/system.rst

This file was deleted.

37 changes: 15 additions & 22 deletions examples/library-tutorial.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
"""
Basic Tutorial for Library functions
====================================
This examples provides an illustration of the functioning of the underlaying library
functions of ``meshlode`` and the construction LODE descriptors (`Grisafi 2019
<https://doi.org/10.1063/1.5128375>`__, `Grisafi 2021
Expand All @@ -18,6 +17,7 @@
import matplotlib.pyplot as plt
import numpy as np
import torch
from metatensor.torch.atomistic import System

import meshlode

Expand Down Expand Up @@ -48,22 +48,19 @@ def sliceplot(mesh, sz=12, cmap="viridis", vmin=None, vmax=None):
# %%
# Builds the structure
# --------------------
#
# Builds a CsCl structure by replicating the primitive cell using ase and convert it to
# a :py:class:`list` of :py:class:`meshlode.System`. We add a bit of noise to make
# it less boring!
# a :py:class:`list` of :py:class:`metatensor.torch.atomistic.System`. We add a bit of
# noise to make it less boring!
#

positions = torch.tensor([[0, 0, 0], [0.5, 0.5, 0.5]]) * 4
atomic_types = torch.tensor([55, 17]) # Cs and Cl
types = torch.tensor([55, 17]) # Cs and Cl
cell = torch.eye(3) * 4
ase_frame = ase.Atoms(positions=positions, cell=cell, numbers=atomic_types).repeat(
[2, 2, 2]
)
ase_frame = ase.Atoms(positions=positions, cell=cell, numbers=types).repeat([2, 2, 2])
ase_frame.positions[:] += np.random.normal(size=ase_frame.positions.shape) * 0.1
charges = torch.tensor([1.0, -1.0] * 8)
frame = meshlode.System(
species=torch.tensor(ase_frame.numbers),
system = System(
types=torch.tensor(ase_frame.numbers),
positions=torch.tensor(np.array(ase_frame.positions)),
cell=torch.tensor(ase_frame.cell),
)
Expand All @@ -85,7 +82,6 @@ def sliceplot(mesh, sz=12, cmap="viridis", vmin=None, vmax=None):
# %%
# MeshInterpolator
# ----------------
#
# ``MeshInterpolator`` serves as a utility class to compute a mesh
# representation of points, and/or to project a function defined on the
# mesh on a set of points. Computing the mesh representation is a two-step
Expand All @@ -95,15 +91,15 @@ def sliceplot(mesh, sz=12, cmap="viridis", vmin=None, vmax=None):
#

interpol = meshlode.lib.mesh_interpolator.MeshInterpolator(
frame.cell, torch.tensor([16, 16, 16]), interpolation_order=3
system.cell, torch.tensor([16, 16, 16]), interpolation_order=3
)

interpol.compute_interpolation_weights(frame.positions)
interpol.compute_interpolation_weights(system.positions)


# %%
# We use two sets of weights: ones (giving the atom density irrespective
# of the species) and charges (giving a smooth representation of the point
# of the types) and charges (giving a smooth representation of the point
# charges).
#

Expand All @@ -126,7 +122,6 @@ def sliceplot(mesh, sz=12, cmap="viridis", vmin=None, vmax=None):
# %%
# Fourier filter
# --------------
#
# This module computes a Fourier-domain filter, that can be used e.g. to
# smear the density and/or compute a 1/r^p potential field. This can also
# be easily extended to compute an arbitrary filter
Expand All @@ -137,15 +132,15 @@ def sliceplot(mesh, sz=12, cmap="viridis", vmin=None, vmax=None):
# %%
# plain atomic_smearing
rho_mesh = fsc.compute(
mesh_values=mesh, cell=frame.cell, potential_exponent=0, atomic_smearing=1
mesh_values=mesh, cell=system.cell, potential_exponent=0, atomic_smearing=1
)

sliceplot(rho_mesh[0, :, :, :5])

# %%
# coulomb-like potential, no atomic_smearing
coulomb_mesh = fsc.compute(
mesh_values=mesh, cell=frame.cell, potential_exponent=1, atomic_smearing=0
mesh_values=mesh, cell=system.cell, potential_exponent=1, atomic_smearing=0
)

sliceplot(coulomb_mesh[1, :, :, :5], cmap="seismic")
Expand All @@ -154,7 +149,6 @@ def sliceplot(mesh, sz=12, cmap="viridis", vmin=None, vmax=None):
# %%
# Back-interpolation (on the same points)
# ---------------------------------------
#
# The same ``MeshInterpolator`` object can be used to compute a field on
# the same points used initially to generate the atom density
#
Expand All @@ -167,21 +161,20 @@ def sliceplot(mesh, sz=12, cmap="viridis", vmin=None, vmax=None):
# %%
# Back-interpolation (on different points)
# ----------------------------------------
#
# In order to compute the field on a different set of points, it is
# sufficient to build another ``MeshInterpolator`` object and to compute
# it with the desired field. One can also use a different
# ``interpolation_order``, if wanted.
#

interpol_slice = meshlode.lib.mesh_interpolator.MeshInterpolator(
frame.cell, torch.tensor([16, 16, 16]), interpolation_order=4
system.cell, torch.tensor([16, 16, 16]), interpolation_order=4
)

# Compute a denser grid on a 2D slice
n_points = 50
x = torch.linspace(0, frame.cell[0, 0], n_points + 1)[:n_points]
y = torch.linspace(0, frame.cell[1, 1], n_points + 1)[:n_points]
x = torch.linspace(0, system.cell[0, 0], n_points + 1)[:n_points]
y = torch.linspace(0, system.cell[1, 1], n_points + 1)[:n_points]
xx, yy = torch.meshgrid(x, y, indexing="ij")

# Flatten xx and yy, and concatenate with a zero column for the z-coordinate
Expand Down
33 changes: 17 additions & 16 deletions examples/madelung.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
"""
Compute Madelung Constants
==========================
In this tutorial we show how to calculate the Madelung constants and total electrostatic
energy of atomic structures using the :py:class:`meshlode.MeshPotential` and
:py:class:`meshlode.metatensor.MeshPotential` calculator.
Expand All @@ -11,24 +10,24 @@
import math

import torch
from metatensor.torch.atomistic import System

import meshlode


# %%
# Define simple example structure having the CsCl structure and compute the reference
# values. MeshPotential by default outputs the species sorted according to the atomic
# values. MeshPotential 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.
atomic_types = torch.tensor([17, 55]) # Cl and Cs
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])
cell = torch.eye(3)
positions = torch.tensor([[0, 0, 0], [0.5, 0.5, 0.5]])
frame = meshlode.System(species=atomic_types, positions=positions, cell=torch.eye(3))

# %%
# Define the expected values of the energy
n_atoms = len(positions)
n_atoms = len(types)
madelung = 2 * 1.7626 / math.sqrt(3)
energies_ref = -madelung * torch.ones((n_atoms, 1))

Expand All @@ -43,7 +42,6 @@
# %%
# Computation using ``meshlode``
# ------------------------------
#
# Compute features using

MP = meshlode.MeshPotential(
Expand All @@ -52,7 +50,7 @@
interpolation_order=interpolation_order,
subtract_self=True,
)
potentials_torch = MP.compute(frame)
potentials_torch = MP.compute(types=types, positions=positions, cell=cell)

# %%
# The "potentials" that have been computed so far are not the actual electrostatic
Expand All @@ -66,7 +64,7 @@
for idx_n in range(n_atoms):
# The coulomb potential between atoms i and j is charge_i * charge_j / d_ij
# The features are simply computing a pure 1/r potential with no prefactors.
# Thus, to compute the energy between atoms of species i and j, we need to
# Thus, to compute the energy between atoms of types i and j, we need to
# multiply by the charges of i and j.
print(charges[idx_c] * charges[idx_n], potentials_torch[idx_n, idx_c])
atomic_energies_torch[idx_c] += (
Expand All @@ -88,33 +86,36 @@
# %%
# Computation using ``meshlode.metatensor``
# -----------------------------------------
#
# We now compute the same constants using the metatensor based calculator
# We now compute the same constants using the metatensor based calculator. To achieve
# this we first store our system parameters like the ``types``, ``positions`` and the
# ``cell`` defined above into a :py:class:`metatensor.torch.atomistic.System` class.

system = System(types=types, positions=positions, cell=cell)

MP = meshlode.metatensor.MeshPotential(
atomic_smearing=atomic_smearing,
mesh_spacing=mesh_spacing,
interpolation_order=interpolation_order,
subtract_self=True,
)
potential_metatensor = MP.compute(frame)
potential_metatensor = MP.compute(system)


# %%
# To get the Madelung constant, we again need to take a linear combination
# of the "potentials" weighted by the charges of the atoms.

atomic_energies_metatensor = torch.zeros((n_atoms, 1))
for idx_c, c in enumerate(atomic_types):
for idx_n, n in enumerate(atomic_types):
# Take the coefficients with the correct center atom and neighbor atom species
for idx_c, c in enumerate(types):
for idx_n, n in enumerate(types):
# Take the coefficients with the correct center atom and neighbor atom types
block = potential_metatensor.block(
{"center_type": int(c), "neighbor_type": int(n)}
)

# The coulomb potential between atoms i and j is charge_i * charge_j / d_ij
# The features are simply computing a pure 1/r potential with no prefactors.
# Thus, to compute the energy between atoms of species i and j, we need to
# Thus, to compute the energy between atoms of types i and j, we need to
# multiply by the charges of i and j.
print(c, n, charges[idx_c] * charges[idx_n], block.values[0, 0])
atomic_energies_metatensor[idx_c] += (
Expand Down
4 changes: 2 additions & 2 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ keywords = [
"Atomistic Simulations",
]
dependencies = [
"torch >= 1.11",
"torch >=1.11",
]
dynamic = ["version"]

Expand All @@ -44,7 +44,7 @@ examples = [
"matplotlib",
]
metatensor = [
"metatensor[torch]",
"metatensor-torch >=0.3",
]

[project.urls]
Expand Down
3 changes: 1 addition & 2 deletions src/meshlode/__init__.py
Original file line number Diff line number Diff line change
@@ -1,11 +1,10 @@
from .calculators.meshpotential import MeshPotential
from .lib.system import System

try:
from . import metatensor # noqa
except ImportError:
pass


__all__ = ["MeshPotential", "System"]
__all__ = ["MeshPotential"]
__version__ = "0.0.0-dev"
Loading

0 comments on commit 31a3c38

Please sign in to comment.