Skip to content

Commit

Permalink
Merge pull request #903 from openforcefield/use-lammps-api
Browse files Browse the repository at this point in the history
Use LAMMPS's Python API in energy driver
  • Loading branch information
mattwthompson authored Feb 15, 2024
2 parents 37162e2 + 8adc4c0 commit e7f0298
Show file tree
Hide file tree
Showing 3 changed files with 32 additions and 61 deletions.
3 changes: 1 addition & 2 deletions openff/interchange/_tests/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@
from openff.utilities.utilities import has_executable, has_package

from openff.interchange.drivers.gromacs import _find_gromacs_executable
from openff.interchange.drivers.lammps import _find_lammps_executable

if sys.version_info >= (3, 10):
from importlib import resources
Expand Down Expand Up @@ -96,7 +95,7 @@ def from_mapped_smiles(self, smiles, name="", **kwargs):


HAS_GROMACS = _find_gromacs_executable() is not None
HAS_LAMMPS = _find_lammps_executable() is not None
HAS_LAMMPS = has_package("lammps")
HAS_SANDER = has_executable("sander")

needs_gmx = pytest.mark.skipif(not HAS_GROMACS, reason="Needs GROMACS")
Expand Down
18 changes: 9 additions & 9 deletions openff/interchange/_tests/unit_tests/drivers/test_all.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@
"""

import math
from shutil import which

import pandas
import pytest
Expand All @@ -12,6 +11,9 @@

from openff.interchange import Interchange
from openff.interchange._tests import (
HAS_GROMACS,
HAS_LAMMPS,
HAS_SANDER,
needs_gmx,
needs_lmp,
needs_not_gmx,
Expand All @@ -20,8 +22,6 @@
needs_sander,
)
from openff.interchange.drivers.all import get_all_energies, get_summary_data
from openff.interchange.drivers.gromacs import _find_gromacs_executable
from openff.interchange.drivers.lammps import _find_lammps_executable


@skip_if_missing("openmm")
Expand Down Expand Up @@ -59,9 +59,9 @@ def test_all_with_minimum(self, basic_interchange):
def test_skipping(self, basic_interchange):
summary = get_all_energies(basic_interchange)

assert ("GROMACS" in summary) == (_find_gromacs_executable() is not None)
assert ("Amber" in summary) == (which("sander") is not None)
assert ("LAMMPS" in summary) == (_find_lammps_executable() is not None)
assert ("GROMACS" in summary) == HAS_GROMACS
assert ("Amber" in summary) == HAS_SANDER
assert ("LAMMPS" in summary) == HAS_LAMMPS

# TODO: Also run all of this with h-bond constraints
def test_summary_data(self, basic_interchange):
Expand All @@ -71,9 +71,9 @@ def test_summary_data(self, basic_interchange):

assert "OpenMM" in summary.index

assert ("GROMACS" in summary.index) == (_find_gromacs_executable() is not None)
assert ("Amber" in summary.index) == (which("sander") is not None)
assert ("LAMMPS" in summary.index) == (_find_lammps_executable() is not None)
assert ("GROMACS" in summary.index) == HAS_GROMACS
assert ("Amber" in summary.index) == HAS_SANDER
assert ("LAMMPS" in summary.index) == HAS_LAMMPS

# Check that (some of) the data is reasonable, this tolerance should be greatly reduced
# See https:/openforcefield/openff-interchange/issues/632
Expand Down
72 changes: 22 additions & 50 deletions openff/interchange/drivers/lammps.py
Original file line number Diff line number Diff line change
@@ -1,32 +1,17 @@
"""Functions for running energy evluations with LAMMPS."""

import subprocess
from shutil import which
from typing import Optional

import numpy
from openff.units import unit
from openff.units import Quantity, unit
from openff.utilities import MissingOptionalDependencyError, requires_package

from openff.interchange import Interchange
from openff.interchange.components.mdconfig import MDConfig
from openff.interchange.drivers.report import EnergyReport
from openff.interchange.exceptions import LAMMPSNotFoundError, LAMMPSRunError


def _find_lammps_executable(raise_exception: bool = False) -> Optional[str]:
"""Attempt to locate a LAMMPS executable based on commonly-used names."""
lammps_executable_names = ["lammps", "lmp_serial", "lmp_mpi"]

for name in lammps_executable_names:
if which(name):
return name

if raise_exception:
raise LAMMPSNotFoundError
else:
return None


def get_lammps_energies(
interchange: Interchange,
round_positions: Optional[int] = None,
Expand Down Expand Up @@ -55,17 +40,21 @@ def get_lammps_energies(
An `EnergyReport` object containing the single-point energies.
"""
return _process(
_get_lammps_energies(interchange, round_positions),
detailed,
)
try:
return _process(
_get_lammps_energies(interchange, round_positions),
detailed,
)
except MissingOptionalDependencyError:
raise LAMMPSNotFoundError


@requires_package("lammps")
def _get_lammps_energies(
interchange: Interchange,
round_positions: Optional[int] = None,
) -> dict[str, unit.Quantity]:
lmp = _find_lammps_executable(raise_exception=True)
import lammps

if round_positions is not None:
interchange.positions = numpy.round(interchange.positions, round_positions)
Expand All @@ -76,24 +65,21 @@ def _get_lammps_energies(
input_file="tmp.in",
)

run_cmd = f"{lmp} -i tmp.in"
runner = lammps.lammps()

proc = subprocess.Popen(
run_cmd,
shell=True,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
universal_newlines=True,
)

_, err = proc.communicate()

if proc.returncode:
raise LAMMPSRunError(err)
try:
runner.file("tmp.in")
# LAMMPS does not raise a custom exception :(
except Exception as error:
raise LAMMPSRunError from error

# thermo_style custom ebond eangle edihed eimp epair evdwl ecoul elong etail pe
parsed_energies = unit.kilocalorie_per_mole * _parse_lammps_log("log.lammps")
parsed_energies = [
Quantity(energy, "kilocalorie_per_mole")
for energy in runner.last_thermo().values()
]

# TODO: Sanely map LAMMPS's energy names to the ones we care about
return {
"Bond": parsed_energies[0],
"Angle": parsed_energies[1],
Expand Down Expand Up @@ -124,17 +110,3 @@ def _process(
),
},
)


def _parse_lammps_log(file_in: str) -> list[float]:
"""Parse a LAMMPS log file for energy components."""
tag = False
with open(file_in) as fi:
for line in fi.readlines():
if tag:
data = [float(val) for val in line.split()]
tag = False
if line.strip().startswith("E_bond"):
tag = True

return data

0 comments on commit e7f0298

Please sign in to comment.