Skip to content

Commit

Permalink
Implemented transformation of lattice to upper-triangular matrix (#68)
Browse files Browse the repository at this point in the history 
* Implemented transformation of lattice to upper-triangular matrix

* Added verbose flag and explanation to upper_triangular_cell

* FIXED: ASE NPT does not accept externalstress=None

* Added tests for Nose-Hoover NVT and NPT MD

* Changed MD test numerical tolerance

* use numpy.testing.assert_allclose over np.isclose().all() for better error messages on failing tests

* tweak upper_triangular_cell() doc str

---------

Co-authored-by: Ziyang HU <[email protected]>
Co-authored-by: Janosh Riebesell <[email protected]>
  • Loading branch information
@tsihyoung @janosh
3 people authored Sep 12, 2023
1 parent 639fcd8 commit 1d1c0af
Show file tree
Hide file tree
Showing 3 changed files with 215 additions and 15 deletions.
35 changes: 32 additions & 3 deletions chgnet/model/dynamics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -479,12 +479,14 @@ def __init__(
Nose-hoover (constant N, V, T) molecular dynamics.
ASE implementation currently only supports upper triangular lattice
"""
self.upper_triangular_cell()
self.dyn = NPT(
atoms=self.atoms,
timestep=timestep * units.fs,
temperature_K=temperature,
externalstress=pressure
* units.GPa, # ase NPT does not like externalstress=None
ttime=taut * units.fs,
externalstress=None,
pfactor=None,
trajectory=trajectory,
logfile=logfile,
Expand Down Expand Up @@ -551,6 +553,7 @@ def __init__(
see: https://gitlab.com/ase/ase/-/blob/master/ase/md/npt.py
ASE implementation currently only supports upper triangular lattice
"""
self.upper_triangular_cell()
ptime = taup * units.fs
self.dyn = NPT(
atoms=self.atoms,
Expand Down Expand Up @@ -628,7 +631,6 @@ def run(self, steps: int):
Args:
steps (int): number of MD steps
Returns:
"""
if self.crystal_feas_logfile:
obs = CrystalFeasObserver(self.atoms)
Expand All @@ -644,13 +646,40 @@ def set_atoms(self, atoms: Atoms):
Args:
atoms (Atoms): new atoms for running MD
Returns:
"""
calculator = self.atoms.calc
self.atoms = atoms
self.dyn.atoms = atoms
self.dyn.atoms.calc = calculator

def upper_triangular_cell(self, verbose: bool | None = False):
"""Transform to upper-triangular cell.
ASE Nose-Hoover implementation only supports upper-triangular cell
while ASE's canonical description is lower-triangular cell.
Args:
verbose (bool): Whether to notify user about upper-triangular cell transformation.
Default = False
"""
if not NPT._isuppertriangular(self.atoms.get_cell()):
a, b, c, alpha, beta, gamma = self.atoms.cell.cellpar()
angles = np.radians((alpha, beta, gamma))
sin_a, sin_b, _sin_g = np.sin(angles)
cos_a, cos_b, cos_g = np.cos(angles)
cos_p = (cos_g - cos_a * cos_b) / (sin_a * sin_b)
cos_p = np.clip(cos_p, -1, 1)
sin_p = (1 - cos_p**2) ** 0.5

new_basis = [
(a * sin_b * sin_p, a * sin_b * cos_p, a * cos_b),
(0, b * sin_a, b * cos_a),
(0, 0, c),
]

self.atoms.set_cell(new_basis, scale_atoms=True)
if verbose:
print("Transformed to upper triangular unit cell.", flush=True)


class EquationOfState:
"""Class to calculate equation of state."""
Expand Down
65 changes: 65 additions & 0 deletions examples/mp-1175469-Li9Co7O16.cif
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,65 @@
# generated using pymatgen
data_Li9Co7O16
_symmetry_space_group_name_H-M 'P 1'
_cell_length_a 6.30240647
_cell_length_b 6.30741300
_cell_length_c 7.41583051
_cell_angle_alpha 75.47000443
_cell_angle_beta 65.23135009
_cell_angle_gamma 72.00717747
_symmetry_Int_Tables_number 1
_chemical_formula_structural Li9Co7O16
_chemical_formula_sum 'Li9 Co7 O16'
_cell_volume 252.05269201
_cell_formula_units_Z 1
loop_
_symmetry_equiv_pos_site_id
_symmetry_equiv_pos_as_xyz
1 'x, y, z'
loop_
_atom_type_symbol
_atom_type_oxidation_number
Li+ 1.0
Co3+ 3.0
Co4+ 4.0
O2- -2.0
loop_
_atom_site_type_symbol
_atom_site_label
_atom_site_symmetry_multiplicity
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
_atom_site_occupancy
Li+ Li0 1 0.31442518 0.68680206 0.56469026 1
Li+ Li1 1 0.80552981 0.18376271 0.06064773 1
Li+ Li2 1 0.43695568 0.57151353 0.18253456 1
Li+ Li3 1 0.92705143 0.06026526 0.69553727 1
Li+ Li4 1 0.56304432 0.42848647 0.81746544 1
Li+ Li5 1 0.07294857 0.93973474 0.30446273 1
Li+ Li6 1 0.68557482 0.31319794 0.43530974 1
Li+ Li7 1 0.19447019 0.81623729 0.93935227 1
Li+ Li8 1 0.50000000 0.00000000 0.00000000 1
Co3+ Co9 1 0.87189169 0.62805253 0.87615769 1
Co3+ Co10 1 0.12810831 0.37194747 0.12384231 1
Co3+ Co11 1 0.37568054 0.12368956 0.37222619 1
Co3+ Co12 1 0.00000000 0.50000000 0.50000000 1
Co3+ Co13 1 0.62431946 0.87631044 0.62777381 1
Co4+ Co14 1 0.25225503 0.24664332 0.75326548 1
Co4+ Co15 1 0.74774497 0.75335668 0.24673452 1
O2- O16 1 0.20298483 0.51557922 0.85294009 1
O2- O17 1 0.70174803 0.01784432 0.35191128 1
O2- O18 1 0.32496647 0.38970329 0.47992830 1
O2- O19 1 0.84540853 0.88153215 0.97729332 1
O2- O20 1 0.44478863 0.28084762 0.11226145 1
O2- O21 1 0.95254917 0.76784039 0.60312012 1
O2- O22 1 0.57889279 0.14834865 0.71233501 1
O2- O23 1 0.07585520 0.64172837 0.22862318 1
O2- O24 1 0.42110721 0.85165135 0.28766499 1
O2- O25 1 0.92414480 0.35827163 0.77137682 1
O2- O26 1 0.55521137 0.71915238 0.88773855 1
O2- O27 1 0.04745083 0.23215961 0.39687988 1
O2- O28 1 0.67503353 0.61029671 0.52007170 1
O2- O29 1 0.15459147 0.11846785 0.02270668 1
O2- O30 1 0.79701517 0.48442078 0.14705991 1
O2- O31 1 0.29825197 0.98215568 0.64808872 1
130 changes: 118 additions & 12 deletions tests/test_md.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,10 +7,14 @@
import numpy as np
import pytest
from ase import Atoms
from ase.md.npt import NPT
from ase.md.nptberendsen import Inhomogeneous_NPTBerendsen
from ase.md.nvtberendsen import NVTBerendsen
from ase.md.verlet import VelocityVerlet
from numpy.testing import assert_allclose
from pymatgen.analysis.structure_matcher import StructureMatcher
from pymatgen.core import Structure
from pymatgen.io.ase import AseAtomsAdaptor
from pytest import MonkeyPatch, approx

from chgnet import ROOT
Expand All @@ -24,6 +28,7 @@

relaxer = StructOptimizer()
structure = Structure.from_file(f"{ROOT}/examples/mp-18767-LiMnO2.cif")
structure_2 = Structure.from_file(f"{ROOT}/examples/mp-1175469-Li9Co7O16.cif")
chgnet = CHGNet.load()


Expand Down Expand Up @@ -66,12 +71,11 @@ def test_md_nvt(
assert isinstance(md.atoms, Atoms)
assert isinstance(md.atoms.calc, CHGNetCalculator)
assert isinstance(md.dyn, NVTBerendsen)
assert os.path.isfile("md_out.traj")
assert os.path.isfile("md_out.log")
assert set(os.listdir()) == {"md_out.log", "md_out.traj"}
with open("md_out.log") as log_file:
next(log_file)
logs = log_file.read()
logs = np.fromstring(logs, dtype=float, sep=" ")
logs = np.fromstring(logs, sep=" ")
ref = np.fromstring(
"0.0000 -58.9727 -58.9727 0.0000 0.0\n"
"0.0200 -58.9723 -58.9731 0.0009 0.8\n"
Expand All @@ -84,10 +88,9 @@ def test_md_nvt(
"0.1600 -58.4724 -58.5531 0.0807 78.1\n"
"0.1800 -58.3891 -58.8077 0.4186 404.8\n"
"0.2000 -58.3398 -58.9244 0.5846 565.4\n",
dtype=float,
sep=" ",
)
assert np.isclose(logs, ref, rtol=1e-2, atol=1e-7).all()
assert_allclose(logs, ref, rtol=2.1e-3, atol=1e-8)


def test_md_nve(tmp_path: Path, monkeypatch: MonkeyPatch):
Expand All @@ -107,8 +110,7 @@ def test_md_nve(tmp_path: Path, monkeypatch: MonkeyPatch):
assert isinstance(md.atoms, Atoms)
assert isinstance(md.atoms.calc, CHGNetCalculator)
assert isinstance(md.dyn, VelocityVerlet)
assert os.path.isfile("md_out.traj")
assert os.path.isfile("md_out.log")
assert set(os.listdir()) == {"md_out.log", "md_out.traj"}
with open("md_out.log") as log_file:
logs = log_file.read()
assert logs == (
Expand Down Expand Up @@ -139,12 +141,11 @@ def test_md_npt_inhomogeneous_berendsen(tmp_path: Path, monkeypatch: MonkeyPatch
assert isinstance(md.dyn, Inhomogeneous_NPTBerendsen)
assert md.bulk_modulus == approx(105.764, rel=1e-2)
assert md.dyn.pressure == approx(6.324e-07, rel=1e-4)
assert os.path.isfile("md_out.traj")
assert os.path.isfile("md_out.log")
assert set(os.listdir()) == {"md_out.log", "md_out.traj"}
with open("md_out.log") as log_file:
next(log_file)
logs = log_file.read()
logs = np.fromstring(logs, dtype=float, sep=" ")
logs = np.fromstring(logs, sep=" ")
ref = np.fromstring(
"0.0000 -58.9727 -58.9727 0.0000 0.0\n"
"0.0200 -58.9723 -58.9731 0.0009 0.8\n"
Expand All @@ -157,10 +158,115 @@ def test_md_npt_inhomogeneous_berendsen(tmp_path: Path, monkeypatch: MonkeyPatch
"0.1600 -58.4731 -58.5533 0.0802 77.6\n"
"0.1800 -58.3897 -58.8064 0.4167 402.9\n"
"0.2000 -58.3404 -58.9253 0.5849 565.6\n",
dtype=float,
sep=" ",
)
assert np.isclose(logs, ref, rtol=1e-2, atol=1e-7).all()
assert_allclose(logs, ref, rtol=2.1e-3, atol=1e-8)


def test_md_nvt_nose_hoover(tmp_path: Path, monkeypatch: MonkeyPatch):
monkeypatch.chdir(tmp_path) # run MD in temporary directory

md = MolecularDynamics(
atoms=structure_2,
model=chgnet,
ensemble="nvt",
thermostat="nose-hoover",
temperature=1000, # in k
timestep=2, # in fs
trajectory="md_out.traj",
logfile="md_out.log",
loginterval=10,
)

new_atoms = AseAtomsAdaptor.get_structure(md.atoms)
assert StructureMatcher(
primitive_cell=False, scale=False, attempt_supercell=False, allow_subset=False
).fit(structure_2, new_atoms)

md.run(100)

assert isinstance(md.atoms, Atoms)
assert isinstance(md.atoms.calc, CHGNetCalculator)
assert isinstance(md.dyn, NPT)
assert_allclose(
md.dyn.externalstress,
[-6.324e-07, -6.324e-07, -6.324e-07, 0.0, 0.0, 0.0],
rtol=1e-4,
atol=1e-8,
)
assert set(os.listdir()) == {"md_out.log", "md_out.traj"}
with open("md_out.log") as log_file:
next(log_file)
logs = log_file.read()
logs = np.fromstring(logs, sep=" ")
ref = np.fromstring(
"0.0200 -199.2479 -199.3994 0.1515 36.6\n"
"0.0400 -199.2459 -199.3440 0.0981 23.7\n"
"0.0600 -199.2394 -199.2669 0.0275 6.6\n"
"0.0800 -199.2348 -199.4143 0.1795 43.4\n"
"0.1000 -199.2274 -199.2774 0.0500 12.1\n"
"0.1200 -199.2123 -199.3001 0.0878 21.2\n"
"0.1400 -199.2040 -199.4000 0.1961 47.4\n"
"0.1600 -199.1856 -199.2181 0.0325 7.9\n"
"0.1800 -199.1603 -199.3266 0.1662 40.2\n"
"0.2000 -199.1455 -199.3490 0.2035 49.2\n",
sep=" ",
)
assert_allclose(logs, ref, rtol=1e-2, atol=1e-7)


def test_md_npt_nose_hoover(tmp_path: Path, monkeypatch: MonkeyPatch):
# https://github.com/CederGroupHub/chgnet/pull/68
monkeypatch.chdir(tmp_path) # run MD in temporary directory

md = MolecularDynamics(
atoms=structure_2,
model=chgnet,
ensemble="npt",
thermostat="nose-hoover",
temperature=1000, # in k
timestep=2, # in fs
trajectory="md_out.traj",
logfile="md_out.log",
loginterval=10,
)

new_atoms = AseAtomsAdaptor.get_structure(md.atoms)
assert StructureMatcher(
primitive_cell=False, scale=False, attempt_supercell=False, allow_subset=False
).fit(structure_2, new_atoms)

md.run(100)

assert isinstance(md.atoms, Atoms)
assert isinstance(md.atoms.calc, CHGNetCalculator)
assert isinstance(md.dyn, NPT)
assert md.bulk_modulus == approx(102.977, rel=1e-2)
assert_allclose(
md.dyn.externalstress,
[-6.324e-07, -6.324e-07, -6.324e-07, 0.0, 0.0, 0.0],
rtol=1e-4,
atol=1e-8,
)
assert set(os.listdir()) == {"md_out.log", "md_out.traj"}
with open("md_out.log") as log_file:
next(log_file)
logs = log_file.read()
logs = np.fromstring(logs, sep=" ")
ref = np.fromstring(
"0.0200 -199.2480 -199.3994 0.1514 36.6\n"
"0.0400 -199.2460 -199.3442 0.0982 23.7\n"
"0.0600 -199.2397 -199.2672 0.0275 6.7\n"
"0.0800 -199.2355 -199.4148 0.1793 43.3\n"
"0.1000 -199.2282 -199.2782 0.0500 12.1\n"
"0.1200 -199.2135 -199.3017 0.0882 21.3\n"
"0.1400 -199.2060 -199.4014 0.1954 47.2\n"
"0.1600 -199.1878 -199.2201 0.0323 7.8\n"
"0.1800 -199.1630 -199.3306 0.1675 40.5\n"
"0.2000 -199.1496 -199.3506 0.2010 48.6\n",
sep=" ",
)
assert_allclose(logs, ref, rtol=1e-2, atol=1e-7)


def test_md_crystal_feas_log(
Expand Down

0 comments on commit 1d1c0af

Please sign in to comment.