Coverage for /builds/ase/ase/ase/calculators/lammpsrun.py: 87.68%
211 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3"""ASE calculator for the LAMMPS classical MD code"""
4# lammps.py (2011/03/29)
5#
6# Copyright (C) 2009 - 2011 Joerg Meyer, joerg.meyer@ch.tum.de
7#
8# This library is free software; you can redistribute it and/or
9# modify it under the terms of the GNU Lesser General Public
10# License as published by the Free Software Foundation; either
11# version 2.1 of the License, or (at your option) any later version.
12#
13# This library is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16# Lesser General Public License for more details.
17#
18# You should have received a copy of the GNU Lesser General Public
19# License along with this file; if not, write to the Free Software
20# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301
21# USA or see <http://www.gnu.org/licenses/>.
24import os
25import shlex
26import shutil
27import subprocess
28import warnings
29from re import IGNORECASE
30from re import compile as re_compile
31from tempfile import NamedTemporaryFile, mkdtemp
32from tempfile import mktemp as uns_mktemp
33from threading import Thread
34from typing import Any, Dict
36import numpy as np
38from ase.calculators.calculator import Calculator, all_changes
39from ase.calculators.lammps import (
40 CALCULATION_END_MARK,
41 Prism,
42 convert,
43 write_lammps_in,
44)
45from ase.data import atomic_masses, chemical_symbols
46from ase.io.lammpsdata import write_lammps_data
47from ase.io.lammpsrun import read_lammps_dump
49__all__ = ["LAMMPS"]
52class LAMMPS(Calculator):
53 """LAMMPS (https://lammps.sandia.gov/) calculator
55 The environment variable :envvar:`ASE_LAMMPSRUN_COMMAND` must be defined to
56 tell ASE how to call a LAMMPS binary. This should contain the path to the
57 lammps binary, or more generally, a command line possibly also including an
58 MPI-launcher command.
60 For example (in a Bourne-shell compatible environment):
62 .. code-block:: bash
64 export ASE_LAMMPSRUN_COMMAND=/path/to/lmp_binary
66 or possibly something similar to
68 .. code-block:: bash
70 export ASE_LAMMPSRUN_COMMAND="/path/to/mpirun --np 4 lmp_binary"
72 Parameters
73 ----------
74 files : list[str]
75 List of files needed by LAMMPS. Typically a list of potential files.
76 parameters : dict[str, Any]
77 Dictionary of settings to be passed into the input file for calculation.
78 specorder : list[str]
79 Within LAAMPS, atoms are identified by an integer value starting from 1.
80 This variable allows the user to define the order of the indices
81 assigned to the atoms in the calculation, with the default
82 if not given being alphabetical
83 keep_tmp_files : bool, default: False
84 Retain any temporary files created. Mostly useful for debugging.
85 tmp_dir : str, default: None
86 path/dirname (default None -> create automatically).
87 Explicitly control where the calculator object should create
88 its files. Using this option implies 'keep_tmp_files=True'.
89 no_data_file : bool, default: False
90 Controls whether an explicit data file will be used for feeding
91 atom coordinates into lammps. Enable it to lessen the pressure on
92 the (tmp) file system. THIS OPTION MIGHT BE UNRELIABLE FOR CERTAIN
93 CORNER CASES (however, if it fails, you will notice...).
94 keep_alive : bool, default: True
95 When using LAMMPS as a spawned subprocess, keep the subprocess
96 alive (but idling when unused) along with the calculator object.
97 always_triclinic : bool, default: False
98 Force LAMMPS to treat the cell as tilted, even if the cell is not
99 tilted, by printing ``xy``, ``xz``, ``yz`` in the data file.
100 reduce_cell : bool, default: False
101 If True, reduce cell to have shorter lattice vectors.
102 write_velocities : bool, default: False
103 If True, forward ASE velocities to LAMMPS.
104 verbose: bool, default: False
105 If True, print additional debugging output to STDOUT.
107 Examples
108 --------
109 Provided that the respective potential file is in the working directory,
110 one can simply run (note that LAMMPS needs to be compiled to work with EAM
111 potentials)
113 ::
115 from ase import Atom, Atoms
116 from ase.build import bulk
117 from ase.calculators.lammpsrun import LAMMPS
119 parameters = {'pair_style': 'eam/alloy',
120 'pair_coeff': ['* * NiAlH_jea.eam.alloy H Ni']}
122 files = ['NiAlH_jea.eam.alloy']
124 Ni = bulk('Ni', cubic=True)
125 H = Atom('H', position=Ni.cell.diagonal()/2)
126 NiH = Ni + H
128 lammps = LAMMPS(files=files, **parameters)
130 NiH.calc = lammps
131 print("Energy ", NiH.get_potential_energy())
132 """
134 name = "lammpsrun"
135 implemented_properties = ["energy", "free_energy", "forces", "stress",
136 "energies"]
138 # parameters to choose options in LAMMPSRUN
139 ase_parameters: Dict[str, Any] = dict(
140 specorder=None,
141 atorder=True,
142 always_triclinic=False,
143 reduce_cell=False,
144 keep_alive=True,
145 keep_tmp_files=False,
146 no_data_file=False,
147 tmp_dir=None,
148 files=[], # usually contains potential parameters
149 verbose=False,
150 write_velocities=False,
151 binary_dump=True, # bool - use binary dump files (full
152 # precision but long long ids are casted to
153 # double)
154 lammps_options="-echo log -screen none -log /dev/stdout",
155 trajectory_out=None, # file object, if is not None the
156 # trajectory will be saved in it
157 )
159 # parameters forwarded to LAMMPS
160 lammps_parameters = dict(
161 boundary=None, # bounadry conditions styles
162 units="metal", # str - Which units used; some potentials
163 # require certain units
164 atom_style="atomic",
165 special_bonds=None,
166 # potential informations
167 pair_style="lj/cut 2.5",
168 pair_coeff=["* * 1 1"],
169 masses=None,
170 pair_modify=None,
171 # variables controlling the output
172 thermo_args=[
173 "step", "temp", "press", "cpu", "pxx", "pyy", "pzz",
174 "pxy", "pxz", "pyz", "ke", "pe", "etotal", "vol", "lx",
175 "ly", "lz", "atoms", ],
176 dump_properties=["id", "type", "x", "y", "z", "vx", "vy",
177 "vz", "fx", "fy", "fz", ],
178 dump_period=1, # period of system snapshot saving (in MD steps)
179 )
181 default_parameters = dict(ase_parameters, **lammps_parameters)
183 def __init__(self, label="lammps", **kwargs):
184 super().__init__(label=label, **kwargs)
186 self.prism = None
187 self.calls = 0
188 self.forces = None
189 # thermo_content contains data "written by" thermo_style.
190 # It is a list of dictionaries, each dict (one for each line
191 # printed by thermo_style) contains a mapping between each
192 # custom_thermo_args-argument and the corresponding
193 # value as printed by lammps. thermo_content will be
194 # re-populated by the read_log method.
195 self.thermo_content = []
197 if self.parameters['tmp_dir'] is not None:
198 # If tmp_dir is pointing somewhere, don't remove stuff!
199 self.parameters['keep_tmp_files'] = True
200 self._lmp_handle = None # To handle the lmp process
202 if self.parameters['tmp_dir'] is None:
203 self.parameters['tmp_dir'] = mkdtemp(prefix="LAMMPS-")
204 else:
205 self.parameters['tmp_dir'] = os.path.realpath(
206 self.parameters['tmp_dir'])
207 if not os.path.isdir(self.parameters['tmp_dir']):
208 os.mkdir(self.parameters['tmp_dir'], 0o755)
210 for f in self.parameters['files']:
211 shutil.copy(
212 f, os.path.join(self.parameters['tmp_dir'],
213 os.path.basename(f)))
215 def get_lammps_command(self):
216 cmd = self.parameters.get('command')
218 if cmd is None:
219 from ase.config import cfg
220 envvar = f'ASE_{self.name.upper()}_COMMAND'
221 cmd = cfg.get(envvar)
223 if cmd is None:
224 # TODO deprecate and remove guesswork
225 cmd = 'lammps'
227 opts = self.parameters.get('lammps_options')
229 if opts is not None:
230 cmd = f'{cmd} {opts}'
232 return cmd
234 def clean(self, force=False):
236 self._lmp_end()
238 if not self.parameters['keep_tmp_files'] or force:
239 shutil.rmtree(self.parameters['tmp_dir'])
241 def check_state(self, atoms, tol=1.0e-10):
242 # Transforming the unit cell to conform to LAMMPS' convention for
243 # orientation (c.f. https://lammps.sandia.gov/doc/Howto_triclinic.html)
244 # results in some precision loss, so we use bit larger tolerance than
245 # machine precision here. Note that there can also be precision loss
246 # related to how many significant digits are specified for things in
247 # the LAMMPS input file.
248 return Calculator.check_state(self, atoms, tol)
250 def calculate(self, atoms=None, properties=None, system_changes=None):
251 if properties is None:
252 properties = self.implemented_properties
253 if system_changes is None:
254 system_changes = all_changes
255 Calculator.calculate(self, atoms, properties, system_changes)
256 self.run()
258 def _lmp_alive(self):
259 # Return True if this calculator is currently handling a running
260 # lammps process
261 return self._lmp_handle and not isinstance(
262 self._lmp_handle.poll(), int
263 )
265 def _lmp_end(self):
266 # Close lammps input and wait for lammps to end. Return process
267 # return value
268 if self._lmp_alive():
269 # !TODO: handle lammps error codes
270 try:
271 self._lmp_handle.communicate(timeout=5)
272 except subprocess.TimeoutExpired:
273 self._lmp_handle.kill()
274 self._lmp_handle.communicate()
275 err = self._lmp_handle.poll()
276 assert err is not None
277 return err
279 def set_missing_parameters(self):
280 """Verify that all necessary variables are set.
281 """
282 symbols = self.atoms.get_chemical_symbols()
283 # If unspecified default to atom types in alphabetic order
284 if not self.parameters.get('specorder'):
285 self.parameters['specorder'] = sorted(set(symbols))
287 # !TODO: handle cases were setting masses actual lead to errors
288 if not self.parameters.get('masses'):
289 self.parameters['masses'] = []
290 for type_id, specie in enumerate(self.parameters['specorder']):
291 mass = atomic_masses[chemical_symbols.index(specie)]
292 self.parameters['masses'] += [
293 f"{type_id + 1:d} {mass:f}"
294 ]
296 # set boundary condtions
297 if not self.parameters.get('boundary'):
298 b_str = " ".join(["fp"[int(x)] for x in self.atoms.pbc])
299 self.parameters['boundary'] = b_str
301 def run(self, set_atoms=False):
302 # !TODO: split this function
303 """Method which explicitly runs LAMMPS."""
304 pbc = self.atoms.get_pbc()
305 if all(pbc):
306 cell = self.atoms.get_cell()
307 elif not any(pbc):
308 # large enough cell for non-periodic calculation -
309 # LAMMPS shrink-wraps automatically via input command
310 # "periodic s s s"
311 # below
312 cell = 2 * np.max(np.abs(self.atoms.get_positions())) * np.eye(3)
313 else:
314 warnings.warn(
315 "semi-periodic ASE cell detected - translation "
316 + "to proper LAMMPS input cell might fail"
317 )
318 cell = self.atoms.get_cell()
319 self.prism = Prism(cell)
321 self.set_missing_parameters()
322 self.calls += 1
324 # change into subdirectory for LAMMPS calculations
325 tempdir = self.parameters['tmp_dir']
327 # setup file names for LAMMPS calculation
328 label = f"{self.label}{self.calls:>06}"
329 lammps_in = uns_mktemp(
330 prefix="in_" + label, dir=tempdir
331 )
332 lammps_log = uns_mktemp(
333 prefix="log_" + label, dir=tempdir
334 )
335 lammps_trj_fd = NamedTemporaryFile(
336 prefix="trj_" + label,
337 suffix=(".bin" if self.parameters['binary_dump'] else ""),
338 dir=tempdir,
339 delete=(not self.parameters['keep_tmp_files']),
340 )
341 lammps_trj = lammps_trj_fd.name
342 if self.parameters['no_data_file']:
343 lammps_data = None
344 else:
345 lammps_data_fd = NamedTemporaryFile(
346 prefix="data_" + label,
347 dir=tempdir,
348 delete=(not self.parameters['keep_tmp_files']),
349 mode='w',
350 encoding='ascii'
351 )
352 write_lammps_data(
353 lammps_data_fd,
354 self.atoms,
355 specorder=self.parameters['specorder'],
356 force_skew=self.parameters['always_triclinic'],
357 reduce_cell=self.parameters['reduce_cell'],
358 velocities=self.parameters['write_velocities'],
359 prismobj=self.prism,
360 units=self.parameters['units'],
361 atom_style=self.parameters['atom_style'],
362 )
363 lammps_data = lammps_data_fd.name
364 lammps_data_fd.flush()
366 # see to it that LAMMPS is started
367 if not self._lmp_alive():
368 command = self.get_lammps_command()
369 # Attempt to (re)start lammps
370 self._lmp_handle = subprocess.Popen(
371 shlex.split(command, posix=(os.name == "posix")),
372 stdin=subprocess.PIPE,
373 stdout=subprocess.PIPE,
374 encoding='ascii',
375 )
376 lmp_handle = self._lmp_handle
378 # Create thread reading lammps stdout (for reference, if requested,
379 # also create lammps_log, although it is never used)
380 if self.parameters['keep_tmp_files']:
381 lammps_log_fd = open(lammps_log, "w")
382 fd = SpecialTee(lmp_handle.stdout, lammps_log_fd)
383 else:
384 fd = lmp_handle.stdout
385 thr_read_log = Thread(target=self.read_lammps_log, args=(fd,))
386 thr_read_log.start()
388 # write LAMMPS input (for reference, also create the file lammps_in,
389 # although it is never used)
390 if self.parameters['keep_tmp_files']:
391 lammps_in_fd = open(lammps_in, "w")
392 fd = SpecialTee(lmp_handle.stdin, lammps_in_fd)
393 else:
394 fd = lmp_handle.stdin
395 write_lammps_in(
396 lammps_in=fd,
397 parameters=self.parameters,
398 atoms=self.atoms,
399 prismobj=self.prism,
400 lammps_trj=lammps_trj,
401 lammps_data=lammps_data,
402 )
404 if self.parameters['keep_tmp_files']:
405 lammps_in_fd.close()
407 # Wait for log output to be read (i.e., for LAMMPS to finish)
408 # and close the log file if there is one
409 thr_read_log.join()
410 if self.parameters['keep_tmp_files']:
411 lammps_log_fd.close()
413 if not self.parameters['keep_alive']:
414 self._lmp_end()
416 exitcode = lmp_handle.poll()
417 if exitcode and exitcode != 0:
418 raise RuntimeError(
419 "LAMMPS exited in {} with exit code: {}."
420 "".format(tempdir, exitcode)
421 )
423 # A few sanity checks
424 if len(self.thermo_content) == 0:
425 raise RuntimeError("Failed to retrieve any thermo_style-output")
426 if int(self.thermo_content[-1]["atoms"]) != len(self.atoms):
427 # This obviously shouldn't happen, but if prism.fold_...() fails,
428 # it could
429 raise RuntimeError("Atoms have gone missing")
431 trj_atoms = read_lammps_dump(
432 infileobj=lammps_trj,
433 order=self.parameters['atorder'],
434 index=-1,
435 prismobj=self.prism,
436 specorder=self.parameters['specorder'],
437 )
439 if set_atoms:
440 self.atoms = trj_atoms.copy()
442 self.forces = trj_atoms.get_forces()
443 # !TODO: trj_atoms is only the last snapshot of the system; Is it
444 # desirable to save also the inbetween steps?
445 if self.parameters['trajectory_out'] is not None:
446 # !TODO: is it advisable to create here temporary atoms-objects
447 self.trajectory_out.write(trj_atoms)
449 tc = self.thermo_content[-1]
450 self.results["energy"] = convert(
451 tc["pe"], "energy", self.parameters["units"], "ASE"
452 )
453 self.results["free_energy"] = self.results["energy"]
454 self.results['forces'] = convert(self.forces.copy(),
455 'force',
456 self.parameters['units'],
457 'ASE')
458 stress = np.array(
459 [-tc[i] for i in ("pxx", "pyy", "pzz", "pyz", "pxz", "pxy")]
460 )
462 # We need to apply the Lammps rotation stuff to the stress:
463 xx, yy, zz, yz, xz, xy = stress
464 stress_tensor = np.array([[xx, xy, xz],
465 [xy, yy, yz],
466 [xz, yz, zz]])
467 stress_atoms = self.prism.tensor2_to_ase(stress_tensor)
468 stress_atoms = stress_atoms[[0, 1, 2, 1, 0, 0],
469 [0, 1, 2, 2, 2, 1]]
470 stress = stress_atoms
472 self.results["stress"] = convert(
473 stress, "pressure", self.parameters["units"], "ASE"
474 )
476 lammps_trj_fd.close()
477 if not self.parameters['no_data_file']:
478 lammps_data_fd.close()
480 def __enter__(self):
481 return self
483 def __exit__(self, *args):
484 self._lmp_end()
486 def read_lammps_log(self, fileobj):
487 # !TODO: somehow communicate 'thermo_content' explicitly
488 """Method which reads a LAMMPS output log file."""
490 # read_log depends on that the first (three) thermo_style custom args
491 # can be capitalized and matched against the log output. I.e.
492 # don't use e.g. 'ke' or 'cpu' which are labeled KinEng and CPU.
493 mark_re = r"^\s*" + r"\s+".join(
494 [x.capitalize() for x in self.parameters['thermo_args'][0:3]]
495 )
496 _custom_thermo_mark = re_compile(mark_re)
498 # !TODO: regex-magic necessary?
499 # Match something which can be converted to a float
500 f_re = r"([+-]?(?:(?:\d+(?:\.\d*)?|\.\d+)(?:e[+-]?\d+)?|nan|inf))"
501 n_args = len(self.parameters["thermo_args"])
502 # Create a re matching exactly N white space separated floatish things
503 _custom_thermo_re = re_compile(
504 r"^\s*" + r"\s+".join([f_re] * n_args) + r"\s*$", flags=IGNORECASE
505 )
507 thermo_content = []
508 line = fileobj.readline()
509 while line and line.strip() != CALCULATION_END_MARK:
510 # check error
511 if 'ERROR:' in line:
512 raise RuntimeError(f'LAMMPS exits with error message: {line}')
514 # get thermo output
515 if _custom_thermo_mark.match(line):
516 while True:
517 line = fileobj.readline()
518 if 'WARNING:' in line:
519 continue
521 bool_match = _custom_thermo_re.match(line)
522 if not bool_match:
523 break
525 # create a dictionary between each of the
526 # thermo_style args and it's corresponding value
527 thermo_content.append(
528 dict(
529 zip(
530 self.parameters['thermo_args'],
531 map(float, bool_match.groups()),
532 )
533 )
534 )
535 else:
536 line = fileobj.readline()
538 self.thermo_content = thermo_content
541class SpecialTee:
542 """A special purpose, with limited applicability, tee-like thing.
544 A subset of stuff read from, or written to, orig_fd,
545 is also written to out_fd.
546 It is used by the lammps calculator for creating file-logs of stuff
547 read from, or written to, stdin and stdout, respectively.
548 """
550 def __init__(self, orig_fd, out_fd):
551 self._orig_fd = orig_fd
552 self._out_fd = out_fd
553 self.name = orig_fd.name
555 def write(self, data):
556 self._orig_fd.write(data)
557 self._out_fd.write(data)
558 self.flush()
560 def read(self, *args, **kwargs):
561 data = self._orig_fd.read(*args, **kwargs)
562 self._out_fd.write(data)
563 return data
565 def readline(self, *args, **kwargs):
566 data = self._orig_fd.readline(*args, **kwargs)
567 self._out_fd.write(data)
568 return data
570 def readlines(self, *args, **kwargs):
571 data = self._orig_fd.readlines(*args, **kwargs)
572 self._out_fd.write("".join(data))
573 return data
575 def flush(self):
576 self._orig_fd.flush()
577 self._out_fd.flush()