Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement lattice transform to upper-triangular matrix for ASE Nose-Hoover thermostat #68

Merged
merged 9 commits into from
Sep 12, 2023
Merged

Implement lattice transform to upper-triangular matrix for ASE Nose-Hoover thermostat #68

Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
0db053c
Implemented transformation of lattice to upper-triangular matrix
Sep 11, 2023
dfd8bc4
Added verbose flag and explanation to upper_triangular_cell
Sep 11, 2023
1beaf1b
FIXED: ASE NPT does not accept externalstress=None
Sep 12, 2023
f22adc1
Added tests for Nose-Hoover NVT and NPT MD
Sep 12, 2023
a6480fd
Changed MD test numerical tolerance
Sep 12, 2023
a976874
use numpy.testing.assert_allclose over np.isclose().all() for better …
janosh Sep 12, 2023
cfba1c9
tweak upper_triangular_cell() doc str
janosh Sep 12, 2023
055b8ed
Merge branch 'main' into upper_triangular_cell
janosh Sep 12, 2023
948a8bd
fix black
janosh Sep 12, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 32 additions & 3 deletions chgnet/model/dynamics.py
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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")
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is there a need to include LCO for Nose-Hoover test?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is there a need to include LCO for Nose-Hoover test?

Because the cell matrix of LiMnO2 is diagonal.

Copy link
Collaborator

@bowen-bd bowen-bd Sep 12, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This makes sense. Another way is to use rotation and non-diagonal supercell expansion to create a non-diagonal LiMnO2 cell

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=2.1e-3, atol=1e-8).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=2.1e-3, atol=1e-8).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