Coverage for ase/calculators/lammpsrun.py: 87.50%

1# fmt: off

3"""ASE calculator for the LAMMPS classical MD code"""

4# lammps.py (2011/03/29)

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 <https://www.gnu.org/licenses/>.

24import os

25import shlex

26import shutil

27import subprocess

28import warnings

29from contextlib import ExitStack

30from re import IGNORECASE

31from re import compile as re_compile

32from tempfile import NamedTemporaryFile, mkdtemp

33from tempfile import mktemp as uns_mktemp

34from threading import Thread

35from typing import Any, Dict

37import numpy as np

39from ase.calculators.calculator import Calculator, all_changes

40from ase.calculators.lammps import (

41 CALCULATION_END_MARK,

42 Prism,

43 convert,

44 write_lammps_in,

45)

46from ase.data import atomic_masses, chemical_symbols

47from ase.io.lammpsdata import write_lammps_data

48from ase.io.lammpsrun import read_lammps_dump

50__all__ = ["LAMMPS"]

53class LAMMPS(Calculator):

54 """LAMMPS (https://lammps.sandia.gov/) calculator

56 The environment variable :envvar:`ASE_LAMMPSRUN_COMMAND` must be defined to

57 tell ASE how to call a LAMMPS binary. This should contain the path to the

58 lammps binary, or more generally, a command line possibly also including an

59 MPI-launcher command.

61 For example (in a Bourne-shell compatible environment):

63 .. code-block:: bash

65 export ASE_LAMMPSRUN_COMMAND=/path/to/lmp_binary

67 or possibly something similar to

69 .. code-block:: bash

71 export ASE_LAMMPSRUN_COMMAND="/path/to/mpirun --np 4 lmp_binary"

73 Parameters

74 ----------

75 files : list[str]

76 List of files needed by LAMMPS. Typically a list of potential files.

77 parameters : dict[str, Any]

78 Dictionary of settings to be passed into the input file for calculation.

79 specorder : list[str]

80 Within LAAMPS, atoms are identified by an integer value starting from 1.

81 This variable allows the user to define the order of the indices

82 assigned to the atoms in the calculation, with the default

83 if not given being alphabetical

84 keep_tmp_files : bool, default: False

85 Retain any temporary files created. Mostly useful for debugging.

86 tmp_dir : str, default: None

87 path/dirname (default None -> create automatically).

88 Explicitly control where the calculator object should create

89 its files. Using this option implies 'keep_tmp_files=True'.

90 no_data_file : bool, default: False

91 Controls whether an explicit data file will be used for feeding

92 atom coordinates into lammps. Enable it to lessen the pressure on

93 the (tmp) file system. THIS OPTION MIGHT BE UNRELIABLE FOR CERTAIN

94 CORNER CASES (however, if it fails, you will notice...).

95 keep_alive : bool, default: True

96 When using LAMMPS as a spawned subprocess, keep the subprocess

97 alive (but idling when unused) along with the calculator object.

98 always_triclinic : bool, default: False

99 Force LAMMPS to treat the cell as tilted, even if the cell is not

100 tilted, by printing ``xy``, ``xz``, ``yz`` in the data file.

101 reduce_cell : bool, default: False

102 If True, reduce cell to have shorter lattice vectors.

103 write_velocities : bool, default: False

104 If True, forward ASE velocities to LAMMPS.

105 verbose: bool, default: False

106 If True, print additional debugging output to STDOUT.

107

108 Examples

109 --------

110 Provided that the respective potential file is in the working directory,

111 one can simply run (note that LAMMPS needs to be compiled to work with EAM

112 potentials)

113

114 ::

115

116 from ase import Atom, Atoms

117 from ase.build import bulk

118 from ase.calculators.lammpsrun import LAMMPS

119

120 parameters = {'pair_style': 'eam/alloy',

121 'pair_coeff': ['* * NiAlH_jea.eam.alloy H Ni']}

122

123 files = ['NiAlH_jea.eam.alloy']

124

125 Ni = bulk('Ni', cubic=True)

126 H = Atom('H', position=Ni.cell.diagonal()/2)

127 NiH = Ni + H

128

129 lammps = LAMMPS(files=files, **parameters)

130

131 NiH.calc = lammps

132 print("Energy ", NiH.get_potential_energy())

133 """

134

135 name = "lammpsrun"

136 implemented_properties = ["energy", "free_energy", "forces", "stress",

137 "energies"]

138

139 # parameters to choose options in LAMMPSRUN

140 ase_parameters: Dict[str, Any] = dict(

141 specorder=None,

142 atorder=True,

143 always_triclinic=False,

144 reduce_cell=False,

145 keep_alive=True,

146 keep_tmp_files=False,

147 no_data_file=False,

148 tmp_dir=None,

149 files=[], # usually contains potential parameters

150 verbose=False,

151 write_velocities=False,

152 binary_dump=True, # bool - use binary dump files (full

153 # precision but long long ids are casted to

154 # double)

155 lammps_options="-echo log -screen none -log /dev/stdout",

156 trajectory_out=None, # file object, if is not None the

157 # trajectory will be saved in it

158 )

159

160 # parameters forwarded to LAMMPS

161 lammps_parameters = dict(

162 boundary=None, # bounadry conditions styles

163 units="metal", # str - Which units used; some potentials

164 # require certain units

165 atom_style="atomic",

166 special_bonds=None,

167 # potential informations

168 pair_style="lj/cut 2.5",

169 pair_coeff=["* * 1 1"],

170 masses=None,

171 pair_modify=None,

172 # variables controlling the output

173 thermo_args=[

174 "step", "temp", "press", "cpu", "pxx", "pyy", "pzz",

175 "pxy", "pxz", "pyz", "ke", "pe", "etotal", "vol", "lx",

176 "ly", "lz", "atoms", ],

177 dump_properties=["id", "type", "x", "y", "z", "vx", "vy",

178 "vz", "fx", "fy", "fz", ],

179 dump_period=1, # period of system snapshot saving (in MD steps)

180 )

181

182 default_parameters = dict(ase_parameters, **lammps_parameters)

183

184 def __init__(self, label="lammps", **kwargs):

185 super().__init__(label=label, **kwargs)

186

187 self.prism = None

188 self.calls = 0

189 self.forces = None

190 # thermo_content contains data "written by" thermo_style.

191 # It is a list of dictionaries, each dict (one for each line

192 # printed by thermo_style) contains a mapping between each

193 # custom_thermo_args-argument and the corresponding

194 # value as printed by lammps. thermo_content will be

195 # re-populated by the read_log method.

196 self.thermo_content = []

197

198 if self.parameters['tmp_dir'] is not None:

199 # If tmp_dir is pointing somewhere, don't remove stuff!

200 self.parameters['keep_tmp_files'] = True

201 self._lmp_handle = None # To handle the lmp process

202

203 if self.parameters['tmp_dir'] is None:

204 self.parameters['tmp_dir'] = mkdtemp(prefix="LAMMPS-")

205 else:

206 self.parameters['tmp_dir'] = os.path.realpath(

207 self.parameters['tmp_dir'])

208 if not os.path.isdir(self.parameters['tmp_dir']):

209 os.mkdir(self.parameters['tmp_dir'], 0o755)

210

211 for f in self.parameters['files']:

212 shutil.copy(

213 f, os.path.join(self.parameters['tmp_dir'],

214 os.path.basename(f)))

215

216 def get_lammps_command(self):

217 cmd = self.parameters.get('command')

218

219 if cmd is None:

220 from ase.config import cfg

221 envvar = f'ASE_{self.name.upper()}_COMMAND'

222 cmd = cfg.get(envvar)

223

224 if cmd is None:

225 # TODO deprecate and remove guesswork

226 cmd = 'lammps'

227

228 opts = self.parameters.get('lammps_options')

229

230 if opts is not None:

231 cmd = f'{cmd} {opts}'

232

233 return cmd

234

235 def clean(self, force=False):

236

237 self._lmp_end()

238

239 if not self.parameters['keep_tmp_files'] or force:

240 shutil.rmtree(self.parameters['tmp_dir'])

241

242 def check_state(self, atoms, tol=1.0e-10):

243 # Transforming the unit cell to conform to LAMMPS' convention for

244 # orientation (c.f. https://lammps.sandia.gov/doc/Howto_triclinic.html)

245 # results in some precision loss, so we use bit larger tolerance than

246 # machine precision here. Note that there can also be precision loss

247 # related to how many significant digits are specified for things in

248 # the LAMMPS input file.

249 return Calculator.check_state(self, atoms, tol)

250

251 def calculate(self, atoms=None, properties=None, system_changes=None):

252 if properties is None:

253 properties = self.implemented_properties

254 if system_changes is None:

255 system_changes = all_changes

256 Calculator.calculate(self, atoms, properties, system_changes)

257 self.run()

258

259 def _lmp_alive(self):

260 # Return True if this calculator is currently handling a running

261 # lammps process

262 return self._lmp_handle and not isinstance(

263 self._lmp_handle.poll(), int

264 )

265

266 def _lmp_end(self):

267 # Close lammps input and wait for lammps to end. Return process

268 # return value

269 if self._lmp_alive():

270 # !TODO: handle lammps error codes

271 try:

272 self._lmp_handle.communicate(timeout=5)

273 except subprocess.TimeoutExpired:

274 self._lmp_handle.kill()

275 self._lmp_handle.communicate()

276 err = self._lmp_handle.poll()

277 assert err is not None

278 return err

279

280 def set_missing_parameters(self):

281 """Verify that all necessary variables are set.

282 """

283 symbols = self.atoms.get_chemical_symbols()

284 # If unspecified default to atom types in alphabetic order

285 if not self.parameters.get('specorder'):

286 self.parameters['specorder'] = sorted(set(symbols))

287

288 # !TODO: handle cases were setting masses actual lead to errors

289 if not self.parameters.get('masses'):

290 self.parameters['masses'] = []

291 for type_id, specie in enumerate(self.parameters['specorder']):

292 mass = atomic_masses[chemical_symbols.index(specie)]

293 self.parameters['masses'] += [

294 f"{type_id + 1:d} {mass:f}"

295 ]

296

297 # set boundary condtions

298 if not self.parameters.get('boundary'):

299 b_str = " ".join(["fp"[int(x)] for x in self.atoms.pbc])

300 self.parameters['boundary'] = b_str

301

302 def run(self, set_atoms=False):

303 # !TODO: split this function

304 """Method which explicitly runs LAMMPS."""

305

306 with ExitStack() as exitstack:

307 self._run(set_atoms, exitstack)

308

309 def _run(self, set_atoms, exitstack):

310 pbc = self.atoms.get_pbc()

311 if all(pbc):

312 cell = self.atoms.get_cell()

313 elif not any(pbc):

314 # large enough cell for non-periodic calculation -

315 # LAMMPS shrink-wraps automatically via input command

316 # "periodic s s s"

317 # below

318 cell = 2 * np.max(np.abs(self.atoms.get_positions())) * np.eye(3)

319 else:

320 warnings.warn(

321 "semi-periodic ASE cell detected - translation "

322 + "to proper LAMMPS input cell might fail"

323 )

324 cell = self.atoms.get_cell()

325 self.prism = Prism(cell)

326

327 self.set_missing_parameters()

328 self.calls += 1

329

330 # change into subdirectory for LAMMPS calculations

331 tempdir = self.parameters['tmp_dir']

332

333 # setup file names for LAMMPS calculation

334 label = f"{self.label}{self.calls:>06}"

335 lammps_in = uns_mktemp(

336 prefix="in_" + label, dir=tempdir

337 )

338 lammps_log = uns_mktemp(

339 prefix="log_" + label, dir=tempdir

340 )

341 lammps_trj_fd = exitstack.enter_context(NamedTemporaryFile(

342 prefix="trj_" + label,

343 suffix=(".bin" if self.parameters['binary_dump'] else ""),

344 dir=tempdir,

345 delete=(not self.parameters['keep_tmp_files']),

346 ))

347 lammps_trj = lammps_trj_fd.name

348 if self.parameters['no_data_file']:

349 lammps_data = None

350 else:

351 lammps_data_fd = exitstack.enter_context(NamedTemporaryFile(

352 prefix="data_" + label,

353 dir=tempdir,

354 delete=(not self.parameters['keep_tmp_files']),

355 mode='w',

356 encoding='ascii'

357 ))

358 write_lammps_data(

359 lammps_data_fd,

360 self.atoms,

361 specorder=self.parameters['specorder'],

362 force_skew=self.parameters['always_triclinic'],

363 reduce_cell=self.parameters['reduce_cell'],

364 velocities=self.parameters['write_velocities'],

365 prismobj=self.prism,

366 units=self.parameters['units'],

367 atom_style=self.parameters['atom_style'],

368 )

369 lammps_data = lammps_data_fd.name

370 lammps_data_fd.flush()

371

372 # see to it that LAMMPS is started

373 if not self._lmp_alive():

374 command = self.get_lammps_command()

375 # Attempt to (re)start lammps

376 self._lmp_handle = subprocess.Popen(

377 shlex.split(command, posix=(os.name == "posix")),

378 stdin=subprocess.PIPE,

379 stdout=subprocess.PIPE,

380 text=True,

381 )

382 lmp_handle = self._lmp_handle

383

384 # Create thread reading lammps stdout (for reference, if requested,

385 # also create lammps_log, although it is never used)

386 if self.parameters['keep_tmp_files']:

387 lammps_log_fd = exitstack.enter_context(open(lammps_log, "w"))

388 fd = SpecialTee(lmp_handle.stdout, lammps_log_fd)

389 else:

390 fd = lmp_handle.stdout

391 thr_read_log = Thread(target=self.read_lammps_log, args=(fd,))

392 thr_read_log.start()

393

394 # write LAMMPS input (for reference, also create the file lammps_in,

395 # although it is never used)

396 if self.parameters['keep_tmp_files']:

397 lammps_in_fd = exitstack.enter_context(open(lammps_in, "w"))

398 fd = SpecialTee(lmp_handle.stdin, lammps_in_fd)

399 else:

400 fd = lmp_handle.stdin

401 write_lammps_in(

402 lammps_in=fd,

403 parameters=self.parameters,

404 atoms=self.atoms,

405 prismobj=self.prism,

406 lammps_trj=lammps_trj,

407 lammps_data=lammps_data,

408 )

409

410 # Wait for log output to be read (i.e., for LAMMPS to finish)

411 # and close the log file if there is one

412 thr_read_log.join()

413

414 if not self.parameters['keep_alive']:

415 self._lmp_end()

416

417 exitcode = lmp_handle.poll()

418 if exitcode and exitcode != 0:

419 raise RuntimeError(

420 "LAMMPS exited in {} with exit code: {}."

421 "".format(tempdir, exitcode)

422 )

423

424 # A few sanity checks

425 if len(self.thermo_content) == 0:

426 raise RuntimeError("Failed to retrieve any thermo_style-output")

427 if int(self.thermo_content[-1]["atoms"]) != len(self.atoms):

428 # This obviously shouldn't happen, but if prism.fold_...() fails,

429 # it could

430 raise RuntimeError("Atoms have gone missing")

431

432 trj_atoms = read_lammps_dump(

433 infileobj=lammps_trj,

434 order=self.parameters['atorder'],

435 index=-1,

436 prismobj=self.prism,

437 specorder=self.parameters['specorder'],

438 )

439

440 if set_atoms:

441 self.atoms = trj_atoms.copy()

442

443 self.forces = trj_atoms.get_forces()

444 # !TODO: trj_atoms is only the last snapshot of the system; Is it

445 # desirable to save also the inbetween steps?

446 if self.parameters['trajectory_out'] is not None:

447 # !TODO: is it advisable to create here temporary atoms-objects

448 self.trajectory_out.write(trj_atoms)

449

450 tc = self.thermo_content[-1]

451 self.results["energy"] = convert(

452 tc["pe"], "energy", self.parameters["units"], "ASE"

453 )

454 self.results["free_energy"] = self.results["energy"]

455 self.results['forces'] = convert(self.forces.copy(),

456 'force',

457 self.parameters['units'],

458 'ASE')

459 stress = np.array(

460 [-tc[i] for i in ("pxx", "pyy", "pzz", "pyz", "pxz", "pxy")]

461 )

462

463 # We need to apply the Lammps rotation stuff to the stress:

464 xx, yy, zz, yz, xz, xy = stress

465 stress_tensor = np.array([[xx, xy, xz],

466 [xy, yy, yz],

467 [xz, yz, zz]])

468 stress_atoms = self.prism.tensor2_to_ase(stress_tensor)

469 stress_atoms = stress_atoms[[0, 1, 2, 1, 0, 0],

470 [0, 1, 2, 2, 2, 1]]

471 stress = stress_atoms

472

473 self.results["stress"] = convert(

474 stress, "pressure", self.parameters["units"], "ASE"

475 )

476

477 def __enter__(self):

478 return self

479

480 def __exit__(self, *args):

481 self._lmp_end()

482

483 def read_lammps_log(self, fileobj):

484 # !TODO: somehow communicate 'thermo_content' explicitly

485 """Method which reads a LAMMPS output log file."""

486

487 # read_log depends on that the first (three) thermo_style custom args

488 # can be capitalized and matched against the log output. I.e.

489 # don't use e.g. 'ke' or 'cpu' which are labeled KinEng and CPU.

490 mark_re = r"^\s*" + r"\s+".join(

491 [x.capitalize() for x in self.parameters['thermo_args'][0:3]]

492 )

493 _custom_thermo_mark = re_compile(mark_re)

494

495 # !TODO: regex-magic necessary?

496 # Match something which can be converted to a float

497 f_re = r"([+-]?(?:(?:\d+(?:\.\d*)?|\.\d+)(?:e[+-]?\d+)?|nan|inf))"

498 n_args = len(self.parameters["thermo_args"])

499 # Create a re matching exactly N white space separated floatish things

500 _custom_thermo_re = re_compile(

501 r"^\s*" + r"\s+".join([f_re] * n_args) + r"\s*$", flags=IGNORECASE

502 )

503

504 thermo_content = []

505 line = fileobj.readline()

506 while line and line.strip() != CALCULATION_END_MARK:

507 # check error

508 if 'ERROR:' in line:

509 raise RuntimeError(f'LAMMPS exits with error message: {line}')

510

511 # get thermo output

512 if _custom_thermo_mark.match(line):

513 while True:

514 line = fileobj.readline()

515 if 'WARNING:' in line:

516 continue

517

518 bool_match = _custom_thermo_re.match(line)

519 if not bool_match:

520 break

521

522 # create a dictionary between each of the

523 # thermo_style args and it's corresponding value

524 thermo_content.append(

525 dict(

526 zip(

527 self.parameters['thermo_args'],

528 map(float, bool_match.groups()),

529 )

530 )

531 )

532 else:

533 line = fileobj.readline()

534

535 self.thermo_content = thermo_content

536

537

538class SpecialTee:

539 """A special purpose, with limited applicability, tee-like thing.

540

541 A subset of stuff read from, or written to, orig_fd,

542 is also written to out_fd.

543 It is used by the lammps calculator for creating file-logs of stuff

544 read from, or written to, stdin and stdout, respectively.

545 """

546

547 def __init__(self, orig_fd, out_fd):

548 self._orig_fd = orig_fd

549 self._out_fd = out_fd

550 self.name = orig_fd.name

551

552 def write(self, data):

553 self._orig_fd.write(data)

554 self._out_fd.write(data)

555 self.flush()

556

557 def read(self, *args, **kwargs):

558 data = self._orig_fd.read(*args, **kwargs)

559 self._out_fd.write(data)

560 return data

561

562 def readline(self, *args, **kwargs):

563 data = self._orig_fd.readline(*args, **kwargs)

564 self._out_fd.write(data)

565 return data

566

567 def readlines(self, *args, **kwargs):

568 data = self._orig_fd.readlines(*args, **kwargs)

569 self._out_fd.write("".join(data))

570 return data

571

572 def flush(self):

573 self._orig_fd.flush()

574 self._out_fd.flush()

Coverage for ase / calculators / lammpsrun.py: 87.50%

208 statements