Note
Go to the end to download the full example code.
The PET-MAD universal potential¶
- Authors:
Philip Loche @PicoCentauri, Michele Ceriotti @ceriottm, Arslan Mazitov @abmazitov
This example demonstrates how to use the PET-MAD model with ASE, i-PI and LAMMPS. PET-MAD is a “universal” machine-learning forcefield trained on a dataset that aims to incorporate a very high degree of structural diversity.
The point-edge transformer (PET) is an unconstrained architecture that achieves a high degree of symmetry compliance through data augmentation during training (see the PET paper for more details). The unconstrained nature of the model simplifies its implementation and structure, making it computationally efficient and very expressive.
The MAD dataset combines “stable” inorganic structures from the MC3D dataset, 2D structures from the MC2D dataset, and molecular crystals from the ShiftML dataset with “Maximum Atomic Diversity” configurations, generated by distorting the composition and structure of these stable templates. By doing so, PET-MAD achieves state-of-the-art accuracy despite the MAD dataset containing fewer than 100k structures. The reference DFT settings are highly converged, but limited to a PBEsol functional, so the accuracy against experimental data depends on how good this level of theory is for a given system. PET-MAD is introduced, and benchmarked for several challenging modeling tasks, in this preprint.
Start by importing the required libraries. To use PET-MAD, and obtain all the necessary dependencies, you can simply use pip to install the PET-MAD package:
pip install git+https://github.com/lab-cosmo/pet-mad.git
import os
import subprocess
from copy import copy, deepcopy
# ASE and i-PI scripting utilities
import ase.units
import chemiscope
import matplotlib.pyplot as plt
import numpy as np
import requests
from ase.optimize import LBFGS
from ipi.utils.mathtools import get_rotation_quadrature_lebedev
from ipi.utils.parsing import read_output, read_trajectory
from ipi.utils.scripting import (
InteractiveSimulation,
forcefield_xml,
motion_nvt_xml,
simulation_xml,
)
# pet-mad ASE calculator
from pet_mad.calculator import PETMADCalculator
if hasattr(__import__("builtins"), "get_ipython"):
get_ipython().run_line_magic("matplotlib", "inline") # noqa: F821
Inference on the MAD test set¶
We begin by using the ase-compatible calculator to evaluate energy and forces for
a test dataset that contains both hold-out structures from the MAD dataset, and a few
structures from popular datasets (MPtrj,
Alexandria,
OC2020,
SPICE,
MD22) re-computed
with consistent DFT settings.
Load the dataset¶
We fetch the dataset, and load only some of the structures, to speed up the example runtime on CPU. The model can also run (much faster) on GPUs if you have some at hand.
filename = "data/mad-test-mad-settings.xyz"
if not os.path.exists(filename):
url = (
"https://huggingface.co/lab-cosmo/pet-mad/resolve/"
"main/benchmarks/mad-test-mad-settings.xyz"
)
response = requests.get(url)
response.raise_for_status()
with open(filename, "wb") as f:
f.write(response.content)
test_structures = ase.io.read(filename, "::16")
# also extract reference energetics and metadata
test_energy = []
test_forces = []
test_natoms = []
test_origin = []
subsets = []
for s in test_structures:
test_energy.append(s.get_potential_energy())
test_natoms.append(len(s))
test_forces.append(s.get_forces())
test_origin.append(s.info["origin"])
if s.info["origin"] not in subsets:
subsets.append(s.info["origin"])
test_natoms = np.array(test_natoms)
test_origin = np.array(test_origin)
test_energy = np.array(test_energy)
test_forces = np.array(test_forces, dtype=object)
Single point energy and forces¶
PET-MAD is compatible with the metatensor atomistic models interface which allows us to run it with ASE and many other MD engines. For more details see the metatensor documentation.
We now load the PET-MAD ASE calculator and calculate energy and forces.
calculator = PETMADCalculator(version="latest", device="cpu")
Here, we run the computation on the CPU. If you have a CUDA GPU you can also set
device="cuda" to speed up the computation.
mad_energy = []
mad_forces = []
mad_structures = []
for structure in test_structures:
tmp = deepcopy(structure)
tmp.calc = copy(calculator) # avoids ovewriting results.
mad_energy.append(tmp.get_potential_energy())
mad_forces.append(tmp.get_forces())
mad_structures.append(tmp)
mad_energy = np.array(mad_energy)
mad_forces = np.array(mad_forces, dtype=object)
A parity plot with the model predictions
tab10 = plt.get_cmap("tab10")
fig, ax = plt.subplots(1, 2, figsize=(6, 3), constrained_layout=True)
ax[0].plot([0, 1], [0, 1], "b:", transform=ax[0].transAxes)
ax[1].plot([0, 1], [0, 1], "b:", transform=ax[1].transAxes)
for i, sub in enumerate(subsets):
sel = np.where(test_origin == sub)[0]
ax[0].plot(
mad_energy[sel] / test_natoms[sel],
test_energy[sel] / test_natoms[sel],
".",
c=tab10(i),
label=sub,
)
ax[1].plot(
np.concatenate(mad_forces[sel]).flatten(),
np.concatenate(test_forces[sel]).flatten(),
".",
c=tab10(i),
)
ax[0].set_xlabel("MAD energy / eV/atom")
ax[0].set_ylabel("Reference energy / eV/atom")
ax[1].set_xlabel("MAD forces / eV/Å")
ax[1].set_ylabel("Refrerence forces / eV/Å")
fig.legend(loc="upper center", bbox_to_anchor=(0.55, 1.20), ncol=3)
<matplotlib.legend.Legend object at 0x7ff792bf53a0>
Explore the dataset using chemiscope
chemiscope.show(
test_structures,
mode="default",
properties={
"origin": test_origin,
"energy_ref": test_energy / test_natoms,
"energy_mad": mad_energy / test_natoms,
"energy_error": np.abs((test_energy - mad_energy) / test_natoms),
"force_error": [
np.linalg.norm(f1 - f2) / n
for (f1, f2, n) in zip(mad_forces, test_forces, test_natoms)
],
},
shapes={
"forces_ref": chemiscope.ase_vectors_to_arrows(
mad_structures, "forces", scale=1.0
),
"forces_mad": chemiscope.ase_vectors_to_arrows(
test_structures, "forces", scale=1.0
),
},
settings=chemiscope.quick_settings(
x="energy_mad",
y="energy_ref",
symbol="origin",
structure_settings={"unitCell": True, "shape": ["forces_ref"]},
),
)
How about equivariance‽¶
The PET architecture does not provide “intrinsically” invariant energy predictions, but learns symmetry from data augmentation. Should you worry? The authors of PET-MAD certainly do, and they have studied extensively whether the symmetry breaking can cause serious artefacts. You can check by yourself following the procedure below, that evaluates a structures over a grid of rotations, estimating the variability in energy (which is around 1meV/atom, much smaller than the test error).
rotations = get_rotation_quadrature_lebedev(3)
rot_test = test_structures[100]
rot_structures = []
rot_weights = []
rot_energies = []
rot_forces = []
rot_angles = []
for rot, w, angles in rotations:
tmp = rot_test.copy()
tmp.positions = tmp.positions @ rot.T
tmp.cell = tmp.cell @ rot.T
tmp.calc = copy(calculator)
rot_weights.append(w)
rot_energies.append(tmp.get_potential_energy() / len(tmp))
rot_forces.append(tmp.get_forces())
rot_structures.append(tmp)
rot_angles.append(angles)
rot_energies = np.array(rot_energies)
rot_weights = np.array(rot_weights)
rot_angles = np.array(rot_angles)
erot_rms = 1e3 * np.sqrt(
np.sum(rot_energies**2 * rot_weights) / np.sum(rot_weights)
- (np.sum(rot_energies * rot_weights) / np.sum(rot_weights)) ** 2
)
erot_max = 1e3 * np.abs(rot_energies.max() - rot_energies.min())
print(
f"""
Symmetry breaking, energy:
RMS: {erot_rms:.3f} meV/at.
Max: {erot_max:.3f} meV/at.
"""
)
Symmetry breaking, energy:
RMS: 0.405 meV/at.
Max: 1.574 meV/at.
You can also inspect the rotational behavior visually
chemiscope.show(
rot_structures,
mode="default",
properties={
"delta_energy": 1e3 * (rot_energies - rot_energies.mean()),
"euler_angles": rot_angles,
},
shapes={
"forces": chemiscope.ase_vectors_to_arrows(rot_structures, "forces", scale=4.0),
},
settings=chemiscope.quick_settings(
x="euler_angles[1]",
y="euler_angles[2]",
z="euler_angles[3]",
color="delta_energy",
structure_settings={"unitCell": True, "shape": ["forces"]},
),
)
Note also that i-PI provides functionalities to do this automatically to obtain MD trajectories with a even higher degree of symmetry-compliance.
Simulating a complex surface¶
PET-MAD is designed to be robust and stable when executing sophisticated modeling workflows. As an example, we consider a slab of an Al-6xxx alloy (aluminum with a few percent Mg and Si) with some oxygen molecules adsorbed at the (111) surface.
Warning
The overall Si+Mg concentration in an Al6xxx alloy is far lower than what depicted here. This is just a demonstrative example and should not be taken as the starting point of a serious study of this system.
al_surface = ase.io.read("data/al6xxx-o2.xyz")
Geometry optimization with ASE¶
As a first example, we use the ase geometry LBFGS optimizer to relax the initial
positions. This leads to the rapid decomposition of the oxygen molecules and the
formation of an oxide layer.
atoms = al_surface.copy()
atoms.calc = calculator
opt = LBFGS(atoms)
traj_atoms = []
traj_energy = []
opt.attach(lambda: traj_atoms.append(atoms.copy()))
opt.attach(lambda: traj_energy.append(atoms.get_potential_energy()))
# stop the optimization early to speed up the example
opt.run(fmax=0.001, steps=30)
Step Time Energy fmax
LBFGS: 0 09:36:43 -953.123230 3.732129
LBFGS: 1 09:36:45 -956.171509 2.793623
LBFGS: 2 09:36:47 -956.024292 14.278265
LBFGS: 3 09:36:49 -959.920532 1.655591
LBFGS: 4 09:36:52 -960.696655 1.516880
LBFGS: 5 09:36:54 -961.391663 2.048758
LBFGS: 6 09:36:56 -961.777954 1.241909
LBFGS: 7 09:36:58 -962.608337 1.428242
LBFGS: 8 09:37:00 -962.762634 1.032299
LBFGS: 9 09:37:02 -963.099060 1.106594
LBFGS: 10 09:37:05 -963.314758 1.698585
LBFGS: 11 09:37:07 -963.743774 1.655965
LBFGS: 12 09:37:09 -964.083374 1.553182
LBFGS: 13 09:37:11 -964.661316 1.819931
LBFGS: 14 09:37:14 -965.175720 1.762308
LBFGS: 15 09:37:16 -965.520203 1.330094
LBFGS: 16 09:37:18 -965.903687 1.197484
LBFGS: 17 09:37:20 -966.199280 1.288282
LBFGS: 18 09:37:23 -966.645874 1.403857
LBFGS: 19 09:37:25 -967.111084 2.378607
LBFGS: 20 09:37:28 -967.621338 2.800260
LBFGS: 21 09:37:30 -968.812317 3.187743
LBFGS: 22 09:37:32 -969.942993 5.129991
LBFGS: 23 09:37:35 -971.183167 7.703203
LBFGS: 24 09:37:37 -973.405884 5.725447
LBFGS: 25 09:37:39 -975.727295 3.588656
LBFGS: 26 09:37:42 -977.515137 3.823957
LBFGS: 27 09:37:44 -979.349792 6.409081
LBFGS: 28 09:37:47 -981.495544 3.498837
LBFGS: 29 09:37:49 -982.948303 4.598636
LBFGS: 30 09:37:51 -985.354797 5.147579
np.False_
Even if the optimization is cut short and far from converged, the decomposition of the oxygen molecules is apparent, and leads to a large energetic stabilization
chemiscope.show(
frames=traj_atoms,
properties={
"index": np.arange(0, len(traj_atoms)),
"energy": traj_energy,
},
mode="default",
settings=chemiscope.quick_settings(trajectory=True),
)
Molecular dynamics with atoms exchange with i-PI¶
The geometry optimization shows the high reactivity of this surface, but does not properly account for finite temperature and does not sample the diffusion of solute atoms in the alloy (which is mediated by vacancies).
We use i-PI to perform a molecular dynamics trajectory at 800K, combined with Monte Carlo steps that swap the nature of atoms, allowing the simulation to reach equilibrium in the solute-atoms distributions without having to introduce vacancies or wait for the very long time scale needed for diffusion.
Before starting the simulations with MD engines, it is important
to export the model to a format that can be used by the engine.
This is done by saving the model to a file, which includes the
model weights and the compiled extensions. We use the collect_extensions
argument to save the compiled extensions to disk. These extensions ensure that
the model remains self-contained and can be executed without requiring the
original Python or C++ source code. In particular, this is necessary for
the LAMMPS interface to work because it has no access to the Python code.
calculator.model.save("pet-mad-latest.pt", collect_extensions="extensions")
/home/runner/work/atomistic-cookbook/atomistic-cookbook/.nox/pet-mad/lib/python3.12/site-packages/metatensor/torch/atomistic/model.py:522: UserWarning: No length unit was provided for the model.
warnings.warn(
The behavior of i-PI is controlled by an XML input file.
The utils.scripting module contains several helper
functions to generate the basic components.
Here we use a <motion mode="multi"> block to combine
a MD run with a <motion mode="atomswap"> block that
attemts swapping atoms, with a Monte Carlo acceptance.
motion_xml = f"""
<motion mode="multi">
{motion_nvt_xml(timestep=5.0 * ase.units.fs)}
<motion mode="atomswap">
<atomswap>
<names> [ Al, Si, Mg, O] </names>
</atomswap>
</motion>
</motion>
"""
input_xml = simulation_xml(
structures=[al_surface],
forcefield=forcefield_xml(
name="pet-mad",
mode="direct",
pes="metatensor",
parameters={"model": "pet-mad-latest.pt", "template": "data/al6xxx-o2.xyz"},
),
motion=motion_xml,
temperature=800,
verbosity="low",
prefix="nvt_atomxc",
)
print(input_xml)
<simulation verbosity='low' safe_stride='20'>
<ffdirect name='pet-mad'>
<pes>metatensor</pes>
<parameters>{model: pet-mad-latest.pt, template: data/al6xxx-o2.xyz}</parameters>
</ffdirect>
<output prefix='nvt_atomxc'>
<properties stride='2' filename='out'>[ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, potential{electronvolt} ]</properties>
<trajectory filename='pos' stride='20' cell_units='angstrom'>positions{angstrom}</trajectory>
<checkpoint stride='200'>
</checkpoint>
</output>
<system>
<beads natoms='241' nbeads='1'>
<q shape='(1, 723)'>
[ 2.70588227e+00, 1.56224186e+00, 1.88972613e+01, 8.11764682e+00, 1.56224186e+00,
1.88972613e+01, 1.35294114e+01, 1.56224186e+00, 1.88972613e+01, 1.89411759e+01,
1.56224186e+00, 1.88972613e+01, 2.43529404e+01, 1.56224186e+00, 1.88972613e+01,
5.41176453e+00, 6.24896742e+00, 1.88972613e+01, 1.08235291e+01, 6.24896742e+00,
1.88972613e+01, 1.62352936e+01, 6.24896742e+00, 1.88972613e+01, 2.16470582e+01,
6.24896742e+00, 1.88972613e+01, 2.70588227e+01, 6.24896742e+00, 1.88972613e+01,
8.11764682e+00, 1.09356930e+01, 1.88972613e+01, 1.35294114e+01, 1.09356930e+01,
1.88972613e+01, 1.89411759e+01, 1.09356930e+01, 1.88972613e+01, 2.43529404e+01,
1.09356930e+01, 1.88972613e+01, 2.97647050e+01, 1.09356930e+01, 1.88972613e+01,
1.08235291e+01, 1.56224186e+01, 1.88972613e+01, 1.62352936e+01, 1.56224186e+01,
1.88972613e+01, 2.16470582e+01, 1.56224186e+01, 1.88972613e+01, 2.70588227e+01,
1.56224186e+01, 1.88972613e+01, 3.24705872e+01, 1.56224186e+01, 1.88972613e+01,
1.35294114e+01, 2.03091441e+01, 1.88972613e+01, 1.89411759e+01, 2.03091441e+01,
1.88972613e+01, 2.43529404e+01, 2.03091441e+01, 1.88972613e+01, 2.97647050e+01,
2.03091441e+01, 1.88972613e+01, 3.51764695e+01, 2.03091441e+01, 1.88972613e+01,
-0.00000000e+00, 3.12448372e+00, 2.33159485e+01, 5.41176453e+00, 3.12448372e+00,
2.33159485e+01, 1.08235291e+01, 3.12448372e+00, 2.33159485e+01, 1.62352936e+01,
3.12448372e+00, 2.33159485e+01, 2.16470582e+01, 3.12448372e+00, 2.33159485e+01,
2.70588227e+00, 7.81120928e+00, 2.33159485e+01, 8.11764682e+00, 7.81120928e+00,
2.33159485e+01, 1.35294114e+01, 7.81120928e+00, 2.33159485e+01, 1.89411759e+01,
7.81120928e+00, 2.33159485e+01, 2.43529404e+01, 7.81120928e+00, 2.33159485e+01,
5.41176453e+00, 1.24979349e+01, 2.33159485e+01, 1.08235291e+01, 1.24979349e+01,
2.33159485e+01, 1.62352936e+01, 1.24979349e+01, 2.33159485e+01, 2.16470582e+01,
1.24979349e+01, 2.33159485e+01, 2.70588227e+01, 1.24979349e+01, 2.33159485e+01,
8.11764682e+00, 1.71846604e+01, 2.33159485e+01, 1.35294114e+01, 1.71846604e+01,
2.33159485e+01, 1.89411759e+01, 1.71846604e+01, 2.33159485e+01, 2.43529404e+01,
1.71846604e+01, 2.33159485e+01, 2.97647050e+01, 1.71846604e+01, 2.33159485e+01,
1.08235291e+01, 2.18713860e+01, 2.33159485e+01, 1.62352936e+01, 2.18713860e+01,
2.33159485e+01, 2.16470582e+01, 2.18713860e+01, 2.33159485e+01, 2.70588227e+01,
2.18713860e+01, 2.33159485e+01, 3.24705872e+01, 2.18713860e+01, 2.33159485e+01,
0.00000000e+00, 0.00000000e+00, 2.77346357e+01, 5.41176453e+00, 0.00000000e+00,
2.77346357e+01, 1.08235291e+01, 0.00000000e+00, 2.77346357e+01, 1.62352936e+01,
0.00000000e+00, 2.77346357e+01, 2.16470582e+01, 0.00000000e+00, 2.77346357e+01,
2.70588227e+00, 4.68672556e+00, 2.77346357e+01, 8.11764682e+00, 4.68672556e+00,
2.77346357e+01, 1.35294114e+01, 4.68672556e+00, 2.77346357e+01, 1.89411759e+01,
4.68672556e+00, 2.77346357e+01, 2.43529404e+01, 4.68672556e+00, 2.77346357e+01,
5.41176453e+00, 9.37345114e+00, 2.77346357e+01, 1.08235291e+01, 9.37345114e+00,
2.77346357e+01, 1.62352936e+01, 9.37345114e+00, 2.77346357e+01, 2.16470582e+01,
9.37345114e+00, 2.77346357e+01, 2.70588227e+01, 9.37345114e+00, 2.77346357e+01,
8.11764682e+00, 1.40601767e+01, 2.77346357e+01, 1.35294114e+01, 1.40601767e+01,
2.77346357e+01, 1.89411759e+01, 1.40601767e+01, 2.77346357e+01, 2.43529404e+01,
1.40601767e+01, 2.77346357e+01, 2.97647050e+01, 1.40601767e+01, 2.77346357e+01,
1.08235291e+01, 1.87469023e+01, 2.77346357e+01, 1.62352936e+01, 1.87469023e+01,
2.77346357e+01, 2.16470582e+01, 1.87469023e+01, 2.77346357e+01, 2.70588227e+01,
1.87469023e+01, 2.77346357e+01, 3.24705872e+01, 1.87469023e+01, 2.77346357e+01,
2.70588227e+00, 1.56224186e+00, 3.21533230e+01, 8.11764682e+00, 1.56224186e+00,
3.21533230e+01, 1.35294114e+01, 1.56224186e+00, 3.21533230e+01, 1.89411759e+01,
1.56224186e+00, 3.21533230e+01, 2.43529404e+01, 1.56224186e+00, 3.21533230e+01,
5.41176453e+00, 6.24896742e+00, 3.21533230e+01, 1.08235291e+01, 6.24896742e+00,
3.21533230e+01, 1.62352936e+01, 6.24896742e+00, 3.21533230e+01, 2.16470582e+01,
6.24896742e+00, 3.21533230e+01, 2.70588227e+01, 6.24896742e+00, 3.21533230e+01,
8.11764682e+00, 1.09356930e+01, 3.21533230e+01, 1.35294114e+01, 1.09356930e+01,
3.21533230e+01, 1.89411759e+01, 1.09356930e+01, 3.21533230e+01, 2.43529404e+01,
1.09356930e+01, 3.21533230e+01, 2.97647050e+01, 1.09356930e+01, 3.21533230e+01,
1.08235291e+01, 1.56224186e+01, 3.21533230e+01, 1.62352936e+01, 1.56224186e+01,
3.21533230e+01, 2.16470582e+01, 1.56224186e+01, 3.21533230e+01, 2.70588227e+01,
1.56224186e+01, 3.21533230e+01, 3.24705872e+01, 1.56224186e+01, 3.21533230e+01,
1.35294114e+01, 2.03091441e+01, 3.21533230e+01, 1.89411759e+01, 2.03091441e+01,
3.21533230e+01, 2.43529404e+01, 2.03091441e+01, 3.21533230e+01, 2.97647050e+01,
2.03091441e+01, 3.21533230e+01, 3.51764695e+01, 2.03091441e+01, 3.21533230e+01,
-0.00000000e+00, 3.12448372e+00, 3.65720102e+01, 5.41176453e+00, 3.12448372e+00,
3.65720102e+01, 1.08235291e+01, 3.12448372e+00, 3.65720102e+01, 1.62352936e+01,
3.12448372e+00, 3.65720102e+01, 2.16470582e+01, 3.12448372e+00, 3.65720102e+01,
2.70588227e+00, 7.81120928e+00, 3.65720102e+01, 8.11764682e+00, 7.81120928e+00,
3.65720102e+01, 1.35294114e+01, 7.81120928e+00, 3.65720102e+01, 1.89411759e+01,
7.81120928e+00, 3.65720102e+01, 2.43529404e+01, 7.81120928e+00, 3.65720102e+01,
5.41176453e+00, 1.24979349e+01, 3.65720102e+01, 1.08235291e+01, 1.24979349e+01,
3.65720102e+01, 1.62352936e+01, 1.24979349e+01, 3.65720102e+01, 2.16470582e+01,
1.24979349e+01, 3.65720102e+01, 2.70588227e+01, 1.24979349e+01, 3.65720102e+01,
8.11764682e+00, 1.71846604e+01, 3.65720102e+01, 1.35294114e+01, 1.71846604e+01,
3.65720102e+01, 1.89411759e+01, 1.71846604e+01, 3.65720102e+01, 2.43529404e+01,
1.71846604e+01, 3.65720102e+01, 2.97647050e+01, 1.71846604e+01, 3.65720102e+01,
1.08235291e+01, 2.18713860e+01, 3.65720102e+01, 1.62352936e+01, 2.18713860e+01,
3.65720102e+01, 2.16470582e+01, 2.18713860e+01, 3.65720102e+01, 2.70588227e+01,
2.18713860e+01, 3.65720102e+01, 3.24705872e+01, 2.18713860e+01, 3.65720102e+01,
0.00000000e+00, 0.00000000e+00, 4.09906975e+01, 5.41176453e+00, 0.00000000e+00,
4.09906975e+01, 1.08235291e+01, 0.00000000e+00, 4.09906975e+01, 1.62352936e+01,
0.00000000e+00, 4.09906975e+01, 2.16470582e+01, 0.00000000e+00, 4.09906975e+01,
2.70588227e+00, 4.68672556e+00, 4.09906975e+01, 8.11764682e+00, 4.68672556e+00,
4.09906975e+01, 1.35294114e+01, 4.68672556e+00, 4.09906975e+01, 1.89411759e+01,
4.68672556e+00, 4.09906975e+01, 2.43529404e+01, 4.68672556e+00, 4.09906975e+01,
5.41176453e+00, 9.37345114e+00, 4.09906975e+01, 1.08235291e+01, 9.37345114e+00,
4.09906975e+01, 1.62352936e+01, 9.37345114e+00, 4.09906975e+01, 2.16470582e+01,
9.37345114e+00, 4.09906975e+01, 2.70588227e+01, 9.37345114e+00, 4.09906975e+01,
8.11764682e+00, 1.40601767e+01, 4.09906975e+01, 1.35294114e+01, 1.40601767e+01,
4.09906975e+01, 1.89411759e+01, 1.40601767e+01, 4.09906975e+01, 2.43529404e+01,
1.40601767e+01, 4.09906975e+01, 2.97647050e+01, 1.40601767e+01, 4.09906975e+01,
1.08235291e+01, 1.87469023e+01, 4.09906975e+01, 1.62352936e+01, 1.87469023e+01,
4.09906975e+01, 2.16470582e+01, 1.87469023e+01, 4.09906975e+01, 2.70588227e+01,
1.87469023e+01, 4.09906975e+01, 3.24705872e+01, 1.87469023e+01, 4.09906975e+01,
2.70588227e+00, 1.56224186e+00, 4.54093847e+01, 8.11764682e+00, 1.56224186e+00,
4.54093847e+01, 1.35294114e+01, 1.56224186e+00, 4.54093847e+01, 1.89411759e+01,
1.56224186e+00, 4.54093847e+01, 2.43529404e+01, 1.56224186e+00, 4.54093847e+01,
5.41176453e+00, 6.24896742e+00, 4.54093847e+01, 1.08235291e+01, 6.24896742e+00,
4.54093847e+01, 1.62352936e+01, 6.24896742e+00, 4.54093847e+01, 2.16470582e+01,
6.24896742e+00, 4.54093847e+01, 2.70588227e+01, 6.24896742e+00, 4.54093847e+01,
8.11764682e+00, 1.09356930e+01, 4.54093847e+01, 1.35294114e+01, 1.09356930e+01,
4.54093847e+01, 1.89411759e+01, 1.09356930e+01, 4.54093847e+01, 2.43529404e+01,
1.09356930e+01, 4.54093847e+01, 2.97647050e+01, 1.09356930e+01, 4.54093847e+01,
1.08235291e+01, 1.56224186e+01, 4.54093847e+01, 1.62352936e+01, 1.56224186e+01,
4.54093847e+01, 2.16470582e+01, 1.56224186e+01, 4.54093847e+01, 2.70588227e+01,
1.56224186e+01, 4.54093847e+01, 3.24705872e+01, 1.56224186e+01, 4.54093847e+01,
1.35294114e+01, 2.03091441e+01, 4.54093847e+01, 1.89411759e+01, 2.03091441e+01,
4.54093847e+01, 2.43529404e+01, 2.03091441e+01, 4.54093847e+01, 2.97647050e+01,
2.03091441e+01, 4.54093847e+01, 3.51764695e+01, 2.03091441e+01, 4.54093847e+01,
-0.00000000e+00, 3.12448372e+00, 4.98280720e+01, 5.41176453e+00, 3.12448372e+00,
4.98280720e+01, 1.08235291e+01, 3.12448372e+00, 4.98280720e+01, 1.62352936e+01,
3.12448372e+00, 4.98280720e+01, 2.16470582e+01, 3.12448372e+00, 4.98280720e+01,
2.70588227e+00, 7.81120928e+00, 4.98280720e+01, 8.11764682e+00, 7.81120928e+00,
4.98280720e+01, 1.35294114e+01, 7.81120928e+00, 4.98280720e+01, 1.89411759e+01,
7.81120928e+00, 4.98280720e+01, 2.43529404e+01, 7.81120928e+00, 4.98280720e+01,
5.41176453e+00, 1.24979349e+01, 4.98280720e+01, 1.08235291e+01, 1.24979349e+01,
4.98280720e+01, 1.62352936e+01, 1.24979349e+01, 4.98280720e+01, 2.16470582e+01,
1.24979349e+01, 4.98280720e+01, 2.70588227e+01, 1.24979349e+01, 4.98280720e+01,
8.11764682e+00, 1.71846604e+01, 4.98280720e+01, 1.35294114e+01, 1.71846604e+01,
4.98280720e+01, 1.89411759e+01, 1.71846604e+01, 4.98280720e+01, 2.43529404e+01,
1.71846604e+01, 4.98280720e+01, 2.97647050e+01, 1.71846604e+01, 4.98280720e+01,
1.08235291e+01, 2.18713860e+01, 4.98280720e+01, 1.62352936e+01, 2.18713860e+01,
4.98280720e+01, 2.16470582e+01, 2.18713860e+01, 4.98280720e+01, 2.70588227e+01,
2.18713860e+01, 4.98280720e+01, 3.24705872e+01, 2.18713860e+01, 4.98280720e+01,
0.00000000e+00, 0.00000000e+00, 5.42467592e+01, 5.41176453e+00, 0.00000000e+00,
5.42467592e+01, 1.08235291e+01, 0.00000000e+00, 5.42467592e+01, 1.62352936e+01,
0.00000000e+00, 5.42467592e+01, 2.16470582e+01, 0.00000000e+00, 5.42467592e+01,
2.70588227e+00, 4.68672556e+00, 5.42467592e+01, 8.11764682e+00, 4.68672556e+00,
5.42467592e+01, 1.35294114e+01, 4.68672556e+00, 5.42467592e+01, 1.89411759e+01,
4.68672556e+00, 5.42467592e+01, 2.43529404e+01, 4.68672556e+00, 5.42467592e+01,
5.41176453e+00, 9.37345114e+00, 5.42467592e+01, 1.08235291e+01, 9.37345114e+00,
5.42467592e+01, 1.62352936e+01, 9.37345114e+00, 5.42467592e+01, 2.16470582e+01,
9.37345114e+00, 5.42467592e+01, 2.70588227e+01, 9.37345114e+00, 5.42467592e+01,
8.11764682e+00, 1.40601767e+01, 5.42467592e+01, 1.35294114e+01, 1.40601767e+01,
5.42467592e+01, 1.89411759e+01, 1.40601767e+01, 5.42467592e+01, 2.43529404e+01,
1.40601767e+01, 5.42467592e+01, 2.97647050e+01, 1.40601767e+01, 5.42467592e+01,
1.08235291e+01, 1.87469023e+01, 5.42467592e+01, 1.62352936e+01, 1.87469023e+01,
5.42467592e+01, 2.16470582e+01, 1.87469023e+01, 5.42467592e+01, 2.70588227e+01,
1.87469023e+01, 5.42467592e+01, 3.24705872e+01, 1.87469023e+01, 5.42467592e+01,
0.00000000e+00, 0.00000000e+00, 5.80262115e+01, -3.02356180e+00, 0.00000000e+00,
5.80262115e+01, 9.44863063e+00, 0.00000000e+00, 5.80262115e+01, 6.42506883e+00,
0.00000000e+00, 5.80262115e+01, 1.88972613e+01, 0.00000000e+00, 5.80262115e+01,
1.58736995e+01, 0.00000000e+00, 5.80262115e+01, 9.44863063e+00, 9.44863063e+00,
5.80262115e+01, 6.42506883e+00, 9.44863063e+00, 5.80262115e+01, 1.88972613e+01,
9.44863063e+00, 5.80262115e+01, 1.58736995e+01, 9.44863063e+00, 5.80262115e+01,
2.83458919e+01, 9.44863063e+00, 5.80262115e+01, 2.53223301e+01, 9.44863063e+00,
5.80262115e+01, 1.88972613e+01, 1.88972613e+01, 5.80262115e+01, 1.58736995e+01,
1.88972613e+01, 5.80262115e+01, 2.83458919e+01, 1.88972613e+01, 5.80262115e+01,
2.53223301e+01, 1.88972613e+01, 5.80262115e+01 ]
</q>
<p shape='(1, 723)'>
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00 ]
</p>
<m shape='(241)'>
[ 4.43053050e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 5.11967350e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.43053050e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 5.11967350e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 5.11967350e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.43053050e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 5.11967350e+04, 4.91843353e+04, 4.43053050e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 5.11967350e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.43053050e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.43053050e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
5.11967350e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
5.11967350e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.43053050e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.43053050e+04, 5.11967350e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 5.11967350e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.43053050e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04,
4.43053050e+04, 4.91843353e+04, 4.91843353e+04, 4.91843353e+04, 5.11967350e+04,
2.91651223e+04, 2.91651223e+04, 2.91651223e+04, 2.91651223e+04, 2.91651223e+04,
2.91651223e+04, 2.91651223e+04, 2.91651223e+04, 2.91651223e+04, 2.91651223e+04,
2.91651223e+04, 2.91651223e+04, 2.91651223e+04, 2.91651223e+04, 2.91651223e+04,
2.91651223e+04 ]
</m>
<names shape='(241)'>
[ Mg, Al, Al, Al, Si,
Al, Al, Al, Al, Al,
Mg, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Si, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Si,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Mg, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Si, Al, Mg, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Si,
Al, Al, Al, Al, Al,
Al, Al, Al, Mg, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Mg, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Si, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Si, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Mg, Al, Al,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Mg, Si, Al, Al, Al,
Al, Al, Al, Al, Si,
Al, Al, Al, Al, Al,
Al, Al, Al, Al, Al,
Al, Mg, Al, Al, Al,
Al, Al, Al, Al, Al,
Mg, Al, Al, Al, Si,
O, O, O, O, O,
O, O, O, O, O,
O, O, O, O, O,
O ]
</names>
</beads>
<cell shape='(3, 3)'>
[ 2.70588227e+01, 1.35294114e+01, 0.00000000e+00, 0.00000000e+00, 2.34336279e+01,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 7.31440205e+01 ]
</cell>
<initialize nbeads='1'><velocities mode='thermal' units='ase'> 800 </velocities></initialize><ensemble><temperature units='ase'> 800 </temperature></ensemble>
<forces>
<force forcefield='pet-mad'> </force>
</forces>
<motion mode="multi">
<motion mode="dynamics">
<dynamics mode="nvt">
<timestep units="ase"> 0.4911347394232032 </timestep>
<thermostat mode='svr'>
<tau units='ase'> 4.911347394232032 </tau>
</thermostat>
</dynamics>
</motion>
<motion mode="atomswap">
<atomswap>
<names> [ Al, Si, Mg, O] </names>
</atomswap>
</motion>
</motion>
</system>
</simulation>
The simulation can be run from a Python script or the command line. By changing the
forcefield interface from direct to the use of a socket, it is also possible to
execute separately i-PI and the metatensor driver.
sim = InteractiveSimulation(input_xml)
sim.run(100)
@system: Initializing system object
@simulation: Initializing simulation object
@ RANDOM SEED: The seed used in this calculation was 1745487473541
@initializer: Resampling velocities at temperature 800.0 ase
@system.bind: Binding the forces
0
@simulation.run: Average timings at MD step 0. t/step: 9.04006e-01
The simulation generates output files that can be parsed and visualized from Python.
data, info = read_output("nvt_atomxc.out")
trj = read_trajectory("nvt_atomxc.pos_0.xyz")
fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True)
ax.plot(data["time"], data["potential"], "b-", label="potential")
ax.plot(data["time"], data["conserved"] - 4, "k-", label="conserved")
ax.set_xlabel("t / ps")
ax.set_ylabel("energy / ev")
ax.legend()
<matplotlib.legend.Legend object at 0x7ff79283e750>
The trajectory (which is started from oxygen molecules placed on top of the surface) shows quick relaxation to an oxide layer. If you look carefully, you’ll also see that Mg and Si atoms tend to cluster together, and accumulate at the surface.
chemiscope.show(
frames=trj,
properties={
"time": data["time"][::10],
"potential": data["potential"][::10],
"temperature": data["temperature"][::10],
},
mode="default",
settings=chemiscope.quick_settings(
map_settings={
"x": {"property": "time", "scale": "linear"},
"y": {"property": "potential", "scale": "linear"},
},
structure_settings={
"unitCell": True,
},
trajectory=True,
),
)
Molecular dynamics with LAMMPS¶
We now run the same MD with LAMMPS. To run a LAMMPS
calculation with a metatomic potential, one needs a LAMMPS build that contains an
appropriate pair style. You can compile it from source, or fetch it from the metatensor channel on
conda. One can then just include in the input a pair_style metatensor that points
to the exported model and a single pair_coeff command that specifies the mapping
from LAMMPS types to the atomic types the model can handle. The first two arguments
must be * * so as to span all LAMMPS atom types. This is followed by a list of N
arguments that specify the mapping of metatensor atomic types to LAMMPS types, where N
is the number of LAMMPS atom types.
with open("data/al6xxx-o2.in", "r") as f:
print(f.read())
units metal # Angstroms, eV, picoseconds
atom_style atomic
read_data al6xxx-o2.data
# loads pet-mad-model
pair_style metatensor &
pet-mad-latest.pt &
device cpu &
extensions extensions/
# define interactions between all atoms and maps the LAMMPS types to elements
pair_coeff * * 13 12 8 14
neighbor 2.0 bin
timestep 0.005
dump myDump all xyz 10 trajectory.xyz
dump_modify myDump element Al Mg O Si
thermo_style multi
thermo 1
velocity all create 800 87287 mom yes rot yes
fix 1 all nvt temp 800 800 0.10
# fix 2 all atom/swap 1 1 12345 800 types 1 2
# fix 2 all atom/swap 1 1 12345 800 types 1 3
# fix 2 all atom/swap 1 1 12345 800 types 1 4
run 100
Warning
Be aware that the extensions are compiled files and depend on your operating system.
Usually you have re-export the extensions for different systems! You can do this
by running the appropriate parts of this file, or using the mtt export
command-line utility.
We also save the geometry to a LAMMPS data file and finally run the simulation.
ase.io.write("al6xxx-o2.data", al_surface, format="lammps-data", masses=True)
subprocess.check_call(["lmp", "-in", "data/al6xxx-o2.in"])
0
The resulting trajectory is qualitatively consistent with what
we observed with i-PI.
lmp_trj = ase.io.read("trajectory.xyz", ":")
chemiscope.show(frames=lmp_trj, mode="structure")
Total running time of the script: (5 minutes 3.935 seconds)