Coverage for /builds/ase/ase/ase/calculators/lammpslib.py : 77.97%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""ASE LAMMPS Calculator Library Version"""
4import ctypes
6import numpy as np
7from numpy.linalg import norm
9from ase.calculators.calculator import Calculator
10from ase.data import (atomic_numbers as ase_atomic_numbers,
11 chemical_symbols as ase_chemical_symbols,
12 atomic_masses as ase_atomic_masses)
13from ase.calculators.lammps import convert
14from ase.geometry import wrap_positions
16# TODO
17# 1. should we make a new lammps object each time ?
18# 4. need a routine to get the model back from lammps
19# 5. if we send a command to lmps directly then the calculator does
20# not know about it and the energy could be wrong.
21# 6. do we need a subroutine generator that converts a lammps string
22# into a python function that can be called
23# 8. make matscipy as fallback
24# 9. keep_alive not needed with no system changes
25# 10. it may be a good idea to unify the cell handling with the one found in
26# lammpsrun.py
29# this one may be moved to some more generic place
30def is_upper_triangular(arr, atol=1e-8):
31 """test for upper triangular matrix based on numpy"""
32 # must be (n x n) matrix
33 assert len(arr.shape) == 2
34 assert arr.shape[0] == arr.shape[1]
35 return np.allclose(np.tril(arr, k=-1), 0., atol=atol) and \
36 np.all(np.diag(arr) >= 0.0)
39def convert_cell(ase_cell):
40 """
41 Convert a parallelepiped (forming right hand basis)
42 to lower triangular matrix LAMMPS can accept. This
43 function transposes cell matrix so the bases are column vectors
44 """
45 cell = ase_cell.T
47 if not is_upper_triangular(cell):
48 # rotate bases into triangular matrix
49 tri_mat = np.zeros((3, 3))
50 A = cell[:, 0]
51 B = cell[:, 1]
52 C = cell[:, 2]
53 tri_mat[0, 0] = norm(A)
54 Ahat = A / norm(A)
55 AxBhat = np.cross(A, B) / norm(np.cross(A, B))
56 tri_mat[0, 1] = np.dot(B, Ahat)
57 tri_mat[1, 1] = norm(np.cross(Ahat, B))
58 tri_mat[0, 2] = np.dot(C, Ahat)
59 tri_mat[1, 2] = np.dot(C, np.cross(AxBhat, Ahat))
60 tri_mat[2, 2] = norm(np.dot(C, AxBhat))
62 # create and save the transformation for coordinates
63 volume = np.linalg.det(ase_cell)
64 trans = np.array([np.cross(B, C), np.cross(C, A), np.cross(A, B)])
65 trans /= volume
66 coord_transform = np.dot(tri_mat, trans)
68 return tri_mat, coord_transform
69 else:
70 return cell, None
73class LAMMPSlib(Calculator):
74 r"""
75**Introduction**
77LAMMPSlib is an interface and calculator for LAMMPS_. LAMMPSlib uses
78the python interface that comes with LAMMPS to solve an atoms model
79for energy, atom forces and cell stress. This calculator creates a
80'.lmp' object which is a running lammps program, so further commands
81can be sent to this object executed until it is explicitly closed. Any
82additional variables calculated by lammps can also be extracted. This
83is still experimental code.
85**Arguments**
87======================= ======================================================
88Keyword Description
89======================= ======================================================
90``lmpcmds`` list of strings of LAMMPS commands. You need to supply
91 enough to define the potential to be used e.g.
93 ["pair_style eam/alloy",
94 "pair_coeff * * potentials/NiAlH_jea.eam.alloy Ni Al"]
96``atom_types`` dictionary of ``atomic_symbol :lammps_atom_type``
97 pairs, e.g. ``{'Cu':1}`` to bind copper to lammps
98 atom type 1. If <None>, autocreated by assigning
99 lammps atom types in order that they appear in the
100 first used atoms object.
102``atom_type_masses`` dictionary of ``atomic_symbol :mass`` pairs, e.g.
103 ``{'Cu':63.546}`` to optionally assign masses that
104 override default ase.data.atomic_masses. Note that
105 since unit conversion is done automatically in this
106 module, these quantities must be given in the
107 standard ase mass units (g/mol)
109``log_file`` string
110 path to the desired LAMMPS log file
112``lammps_header`` string to use for lammps setup. Default is to use
113 metal units and simple atom simulation.
115 lammps_header=['units metal',
116 'atom_style atomic',
117 'atom_modify map array sort 0 0'])
119``amendments`` extra list of strings of LAMMPS commands to be run
120 post initialization. (Use: Initialization amendments)
121 e.g.
123 ["mass 1 58.6934"]
125``post_changebox_cmds`` extra list of strings of LAMMPS commands to be run
126 after any LAMMPS 'change_box' command is performed by
127 the calculator. This is relevant because some
128 potentials either themselves depend on the geometry
129 and boundary conditions of the simulation box, or are
130 frequently coupled with other LAMMPS commands that
131 do, e.g. the 'buck/coul/long' pair style is often
132 used with the kspace_* commands, which are sensitive
133 to the periodicity of the simulation box.
135``keep_alive`` Boolean
136 whether to keep the lammps routine alive for more
137 commands
139======================= ======================================================
142**Requirements**
144To run this calculator you must have LAMMPS installed and compiled to
145enable the python interface. See the LAMMPS manual.
147If the following code runs then lammps is installed correctly.
149 >>> from lammps import lammps
150 >>> lmp = lammps()
152The version of LAMMPS is also important. LAMMPSlib is suitable for
153versions after approximately 2011. Prior to this the python interface
154is slightly different from that used by LAMMPSlib. It is not difficult
155to change to the earlier format.
157**LAMMPS and LAMMPSlib**
159The LAMMPS calculator is another calculator that uses LAMMPS (the
160program) to calculate the energy by generating input files and running
161a separate LAMMPS job to perform the analysis. The output data is then
162read back into python. LAMMPSlib makes direct use of the LAMMPS (the
163program) python interface. As well as directly running any LAMMPS
164command line it allows the values of any of LAMMPS variables to be
165extracted and returned to python.
167**Example**
169Provided that the respective potential file is in the working directory, one
170can simply run (note that LAMMPS needs to be compiled to work with EAM
171potentials)
173::
175 from ase import Atom, Atoms
176 from ase.build import bulk
177 from ase.calculators.lammpslib import LAMMPSlib
179 cmds = ["pair_style eam/alloy",
180 "pair_coeff * * NiAlH_jea.eam.alloy Ni H"]
182 Ni = bulk('Ni', cubic=True)
183 H = Atom('H', position=Ni.cell.diagonal()/2)
184 NiH = Ni + H
186 lammps = LAMMPSlib(lmpcmds=cmds, log_file='test.log')
188 NiH.calc = lammps
189 print("Energy ", NiH.get_potential_energy())
192**Implementation**
194LAMMPS provides a set of python functions to allow execution of the
195underlying C++ LAMMPS code. The functions used by the LAMMPSlib
196interface are::
198 from lammps import lammps
200 lmp = lammps(cmd_args) # initiate LAMMPS object with command line args
202 lmp.scatter_atoms('x',1,3,positions) # atom coords to LAMMPS C array
203 lmp.command(cmd) # executes a one line cmd string
204 lmp.extract_variable(...) # extracts a per atom variable
205 lmp.extract_global(...) # extracts a global variable
206 lmp.close() # close the lammps object
208For a single Ni atom model the following lammps file commands would be run
209by invoking the get_potential_energy() method::
211 units metal
212 atom_style atomic
213 atom_modify map array sort 0 0
215 region cell prism 0 xhi 0 yhi 0 zhi xy xz yz units box
216 create_box 1 cell
217 create_atoms 1 single 0 0 0 units box
218 mass * 1.0
220 ## user lmpcmds get executed here
221 pair_style eam/alloy
222 pair_coeff * * NiAlH_jea.eam.alloy Ni
223 ## end of user lmmpcmds
225 run 0
227where xhi, yhi and zhi are the lattice vector lengths and xy,
228xz and yz are the tilt of the lattice vectors, all to be edited.
231**Notes**
233.. _LAMMPS: http://lammps.sandia.gov/
235* Units: The default lammps_header sets the units to Angstrom and eV
236 and for compatibility with ASE Stress is in GPa.
238* The global energy is currently extracted from LAMMPS using
239 extract_variable since lammps.lammps currently extract_global only
240 accepts the following ['dt', 'boxxlo', 'boxxhi', 'boxylo', 'boxyhi',
241 'boxzlo', 'boxzhi', 'natoms', 'nlocal'].
243* If an error occurs while lammps is in control it will crash
244 Python. Check the output of the log file to find the lammps error.
246* If the are commands directly sent to the LAMMPS object this may
247 change the energy value of the model. However the calculator will not
248 know of it and still return the original energy value.
250"""
252 implemented_properties = ['energy', 'free_energy', 'forces', 'stress',
253 'energies']
255 started = False
256 initialized = False
258 default_parameters = dict(
259 atom_types=None,
260 atom_type_masses=None,
261 log_file=None,
262 lammps_name='',
263 keep_alive=False,
264 lammps_header=['units metal',
265 'atom_style atomic',
266 'atom_modify map array sort 0 0'],
267 amendments=None,
268 post_changebox_cmds=None,
269 boundary=True,
270 create_box=True,
271 create_atoms=True,
272 read_molecular_info=False,
273 comm=None)
275 def __init__(self, *args, **kwargs):
276 Calculator.__init__(self, *args, **kwargs)
277 self.lmp = None
279 def __del__(self):
280 if self.started:
281 self.lmp.close()
282 self.started = False
283 self.lmp = None
285 def set_cell(self, atoms, change=False):
286 lammps_cell, self.coord_transform = convert_cell(atoms.get_cell())
288 xhi, xy, xz, _, yhi, yz, _, _, zhi = convert(
289 lammps_cell.flatten(order='C'), "distance", "ASE", self.units)
290 box_hi = [xhi, yhi, zhi]
292 if change:
293 cell_cmd = ('change_box all '
294 'x final 0 {} y final 0 {} z final 0 {} '
295 'xy final {} xz final {} yz final {} units box'
296 ''.format(xhi, yhi, zhi, xy, xz, yz))
297 if self.parameters.post_changebox_cmds is not None:
298 for cmd in self.parameters.post_changebox_cmds:
299 self.lmp.command(cmd)
300 else:
301 # just in case we'll want to run with a funny shape box,
302 # and here command will only happen once, and before
303 # any calculation
304 if self.parameters.create_box:
305 self.lmp.command('box tilt large')
307 # Check if there are any indefinite boundaries. If so,
308 # shrink-wrapping will end up being used, but we want to
309 # define the LAMMPS region and box fairly tight around the
310 # atoms to avoid losing any
311 lammps_boundary_conditions = self.lammpsbc(atoms).split()
312 if 's' in lammps_boundary_conditions:
313 pos = atoms.get_positions()
314 if self.coord_transform is not None:
315 pos = np.dot(self.coord_transform, pos.transpose())
316 pos = pos.transpose()
317 posmin = np.amin(pos, axis=0)
318 posmax = np.amax(pos, axis=0)
320 for i in range(0, 3):
321 if lammps_boundary_conditions[i] == 's':
322 box_hi[i] = 1.05 * abs(posmax[i] - posmin[i])
324 cell_cmd = ('region cell prism '
325 '0 {} 0 {} 0 {} '
326 '{} {} {} units box'
327 ''.format(*box_hi, xy, xz, yz))
329 self.lmp.command(cell_cmd)
331 def set_lammps_pos(self, atoms):
332 # Create local copy of positions that are wrapped along any periodic
333 # directions
334 cell = convert(atoms.cell, "distance", "ASE", self.units)
335 pos = convert(atoms.positions, "distance", "ASE", self.units)
337 # If necessary, transform the positions to new coordinate system
338 if self.coord_transform is not None:
339 pos = np.dot(pos, self.coord_transform.T)
340 cell = np.dot(cell, self.coord_transform.T)
342 # wrap only after scaling and rotating to reduce chances of
343 # lammps neighbor list bugs.
344 pos = wrap_positions(pos, cell, atoms.get_pbc())
346 # Convert ase position matrix to lammps-style position array
347 # contiguous in memory
348 lmp_positions = list(pos.ravel())
350 # Convert that lammps-style array into a C object
351 c_double_array = (ctypes.c_double * len(lmp_positions))
352 lmp_c_positions = c_double_array(*lmp_positions)
353 # self.lmp.put_coosrds(lmp_c_positions)
354 self.lmp.scatter_atoms('x', 1, 3, lmp_c_positions)
356 def calculate(self, atoms, properties, system_changes):
357 self.propagate(atoms, properties, system_changes, 0)
359 def propagate(self, atoms, properties, system_changes, n_steps, dt=None,
360 dt_not_real_time=False, velocity_field=None):
361 """"atoms: Atoms object
362 Contains positions, unit-cell, ...
363 properties: list of str
364 List of what needs to be calculated. Can be any combination
365 of 'energy', 'forces', 'stress', 'dipole', 'charges', 'magmom'
366 and 'magmoms'.
367 system_changes: list of str
368 List of what has changed since last calculation. Can be
369 any combination of these five: 'positions', 'numbers', 'cell',
370 'pbc', 'charges' and 'magmoms'.
371 """
372 if len(system_changes) == 0:
373 return
375 self.coord_transform = None
377 if not self.started:
378 self.start_lammps()
379 if not self.initialized:
380 self.initialise_lammps(atoms)
381 else: # still need to reset cell
382 # NOTE: The whole point of ``post_changebox_cmds`` is that they're
383 # executed after any call to LAMMPS' change_box command. Here, we
384 # rely on the fact that self.set_cell(), where we have currently
385 # placed the execution of ``post_changebox_cmds``, gets called
386 # after this initial change_box call.
388 # Apply only requested boundary condition changes. Note this needs
389 # to happen before the call to set_cell since 'change_box' will
390 # apply any shrink-wrapping *after* it's updated the cell
391 # dimensions
392 if 'pbc' in system_changes:
393 change_box_str = 'change_box all boundary {}'
394 change_box_cmd = change_box_str.format(self.lammpsbc(atoms))
395 self.lmp.command(change_box_cmd)
397 # Reset positions so that if they are crazy from last
398 # propagation, change_box (in set_cell()) won't hang.
399 # Could do this only after testing for crazy positions?
400 # Could also use scatter_atoms() to set values (requires
401 # MPI comm), or extra_atoms() to get pointers to local
402 # data structures to zero, but then we would have to be
403 # careful with parallelism.
404 self.lmp.command("set atom * x 0.0 y 0.0 z 0.0")
405 self.set_cell(atoms, change=True)
407 if self.parameters.atom_types is None:
408 raise NameError("atom_types are mandatory.")
410 do_rebuild = (not np.array_equal(atoms.numbers,
411 self.previous_atoms_numbers)
412 or ("numbers" in system_changes))
413 if not do_rebuild:
414 do_redo_atom_types = not np.array_equal(
415 atoms.numbers, self.previous_atoms_numbers)
416 else:
417 do_redo_atom_types = False
419 self.lmp.command('echo none') # don't echo the atom positions
420 if do_rebuild:
421 self.rebuild(atoms)
422 elif do_redo_atom_types:
423 self.redo_atom_types(atoms)
424 self.lmp.command('echo log') # switch back log
426 self.set_lammps_pos(atoms)
428 if self.parameters.amendments is not None:
429 for cmd in self.parameters.amendments:
430 self.lmp.command(cmd)
432 if n_steps > 0:
433 if velocity_field is None:
434 vel = convert(
435 atoms.get_velocities(),
436 "velocity",
437 "ASE",
438 self.units)
439 else:
440 # FIXME: Do we need to worry about converting to lammps units
441 # here?
442 vel = atoms.arrays[velocity_field]
444 # If necessary, transform the velocities to new coordinate system
445 if self.coord_transform is not None:
446 vel = np.dot(self.coord_transform, vel.T).T
448 # Convert ase velocities matrix to lammps-style velocities array
449 lmp_velocities = list(vel.ravel())
451 # Convert that lammps-style array into a C object
452 c_double_array = (ctypes.c_double * len(lmp_velocities))
453 lmp_c_velocities = c_double_array(*lmp_velocities)
454 self.lmp.scatter_atoms('v', 1, 3, lmp_c_velocities)
456 # Run for 0 time to calculate
457 if dt is not None:
458 if dt_not_real_time:
459 self.lmp.command('timestep %.30f' % dt)
460 else:
461 self.lmp.command('timestep %.30f' %
462 convert(dt, "time", "ASE", self.units))
463 self.lmp.command('run %d' % n_steps)
465 if n_steps > 0:
466 # TODO this must be slower than native copy, but why is it broken?
467 pos = np.array(
468 [x for x in self.lmp.gather_atoms("x", 1, 3)]).reshape(-1, 3)
469 if self.coord_transform is not None:
470 pos = np.dot(pos, self.coord_transform)
472 # Convert from LAMMPS units to ASE units
473 pos = convert(pos, "distance", self.units, "ASE")
475 atoms.set_positions(pos)
477 vel = np.array(
478 [v for v in self.lmp.gather_atoms("v", 1, 3)]).reshape(-1, 3)
479 if self.coord_transform is not None:
480 vel = np.dot(vel, self.coord_transform)
481 if velocity_field is None:
482 atoms.set_velocities(convert(vel, 'velocity', self.units,
483 'ASE'))
485 # Extract the forces and energy
486 self.results['energy'] = convert(
487 self.lmp.extract_variable('pe', None, 0),
488 "energy", self.units, "ASE"
489 )
490 self.results['free_energy'] = self.results['energy']
492 ids = self.lmp.numpy.extract_atom("id")
493 # if ids doesn't match atoms then data is MPI distributed, which
494 # we can't handle
495 assert len(ids) == len(atoms)
496 self.results["energies"] = convert(
497 self.lmp.numpy.extract_compute('pe_peratom',
498 self.LMP_STYLE_ATOM,
499 self.LMP_TYPE_VECTOR),
500 "energy", self.units, "ASE"
501 )
502 self.results["energies"][ids - 1] = self.results["energies"]
504 stress = np.empty(6)
505 stress_vars = ['pxx', 'pyy', 'pzz', 'pyz', 'pxz', 'pxy']
507 for i, var in enumerate(stress_vars):
508 stress[i] = self.lmp.extract_variable(var, None, 0)
510 stress_mat = np.zeros((3, 3))
511 stress_mat[0, 0] = stress[0]
512 stress_mat[1, 1] = stress[1]
513 stress_mat[2, 2] = stress[2]
514 stress_mat[1, 2] = stress[3]
515 stress_mat[2, 1] = stress[3]
516 stress_mat[0, 2] = stress[4]
517 stress_mat[2, 0] = stress[4]
518 stress_mat[0, 1] = stress[5]
519 stress_mat[1, 0] = stress[5]
520 if self.coord_transform is not None:
521 stress_mat = np.dot(self.coord_transform.T,
522 np.dot(stress_mat, self.coord_transform))
523 stress[0] = stress_mat[0, 0]
524 stress[1] = stress_mat[1, 1]
525 stress[2] = stress_mat[2, 2]
526 stress[3] = stress_mat[1, 2]
527 stress[4] = stress_mat[0, 2]
528 stress[5] = stress_mat[0, 1]
530 self.results['stress'] = convert(-stress, "pressure", self.units, "ASE")
532 # definitely yields atom-id ordered force array
533 f = convert(np.array(self.lmp.gather_atoms("f", 1, 3)).reshape(-1, 3),
534 "force", self.units, "ASE")
536 if self.coord_transform is not None:
537 self.results['forces'] = np.dot(f, self.coord_transform)
538 else:
539 self.results['forces'] = f.copy()
541 # otherwise check_state will always trigger a new calculation
542 self.atoms = atoms.copy()
544 if not self.parameters.keep_alive:
545 self.lmp.close()
547 def lammpsbc(self, atoms):
548 """Determine LAMMPS boundary types based on ASE pbc settings. For
549 non-periodic dimensions, if the cell length is finite then
550 fixed BCs ('f') are used; if the cell length is approximately
551 zero, shrink-wrapped BCs ('s') are used."""
553 retval = ''
554 pbc = atoms.get_pbc()
555 if np.all(pbc):
556 retval = 'p p p'
557 else:
558 cell = atoms.get_cell()
559 for i in range(0, 3):
560 if pbc[i]:
561 retval += 'p '
562 else:
563 # See if we're using indefinite ASE boundaries along this
564 # direction
565 if np.linalg.norm(cell[i]) < np.finfo(cell[i][0]).tiny:
566 retval += 's '
567 else:
568 retval += 'f '
570 return retval.strip()
572 def rebuild(self, atoms):
573 try:
574 n_diff = len(atoms.numbers) - len(self.previous_atoms_numbers)
575 except Exception: # XXX Which kind of exception?
576 n_diff = len(atoms.numbers)
578 if n_diff > 0:
579 if any([("reax/c" in cmd) for cmd in self.parameters.lmpcmds]):
580 self.lmp.command("pair_style lj/cut 2.5")
581 self.lmp.command("pair_coeff * * 1 1")
583 for cmd in self.parameters.lmpcmds:
584 if (("pair_style" in cmd) or ("pair_coeff" in cmd) or
585 ("qeq/reax" in cmd)):
586 self.lmp.command(cmd)
588 cmd = "create_atoms 1 random {} 1 NULL".format(n_diff)
589 self.lmp.command(cmd)
590 elif n_diff < 0:
591 cmd = "group delatoms id {}:{}".format(
592 len(atoms.numbers) + 1, len(self.previous_atoms_numbers))
593 self.lmp.command(cmd)
594 cmd = "delete_atoms group delatoms"
595 self.lmp.command(cmd)
597 self.redo_atom_types(atoms)
599 def redo_atom_types(self, atoms):
600 current_types = set(
601 (i + 1, self.parameters.atom_types[sym]) for i, sym
602 in enumerate(atoms.get_chemical_symbols()))
604 try:
605 previous_types = set(
606 (i + 1, self.parameters.atom_types[ase_chemical_symbols[Z]])
607 for i, Z in enumerate(self.previous_atoms_numbers))
608 except Exception: # XXX which kind of exception?
609 previous_types = set()
611 for (i, i_type) in current_types - previous_types:
612 cmd = "set atom {} type {}".format(i, i_type)
613 self.lmp.command(cmd)
615 self.previous_atoms_numbers = atoms.numbers.copy()
617 def restart_lammps(self, atoms):
618 if self.started:
619 self.lmp.command("clear")
620 # hope there's no other state to be reset
621 self.started = False
622 self.initialized = False
623 self.previous_atoms_numbers = []
624 self.start_lammps()
625 self.initialise_lammps(atoms)
627 def start_lammps(self):
628 # Only import lammps when running a calculation
629 # so it is not required to use other parts of the
630 # module
631 from lammps import lammps, LMP_STYLE_ATOM, LMP_TYPE_VECTOR
633 self.LMP_STYLE_ATOM = LMP_STYLE_ATOM
634 self.LMP_TYPE_VECTOR = LMP_TYPE_VECTOR
636 # start lammps process
637 if self.parameters.log_file is None:
638 cmd_args = ['-echo', 'log', '-log', 'none', '-screen', 'none',
639 '-nocite']
640 else:
641 cmd_args = ['-echo', 'log', '-log', self.parameters.log_file,
642 '-screen', 'none', '-nocite']
644 self.cmd_args = cmd_args
646 if self.lmp is None:
647 self.lmp = lammps(self.parameters.lammps_name, self.cmd_args,
648 comm=self.parameters.comm)
650 # Run header commands to set up lammps (units, etc.)
651 for cmd in self.parameters.lammps_header:
652 self.lmp.command(cmd)
654 for cmd in self.parameters.lammps_header:
655 if "units" in cmd:
656 self.units = cmd.split()[1]
658 if 'lammps_header_extra' in self.parameters:
659 if self.parameters.lammps_header_extra is not None:
660 for cmd in self.parameters.lammps_header_extra:
661 self.lmp.command(cmd)
663 self.started = True
665 def initialise_lammps(self, atoms):
666 # Initialising commands
667 if self.parameters.boundary:
668 # if the boundary command is in the supplied commands use that
669 # otherwise use atoms pbc
670 for cmd in self.parameters.lmpcmds:
671 if 'boundary' in cmd:
672 break
673 else:
674 self.lmp.command('boundary ' + self.lammpsbc(atoms))
676 # Initialize cell
677 self.set_cell(atoms, change=not self.parameters.create_box)
679 if self.parameters.atom_types is None:
680 # if None is given, create from atoms object in order of appearance
681 s = atoms.get_chemical_symbols()
682 _, idx = np.unique(s, return_index=True)
683 s_red = np.array(s)[np.sort(idx)].tolist()
684 self.parameters.atom_types = {j: i + 1 for i, j in enumerate(s_red)}
686 # Initialize box
687 if self.parameters.create_box:
688 # count number of known types
689 n_types = len(self.parameters.atom_types)
690 create_box_command = 'create_box {} cell'.format(n_types)
691 self.lmp.command(create_box_command)
693 # Initialize the atoms with their types
694 # positions do not matter here
695 if self.parameters.create_atoms:
696 self.lmp.command('echo none') # don't echo the atom positions
697 self.rebuild(atoms)
698 self.lmp.command('echo log') # turn back on
699 else:
700 self.previous_atoms_numbers = atoms.numbers.copy()
702 # execute the user commands
703 for cmd in self.parameters.lmpcmds + ["compute pe_peratom all pe/atom"]:
704 self.lmp.command(cmd)
706 # Set masses after user commands, e.g. to override
707 # EAM-provided masses
708 for sym in self.parameters.atom_types:
709 if self.parameters.atom_type_masses is None:
710 mass = ase_atomic_masses[ase_atomic_numbers[sym]]
711 else:
712 mass = self.parameters.atom_type_masses[sym]
713 self.lmp.command('mass %d %.30f' % (
714 self.parameters.atom_types[sym],
715 convert(mass, "mass", "ASE", self.units)))
717 # Define force & energy variables for extraction
718 self.lmp.command('variable pxx equal pxx')
719 self.lmp.command('variable pyy equal pyy')
720 self.lmp.command('variable pzz equal pzz')
721 self.lmp.command('variable pxy equal pxy')
722 self.lmp.command('variable pxz equal pxz')
723 self.lmp.command('variable pyz equal pyz')
725 # I am not sure why we need this next line but LAMMPS will
726 # raise an error if it is not there. Perhaps it is needed to
727 # ensure the cell stresses are calculated
728 self.lmp.command('thermo_style custom pe pxx emol ecoul')
730 self.lmp.command('variable fx atom fx')
731 self.lmp.command('variable fy atom fy')
732 self.lmp.command('variable fz atom fz')
734 # do we need this if we extract from a global ?
735 self.lmp.command('variable pe equal pe')
737 self.lmp.command("neigh_modify delay 0 every 1 check yes")
739 self.initialized = True