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