Periodic Boundaries not recognized during molecular dynamics simulation

[VerySpeedyPhoton]

I would like to use MatGL on a system with periodic boundary conditions, but have been unable to figure out how to do this. The issue is that even when periodic boundaries are specified in the system definitions, MatGL doesn't seem to recognize them.

To illustrate the issue, here are two python scripts that simulate 4 atoms of copper at 2400K. One uses the standard ASE package (with the EMT potential); one uses MatGL.

(1) ASE script: (produces a .traj file)

from asap3 import EMT  
from ase import units
from ase.io.trajectory import Trajectory
from ase.lattice.cubic import FaceCenteredCubic
from ase.md.langevin import Langevin
from pymatgen.core import Lattice, Structure
from pymatgen.io.ase import AseAtomsAdaptor
import warnings
warnings.simplefilter("ignore")

# Define the lattice geometry and atom coordinates
lattice = Lattice.cubic(3.61, (True, True, True))       # lattice constant in units of Angstrom ; (True, True, True) sets periodic boundaries
fractcoords = [[0, 0, 0], [0, 0.5, 0.5], [0.5, 0, 0.5], [0.5, 0.5, 0]]
struct = Structure(lattice, ["Cu", "Cu", "Cu", "Cu"], fractcoords)

ase_adaptor = AseAtomsAdaptor()
atoms = ase_adaptor.get_atoms(struct)

# Describe the interatomic interactions with the Effective Medium Theory
atoms.calc = EMT()

T = 2400  # Kelvin
dyn = Langevin(atoms, 5 * units.fs, T * units.kB, 0.002)

def printenergy(a=atoms):  # store a reference to atoms in the definition.
    """Function to print the potential, kinetic and total energy."""
    epot = a.get_potential_energy() / len(a)
    ekin = a.get_kinetic_energy() / len(a)
    print('Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)  '
          'Etot = %.3feV' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))

dyn.attach(printenergy, interval=50)
traj = Trajectory('Cu_ASE.traj', 'w', atoms)
dyn.attach(traj.write, interval=50)
# Now run the dynamics
printenergy()
dyn.run(2000)

(2) MatGL script: (produces .traj and .log files)

from __future__ import annotations
import warnings
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from pymatgen.core import Lattice, Structure
from pymatgen.io.ase import AseAtomsAdaptor
import matgl
from matgl.ext.ase import M3GNetCalculator, MolecularDynamics, Relaxer

warnings.simplefilter("ignore")
pot = matgl.load_model("M3GNet-MP-2021.2.8-PES")

# Define the lattice geometry and atom coordinates
lattice = Lattice.cubic(3.61, (True, True, True))       # lattice constant in units of Angstrom ; (True, True, True) sets periodic boundaries
fractcoords = [[0, 0, 0], [0, 0.5, 0.5], [0.5, 0, 0.5], [0.5, 0.5, 0]]
struct = Structure(lattice, ["Cu", "Cu", "Cu", "Cu"], fractcoords)

# Prepare atoms for molecular dynamics
ase_adaptor = AseAtomsAdaptor()
atoms = ase_adaptor.get_atoms(struct)

# Initiate temperature distribution
MaxwellBoltzmannDistribution(atoms, temperature_K=2400)

# Define molecular dynamics settings
driver = MolecularDynamics(
        atoms,
        potential=pot,                # uses the M3GNet interatomic potential
        temperature=2400,
        timestep=1, # 1fs,
        logfile="Cu_MatGL.log",
        trajectory="Cu_MatGL.traj",       # save trajectory
        loginterval=50,               # interval for recording the log
        ensemble='nvt'                # NVT ensemble
)

# Run molecular dynamics
driver.run(2000)

Also, to convert one of the .traj files to human-readable format, here is a little script: (produces .xyz trajectory file)

import ase
from ase.io import read, write
from ase.io.trajectory import Trajectory

traj = Trajectory("Cu_MatGL.traj")
#traj = Trajectory("Cu_ASE.traj")

atoms=traj[:]

writeme = ase.io.write("Cu_MatGL.xyz", atoms, "xyz")
#writeme = ase.io.write("Cu_ASE.xyz", atoms, "xyz")

Looking at the .xyz files, you'll notice that in the basic ASE simulation, none of the atom positions are larger than the specified lattice constant (3.61 Angstroms). This is as expected with periodic boundary conditions. However, in the MatGL simulation, the atom positions readily exceed the lattice constant value, and the atoms appear to drift through space.

For example, here are excerpts of the last recorded frame from the .xyz trajectory files: (1) ASE trajectory excerpt

4

Cu      0.344162968499945     -0.148546966782993     -0.009870944387788
Cu     -0.177349810618174      1.722229170919916      1.765656855202022
Cu      1.736425681826978      0.345596431729686      1.793462024164234
Cu      1.727657879060445      1.651865872484211      0.104391311992427

(2) MatGL trajectory excerpt

4

Cu     10.754732275176222      1.779757546738823      7.457469460754528
Cu     10.395401224745564      3.745542174112556      9.486116801335781
Cu     12.303262940663885      1.970363284378915      9.334171650233468
Cu     12.388181989458412      3.671734565914415      7.549590636835549

Do you have any recommendations for how to resolve this issue? How should one go about implementing periodic boundary conditions in a molecular dynamics simulation that leverages MatGL? Please advise, thank you.

[shyuep]

It is not true that MatGL does not support PbC. In all instances, the potentials were trained on crystals with PBC and that is the predominant use case. I am more puzzled why you'd think having atoms staying within the lattice parameter is a sign of PBC. By definition, periodic boundary means that a coordinate of (1000, 1000, 1000) is equivalent to (0,0,0) if the lattice constant is 1.

[kenko911]

Hi @VerySpeedyPhoton, thank you very much for your questions. If you have already had a look at the log file and visualized the MD trajectory using ase-gui. There is no problem for MatGL to recognize the periodic systems. In other words, the way you convert ASE trajectory into other formats cause the problem. I suggest you read the ASE documentation carefully to understand how to convert the ASE MD trajectory into some formats, which can visualize periodic systems properly. In my opinion, xyz format is not appropriate for visualizing periodic systems since they don't support periodic boundary condition and you should try other formats instead like extended xyz or vasp format.