Coverage for ase / calculators / lammpslib.py: 79.85%
263 statements
« prev ^ index » next coverage.py v7.13.3, created at 2026-02-04 10:20 +0000
1# fmt: off
3"""ASE LAMMPS Calculator Library Version"""
6import ctypes
8import numpy as np
10from ase import Atoms
11from ase.calculators.calculator import Calculator
12from ase.calculators.lammps import Prism, convert
13from ase.data import atomic_masses as ase_atomic_masses
14from ase.data import atomic_numbers as ase_atomic_numbers
15from ase.data import chemical_symbols as ase_chemical_symbols
17# TODO
18# 1. should we make a new lammps object each time ?
19# 4. need a routine to get the model back from lammps
20# 5. if we send a command to lmps directly then the calculator does
21# not know about it and the energy could be wrong.
22# 6. do we need a subroutine generator that converts a lammps string
23# into a python function that can be called
24# 8. make matscipy as fallback
25# 9. keep_alive not needed with no system changes
28class LAMMPSlib(Calculator):
29 r"""
30**Introduction**
32LAMMPSlib is an interface and calculator for LAMMPS_. LAMMPSlib uses
33the python interface that comes with LAMMPS to solve an atoms model
34for energy, atom forces and cell stress. This calculator creates a
35'.lmp' object which is a running lammps program, so further commands
36can be sent to this object executed until it is explicitly closed. Any
37additional variables calculated by lammps can also be extracted. This
38is still experimental code.
40**Arguments**
42======================= ======================================================
43Keyword Description
44======================= ======================================================
45``lmpcmds`` list of strings of LAMMPS commands. You need to supply
46 enough to define the potential to be used e.g.
48 ["pair_style eam/alloy",
49 "pair_coeff * * potentials/NiAlH_jea.eam.alloy Ni Al"]
51``atom_types`` dictionary of ``atomic_symbol :lammps_atom_type``
52 pairs, e.g. ``{'Cu':1}`` to bind copper to lammps
53 atom type 1. If <None>, autocreated by assigning
54 lammps atom types in order that they appear in the
55 first used atoms object.
57``atom_type_masses`` dictionary of ``atomic_symbol :mass`` pairs, e.g.
58 ``{'Cu':63.546}`` to optionally assign masses that
59 override default ase.data.atomic_masses. Note that
60 since unit conversion is done automatically in this
61 module, these quantities must be given in the
62 standard ase mass units (g/mol)
64``log_file`` string
65 path to the desired LAMMPS log file
67``lammps_header`` string to use for lammps setup. Default is to use
68 metal units and simple atom simulation.
70 lammps_header=['units metal',
71 'atom_style atomic',
72 'atom_modify map array sort 0 0'])
74``amendments`` extra list of strings of LAMMPS commands to be run
75 post initialization. (Use: Initialization amendments)
76 e.g.
78 ["mass 1 58.6934"]
80``post_changebox_cmds`` extra list of strings of LAMMPS commands to be run
81 after any LAMMPS 'change_box' command is performed by
82 the calculator. This is relevant because some
83 potentials either themselves depend on the geometry
84 and boundary conditions of the simulation box, or are
85 frequently coupled with other LAMMPS commands that
86 do, e.g. the 'buck/coul/long' pair style is often
87 used with the kspace_* commands, which are sensitive
88 to the periodicity of the simulation box.
90``extra_cmd_args`` list of extra arguments for
91 `lammps.lammps(cmd_args=...)`, e.g. for kokkos mliap on
92 a gpu
93 ```
94 ("-k on g 1 -sf kk "
95 "-pk kokkos neigh half newton on").split()
96 ```
98``intializer`` callback function that does arbitrary LAMMPS python
99 API initializtion tasks (e.g. calling
100 `lammps.mliap.activate_mliapy`) and accepts a single
101 positional argument, `self.lmp`.
103``keep_alive`` Boolean
104 whether to keep the lammps routine alive for more
105 commands. Default is True.
107======================= ======================================================
110**Requirements**
112To run this calculator you must have LAMMPS installed and compiled to
113enable the python interface. See the LAMMPS manual.
115If the following code runs then lammps is installed correctly.
117 >>> from lammps import lammps
118 >>> lmp = lammps()
120The version of LAMMPS is also important. LAMMPSlib is suitable for
121versions after approximately 2011. Prior to this the python interface
122is slightly different from that used by LAMMPSlib. It is not difficult
123to change to the earlier format.
125**LAMMPS and LAMMPSlib**
127The LAMMPS calculator is another calculator that uses LAMMPS (the
128program) to calculate the energy by generating input files and running
129a separate LAMMPS job to perform the analysis. The output data is then
130read back into python. LAMMPSlib makes direct use of the LAMMPS (the
131program) python interface. As well as directly running any LAMMPS
132command line it allows the values of any of LAMMPS variables to be
133extracted and returned to python.
135**Example**
137Provided that the respective potential file is in the working directory, one
138can simply run (note that LAMMPS needs to be compiled to work with EAM
139potentials)
141::
143 from ase import Atom, Atoms
144 from ase.build import bulk
145 from ase.calculators.lammpslib import LAMMPSlib
147 cmds = ["pair_style eam/alloy",
148 "pair_coeff * * NiAlH_jea.eam.alloy Ni H"]
150 Ni = bulk('Ni', cubic=True)
151 H = Atom('H', position=Ni.cell.diagonal()/2)
152 NiH = Ni + H
154 lammps = LAMMPSlib(lmpcmds=cmds, log_file='test.log')
156 NiH.calc = lammps
157 print("Energy ", NiH.get_potential_energy())
160**Implementation**
162LAMMPS provides a set of python functions to allow execution of the
163underlying C++ LAMMPS code. The functions used by the LAMMPSlib
164interface are::
166 from lammps import lammps
168 lmp = lammps(cmd_args) # initiate LAMMPS object with command line args
170 lmp.scatter_atoms('x',1,3,positions) # atom coords to LAMMPS C array
171 lmp.command(cmd) # executes a one line cmd string
172 lmp.extract_variable(...) # extracts a per atom variable
173 lmp.extract_global(...) # extracts a global variable
174 lmp.close() # close the lammps object
176For a single Ni atom model the following lammps file commands would be run
177by invoking the get_potential_energy() method::
179 units metal
180 atom_style atomic
181 atom_modify map array sort 0 0
183 region cell prism 0 xhi 0 yhi 0 zhi xy xz yz units box
184 create_box 1 cell
185 create_atoms 1 single 0 0 0 units box
186 mass * 1.0
188 ## user lmpcmds get executed here
189 pair_style eam/alloy
190 pair_coeff * * NiAlH_jea.eam.alloy Ni
191 ## end of user lmmpcmds
193 run 0
195where xhi, yhi and zhi are the lattice vector lengths and xy,
196xz and yz are the tilt of the lattice vectors, all to be edited.
199**Notes**
201.. _LAMMPS: https://lammps.org/
203* Units: The default lammps_header sets the units to Angstrom and eV
204 and for compatibility with ASE Stress is in GPa.
206* The global energy is currently extracted from LAMMPS using
207 extract_variable since lammps.lammps currently extract_global only
208 accepts the following ['dt', 'boxxlo', 'boxxhi', 'boxylo', 'boxyhi',
209 'boxzlo', 'boxzhi', 'natoms', 'nlocal'].
211* If an error occurs while lammps is in control it will crash
212 Python. Check the output of the log file to find the lammps error.
214* If the are commands directly sent to the LAMMPS object this may
215 change the energy value of the model. However the calculator will not
216 know of it and still return the original energy value.
218"""
220 implemented_properties = ['energy', 'free_energy', 'forces', 'stress',
221 'energies']
223 started = False
224 initialized = False
226 default_parameters = dict(
227 atom_types=None,
228 atom_type_masses=None,
229 log_file=None,
230 lammps_name='',
231 keep_alive=True,
232 lammps_header=['units metal',
233 'atom_style atomic',
234 'atom_modify map array sort 0 0'],
235 amendments=None,
236 post_changebox_cmds=None,
237 extra_cmd_args=(),
238 initializer=None,
239 boundary=True,
240 create_box=True,
241 create_atoms=True,
242 read_molecular_info=False,
243 comm=None)
245 def __init__(self, *args, **kwargs):
246 Calculator.__init__(self, *args, **kwargs)
247 self.lmp = None
249 def __enter__(self):
250 return self
252 def __exit__(self, *args):
253 self.clean()
255 def clean(self):
256 if self.started:
257 self.lmp.close()
258 self.started = False
259 self.initialized = False
260 self.lmp = None
262 def set_cell(self, atoms: Atoms, change: bool = False):
263 self.prism = Prism(atoms.cell, atoms.pbc)
264 _ = self.prism.get_lammps_prism()
265 xhi, yhi, zhi, xy, xz, yz = convert(_, "distance", "ASE", self.units)
266 box_hi = [xhi, yhi, zhi]
268 if change:
269 cell_cmd = ('change_box all '
270 'x final 0 {} y final 0 {} z final 0 {} '
271 'xy final {} xz final {} yz final {} units box'
272 ''.format(xhi, yhi, zhi, xy, xz, yz))
273 if self.parameters.post_changebox_cmds is not None:
274 for cmd in self.parameters.post_changebox_cmds:
275 self.lmp.command(cmd)
276 else:
277 # just in case we'll want to run with a funny shape box,
278 # and here command will only happen once, and before
279 # any calculation
280 if self.parameters.create_box:
281 self.lmp.command('box tilt large')
283 # Check if there are any indefinite boundaries. If so,
284 # shrink-wrapping will end up being used, but we want to
285 # define the LAMMPS region and box fairly tight around the
286 # atoms to avoid losing any
287 lammps_boundary_conditions = self.lammpsbc(atoms).split()
288 if 's' in lammps_boundary_conditions:
289 pos = self.prism.vector_to_lammps(atoms.positions)
290 posmin = np.amin(pos, axis=0)
291 posmax = np.amax(pos, axis=0)
293 for i in range(3):
294 if lammps_boundary_conditions[i] == 's':
295 box_hi[i] = 1.05 * abs(posmax[i] - posmin[i])
297 cell_cmd = ('region cell prism '
298 '0 {} 0 {} 0 {} '
299 '{} {} {} units box'
300 ''.format(*box_hi, xy, xz, yz))
302 self.lmp.command(cell_cmd)
304 def set_lammps_pos(self, atoms: Atoms):
305 # Create local copy of positions that are wrapped along any periodic
306 # directions
307 pos = convert(atoms.positions, "distance", "ASE", self.units)
309 # wrap only after scaling and rotating to reduce chances of
310 # lammps neighbor list bugs.
311 pos = self.prism.vector_to_lammps(pos, wrap=True)
313 # Convert ase position matrix to lammps-style position array
314 # contiguous in memory
315 lmp_positions = list(pos.ravel())
317 # Convert that lammps-style array into a C object
318 c_double_array = (ctypes.c_double * len(lmp_positions))
319 lmp_c_positions = c_double_array(*lmp_positions)
320 # self.lmp.put_coosrds(lmp_c_positions)
321 self.lmp.scatter_atoms('x', 1, 3, lmp_c_positions)
323 def calculate(self, atoms, properties, system_changes):
324 self.propagate(atoms, properties, system_changes, 0)
326 def propagate(self, atoms, properties, system_changes, n_steps, dt=None,
327 dt_not_real_time=False, velocity_field=None):
328 """"atoms: Atoms object
329 Contains positions, unit-cell, ...
330 properties: list of str
331 List of what needs to be calculated. Can be any combination
332 of 'energy', 'forces', 'stress', 'dipole', 'charges', 'magmom'
333 and 'magmoms'.
334 system_changes: list of str
335 List of what has changed since last calculation. Can be
336 any combination of these five: 'positions', 'numbers', 'cell',
337 'pbc', 'initial_charges' and 'initial_magmoms'.
338 """
339 if len(system_changes) == 0:
340 return
342 if not self.started:
343 self.start_lammps()
344 if not self.initialized:
345 self.initialise_lammps(atoms)
346 else: # still need to reset cell
347 # NOTE: The whole point of ``post_changebox_cmds`` is that they're
348 # executed after any call to LAMMPS' change_box command. Here, we
349 # rely on the fact that self.set_cell(), where we have currently
350 # placed the execution of ``post_changebox_cmds``, gets called
351 # after this initial change_box call.
353 # Apply only requested boundary condition changes. Note this needs
354 # to happen before the call to set_cell since 'change_box' will
355 # apply any shrink-wrapping *after* it's updated the cell
356 # dimensions
357 if 'pbc' in system_changes:
358 change_box_str = 'change_box all boundary {}'
359 change_box_cmd = change_box_str.format(self.lammpsbc(atoms))
360 self.lmp.command(change_box_cmd)
362 # Reset positions so that if they are crazy from last
363 # propagation, change_box (in set_cell()) won't hang.
364 # Could do this only after testing for crazy positions?
365 # Could also use scatter_atoms() to set values (requires
366 # MPI comm), or extra_atoms() to get pointers to local
367 # data structures to zero, but then we would have to be
368 # careful with parallelism.
369 self.lmp.command("set atom * x 0.0 y 0.0 z 0.0")
370 self.set_cell(atoms, change=True)
372 if self.parameters.atom_types is None:
373 raise NameError("atom_types are mandatory.")
375 do_rebuild = (
376 not np.array_equal(atoms.numbers, self.previous_atoms_numbers)
377 or any(_ in system_changes for _ in ('numbers', 'initial_charges'))
378 )
379 if not do_rebuild:
380 do_redo_atom_types = not np.array_equal(
381 atoms.numbers, self.previous_atoms_numbers)
382 else:
383 do_redo_atom_types = False
385 self.lmp.command('echo none') # don't echo the atom positions
386 if do_rebuild:
387 self.rebuild(atoms)
388 elif do_redo_atom_types:
389 self.redo_atom_types(atoms)
390 self.lmp.command('echo log') # switch back log
392 self.set_lammps_pos(atoms)
394 if self.parameters.amendments is not None:
395 for cmd in self.parameters.amendments:
396 self.lmp.command(cmd)
398 if n_steps > 0:
399 if velocity_field is None:
400 vel = convert(
401 atoms.get_velocities(),
402 "velocity",
403 "ASE",
404 self.units)
405 else:
406 # FIXME: Do we need to worry about converting to lammps units
407 # here?
408 vel = atoms.arrays[velocity_field]
410 # Transform the velocities to new coordinate system
411 vel = self.prism.vector_to_lammps(vel)
413 # Convert ase velocities matrix to lammps-style velocities array
414 lmp_velocities = list(vel.ravel())
416 # Convert that lammps-style array into a C object
417 c_double_array = (ctypes.c_double * len(lmp_velocities))
418 lmp_c_velocities = c_double_array(*lmp_velocities)
419 self.lmp.scatter_atoms('v', 1, 3, lmp_c_velocities)
421 # Run for 0 time to calculate
422 if dt is not None:
423 if dt_not_real_time:
424 self.lmp.command('timestep %.30f' % dt)
425 else:
426 self.lmp.command('timestep %.30f' %
427 convert(dt, "time", "ASE", self.units))
428 self.lmp.command('run %d' % n_steps)
430 if n_steps > 0:
431 # TODO this must be slower than native copy, but why is it broken?
432 pos = np.array(
433 [x for x in self.lmp.gather_atoms("x", 1, 3)]).reshape(-1, 3)
434 pos = self.prism.vector_to_ase(pos)
436 # Convert from LAMMPS units to ASE units
437 pos = convert(pos, "distance", self.units, "ASE")
439 atoms.set_positions(pos)
441 vel = np.array(
442 [v for v in self.lmp.gather_atoms("v", 1, 3)]).reshape(-1, 3)
443 vel = self.prism.vector_to_lammps(vel)
444 if velocity_field is None:
445 atoms.set_velocities(convert(vel, 'velocity', self.units,
446 'ASE'))
448 # Extract the forces and energy
449 self.results['energy'] = convert(
450 self.lmp.extract_variable('pe', None, 0),
451 "energy", self.units, "ASE"
452 )
453 self.results['free_energy'] = self.results['energy']
455 # check for MPI active as per https://matsci.org/t/lammps-ids-in-python/60509/5
456 world_size = self.lmp.extract_setting('world_size')
457 if world_size != 1:
458 raise RuntimeError('Unsupported MPI active as indicated by '
459 f'world_size == {world_size} != 1')
460 # select just n_local which we assume is equal to len(atoms)
461 ids = self.lmp.numpy.extract_atom("id")[0:len(atoms)]
462 self.results["energies"] = convert(
463 self.lmp.numpy.extract_compute('pe_peratom',
464 self.LMP_STYLE_ATOM,
465 self.LMP_TYPE_VECTOR),
466 "energy", self.units, "ASE"
467 )
468 self.results["energies"][ids - 1] = self.results["energies"]
470 stress = np.empty(6)
471 stress_vars = ['pxx', 'pyy', 'pzz', 'pyz', 'pxz', 'pxy']
473 for i, var in enumerate(stress_vars):
474 stress[i] = self.lmp.extract_variable(var, None, 0)
476 stress_mat = np.zeros((3, 3))
477 stress_mat[0, 0] = stress[0]
478 stress_mat[1, 1] = stress[1]
479 stress_mat[2, 2] = stress[2]
480 stress_mat[1, 2] = stress[3]
481 stress_mat[2, 1] = stress[3]
482 stress_mat[0, 2] = stress[4]
483 stress_mat[2, 0] = stress[4]
484 stress_mat[0, 1] = stress[5]
485 stress_mat[1, 0] = stress[5]
487 stress_mat = self.prism.tensor2_to_ase(stress_mat)
489 stress[0] = stress_mat[0, 0]
490 stress[1] = stress_mat[1, 1]
491 stress[2] = stress_mat[2, 2]
492 stress[3] = stress_mat[1, 2]
493 stress[4] = stress_mat[0, 2]
494 stress[5] = stress_mat[0, 1]
496 self.results['stress'] = convert(-stress, "pressure", self.units, "ASE")
498 # definitely yields atom-id ordered force array
499 f = convert(np.array(self.lmp.gather_atoms("f", 1, 3)).reshape(-1, 3),
500 "force", self.units, "ASE")
501 self.results['forces'] = self.prism.vector_to_ase(f)
503 # otherwise check_state will always trigger a new calculation
504 self.atoms = atoms.copy()
506 if not self.parameters["keep_alive"]:
507 self.clean()
509 def lammpsbc(self, atoms):
510 """Determine LAMMPS boundary types based on ASE pbc settings. For
511 non-periodic dimensions, if the cell length is finite then
512 fixed BCs ('f') are used; if the cell length is approximately
513 zero, shrink-wrapped BCs ('s') are used."""
515 retval = ''
516 pbc = atoms.get_pbc()
517 if np.all(pbc):
518 retval = 'p p p'
519 else:
520 cell = atoms.get_cell()
521 for i in range(3):
522 if pbc[i]:
523 retval += 'p '
524 else:
525 # See if we're using indefinite ASE boundaries along this
526 # direction
527 if np.linalg.norm(cell[i]) < np.finfo(cell[i][0]).tiny:
528 retval += 's '
529 else:
530 retval += 'f '
532 return retval.strip()
534 def rebuild(self, atoms):
535 try:
536 n_diff = len(atoms.numbers) - len(self.previous_atoms_numbers)
537 except Exception: # XXX Which kind of exception?
538 n_diff = len(atoms.numbers)
540 if n_diff > 0:
541 if any(("reax/c" in cmd) for cmd in self.parameters.lmpcmds):
542 self.lmp.command("pair_style lj/cut 2.5")
543 self.lmp.command("pair_coeff * * 1 1")
545 for cmd in self.parameters.lmpcmds:
546 if (("pair_style" in cmd) or ("pair_coeff" in cmd) or
547 ("qeq/reax" in cmd)):
548 self.lmp.command(cmd)
550 cmd = f"create_atoms 1 random {n_diff} 1 NULL"
551 self.lmp.command(cmd)
552 elif n_diff < 0:
553 cmd = "group delatoms id {}:{}".format(
554 len(atoms.numbers) + 1, len(self.previous_atoms_numbers))
555 self.lmp.command(cmd)
556 cmd = "delete_atoms group delatoms"
557 self.lmp.command(cmd)
559 self.redo_atom_types(atoms)
561 def redo_atom_types(self, atoms):
562 current_types = {
563 (i + 1, self.parameters.atom_types[sym]) for i, sym
564 in enumerate(atoms.get_chemical_symbols())}
566 try:
567 previous_types = {
568 (i + 1, self.parameters.atom_types[ase_chemical_symbols[Z]])
569 for i, Z in enumerate(self.previous_atoms_numbers)}
570 except Exception: # XXX which kind of exception?
571 previous_types = set()
573 for (i, i_type) in current_types - previous_types:
574 cmd = f"set atom {i} type {i_type}"
575 self.lmp.command(cmd)
577 # set charges only when LAMMPS `atom_style` permits charges
578 # https://docs.lammps.org/Library_properties.html#extract-atom-flags
579 if self.lmp.extract_setting('q_flag') == 1:
580 charges = atoms.get_initial_charges()
581 if np.any(charges != 0.0):
582 for i, q in enumerate(charges):
583 self.lmp.command(f'set atom {i + 1} charge {q}')
585 self.previous_atoms_numbers = atoms.numbers.copy()
587 def restart_lammps(self, atoms):
588 if self.started:
589 self.lmp.command("clear")
590 # hope there's no other state to be reset
591 self.started = False
592 self.initialized = False
593 self.previous_atoms_numbers = []
594 self.start_lammps()
595 self.initialise_lammps(atoms)
597 def start_lammps(self):
598 # Only import lammps when running a calculation
599 # so it is not required to use other parts of the
600 # module
601 from lammps import LMP_STYLE_ATOM, LMP_TYPE_VECTOR, lammps
603 self.LMP_STYLE_ATOM = LMP_STYLE_ATOM
604 self.LMP_TYPE_VECTOR = LMP_TYPE_VECTOR
606 # start lammps process
607 if self.parameters.log_file is None:
608 cmd_args = ['-echo', 'log', '-log', 'none', '-screen', 'none',
609 '-nocite']
610 else:
611 cmd_args = ['-echo', 'log', '-log', self.parameters.log_file,
612 '-screen', 'none', '-nocite']
614 self.cmd_args = cmd_args + list(self.parameters.extra_cmd_args)
616 if self.lmp is None:
617 self.lmp = lammps(self.parameters.lammps_name, self.cmd_args,
618 comm=self.parameters.comm)
620 if self.parameters.initializer is not None:
621 self.parameters.initializer(self.lmp)
623 # Run header commands to set up lammps (units, etc.)
624 for cmd in self.parameters.lammps_header:
625 self.lmp.command(cmd)
627 for cmd in self.parameters.lammps_header:
628 if "units" in cmd:
629 self.units = cmd.split()[1]
631 if 'lammps_header_extra' in self.parameters:
632 if self.parameters.lammps_header_extra is not None:
633 for cmd in self.parameters.lammps_header_extra:
634 self.lmp.command(cmd)
636 self.started = True
638 def initialise_lammps(self, atoms):
639 # Initialising commands
640 if self.parameters.boundary:
641 # if the boundary command is in the supplied commands use that
642 # otherwise use atoms pbc
643 for cmd in self.parameters.lmpcmds:
644 if 'boundary' in cmd:
645 break
646 else:
647 self.lmp.command('boundary ' + self.lammpsbc(atoms))
649 # Initialize cell
650 self.set_cell(atoms, change=not self.parameters.create_box)
652 if self.parameters.atom_types is None:
653 # if None is given, create from atoms object in order of appearance
654 s = atoms.get_chemical_symbols()
655 _, idx = np.unique(s, return_index=True)
656 s_red = np.array(s)[np.sort(idx)].tolist()
657 self.parameters.atom_types = {j: i + 1 for i, j in enumerate(s_red)}
659 # Initialize box
660 if self.parameters.create_box:
661 # count number of known types
662 n_types = len(self.parameters.atom_types)
663 create_box_command = f'create_box {n_types} cell'
664 self.lmp.command(create_box_command)
666 # Initialize the atoms with their types
667 # positions do not matter here
668 if self.parameters.create_atoms:
669 self.lmp.command('echo none') # don't echo the atom positions
670 self.rebuild(atoms)
671 self.lmp.command('echo log') # turn back on
672 else:
673 self.previous_atoms_numbers = atoms.numbers.copy()
675 # execute the user commands
676 for cmd in self.parameters.lmpcmds + ["compute pe_peratom all pe/atom"]:
677 self.lmp.command(cmd)
679 # Set masses after user commands, e.g. to override
680 # EAM-provided masses
681 for sym in self.parameters.atom_types:
682 if self.parameters.atom_type_masses is None:
683 mass = ase_atomic_masses[ase_atomic_numbers[sym]]
684 else:
685 mass = self.parameters.atom_type_masses[sym]
686 self.lmp.command('mass %d %.30f' % (
687 self.parameters.atom_types[sym],
688 convert(mass, "mass", "ASE", self.units)))
690 # Define force & energy variables for extraction
691 self.lmp.command('variable pxx equal pxx')
692 self.lmp.command('variable pyy equal pyy')
693 self.lmp.command('variable pzz equal pzz')
694 self.lmp.command('variable pxy equal pxy')
695 self.lmp.command('variable pxz equal pxz')
696 self.lmp.command('variable pyz equal pyz')
698 # I am not sure why we need this next line but LAMMPS will
699 # raise an error if it is not there. Perhaps it is needed to
700 # ensure the cell stresses are calculated
701 self.lmp.command('thermo_style custom pe pxx emol ecoul')
703 self.lmp.command('variable fx atom fx')
704 self.lmp.command('variable fy atom fy')
705 self.lmp.command('variable fz atom fz')
707 # do we need this if we extract from a global ?
708 self.lmp.command('variable pe equal pe')
710 self.lmp.command("neigh_modify delay 0 every 1 check yes")
712 self.initialized = True