Coverage for /builds/debichem-team/python-ase/ase/calculators/lammpslib.py: 73.96%
288 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +0000
1"""ASE LAMMPS Calculator Library Version"""
4import ctypes
6import numpy as np
7from numpy.linalg import norm
9from ase import Atoms
10from ase.calculators.calculator import Calculator
11from ase.calculators.lammps import Prism, convert
12from ase.data import atomic_masses as ase_atomic_masses
13from ase.data import atomic_numbers as ase_atomic_numbers
14from ase.data import chemical_symbols as ase_chemical_symbols
15from ase.utils import deprecated
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
28# this one may be moved to some more generic place
29@deprecated("Please use the technique in https://stackoverflow.com/a/26912166")
30def is_upper_triangular(arr, atol=1e-8):
31 """test for upper triangular matrix based on numpy
32 .. deprecated:: 3.23.0
33 Please use the technique in https://stackoverflow.com/a/26912166
34 """
35 # must be (n x n) matrix
36 assert len(arr.shape) == 2
37 assert arr.shape[0] == arr.shape[1]
38 return np.allclose(np.tril(arr, k=-1), 0., atol=atol) and \
39 np.all(np.diag(arr) >= 0.0)
42@deprecated(
43 "Please use "
44 "`ase.calculators.lammps.coordinatetransform.calc_rotated_cell`. "
45 "Note that the new function returns the ASE lower trianglar cell and does "
46 "not return the conversion matrix."
47)
48def convert_cell(ase_cell):
49 """
50 Convert a parallelepiped (forming right hand basis)
51 to lower triangular matrix LAMMPS can accept. This
52 function transposes cell matrix so the bases are column vectors
54 .. deprecated:: 3.23.0
55 Please use
56 :func:`~ase.calculators.lammps.coordinatetransform.calc_rotated_cell`.
57 """
58 cell = ase_cell.T
60 if not is_upper_triangular(cell):
61 # rotate bases into triangular matrix
62 tri_mat = np.zeros((3, 3))
63 A = cell[:, 0]
64 B = cell[:, 1]
65 C = cell[:, 2]
66 tri_mat[0, 0] = norm(A)
67 Ahat = A / norm(A)
68 AxBhat = np.cross(A, B) / norm(np.cross(A, B))
69 tri_mat[0, 1] = np.dot(B, Ahat)
70 tri_mat[1, 1] = norm(np.cross(Ahat, B))
71 tri_mat[0, 2] = np.dot(C, Ahat)
72 tri_mat[1, 2] = np.dot(C, np.cross(AxBhat, Ahat))
73 tri_mat[2, 2] = norm(np.dot(C, AxBhat))
75 # create and save the transformation for coordinates
76 volume = np.linalg.det(ase_cell)
77 trans = np.array([np.cross(B, C), np.cross(C, A), np.cross(A, B)])
78 trans /= volume
79 coord_transform = np.dot(tri_mat, trans)
81 return tri_mat, coord_transform
82 else:
83 return cell, None
86class LAMMPSlib(Calculator):
87 r"""
88**Introduction**
90LAMMPSlib is an interface and calculator for LAMMPS_. LAMMPSlib uses
91the python interface that comes with LAMMPS to solve an atoms model
92for energy, atom forces and cell stress. This calculator creates a
93'.lmp' object which is a running lammps program, so further commands
94can be sent to this object executed until it is explicitly closed. Any
95additional variables calculated by lammps can also be extracted. This
96is still experimental code.
98**Arguments**
100======================= ======================================================
101Keyword Description
102======================= ======================================================
103``lmpcmds`` list of strings of LAMMPS commands. You need to supply
104 enough to define the potential to be used e.g.
106 ["pair_style eam/alloy",
107 "pair_coeff * * potentials/NiAlH_jea.eam.alloy Ni Al"]
109``atom_types`` dictionary of ``atomic_symbol :lammps_atom_type``
110 pairs, e.g. ``{'Cu':1}`` to bind copper to lammps
111 atom type 1. If <None>, autocreated by assigning
112 lammps atom types in order that they appear in the
113 first used atoms object.
115``atom_type_masses`` dictionary of ``atomic_symbol :mass`` pairs, e.g.
116 ``{'Cu':63.546}`` to optionally assign masses that
117 override default ase.data.atomic_masses. Note that
118 since unit conversion is done automatically in this
119 module, these quantities must be given in the
120 standard ase mass units (g/mol)
122``log_file`` string
123 path to the desired LAMMPS log file
125``lammps_header`` string to use for lammps setup. Default is to use
126 metal units and simple atom simulation.
128 lammps_header=['units metal',
129 'atom_style atomic',
130 'atom_modify map array sort 0 0'])
132``amendments`` extra list of strings of LAMMPS commands to be run
133 post initialization. (Use: Initialization amendments)
134 e.g.
136 ["mass 1 58.6934"]
138``post_changebox_cmds`` extra list of strings of LAMMPS commands to be run
139 after any LAMMPS 'change_box' command is performed by
140 the calculator. This is relevant because some
141 potentials either themselves depend on the geometry
142 and boundary conditions of the simulation box, or are
143 frequently coupled with other LAMMPS commands that
144 do, e.g. the 'buck/coul/long' pair style is often
145 used with the kspace_* commands, which are sensitive
146 to the periodicity of the simulation box.
148``keep_alive`` Boolean
149 whether to keep the lammps routine alive for more
150 commands. Default is True.
152======================= ======================================================
155**Requirements**
157To run this calculator you must have LAMMPS installed and compiled to
158enable the python interface. See the LAMMPS manual.
160If the following code runs then lammps is installed correctly.
162 >>> from lammps import lammps
163 >>> lmp = lammps()
165The version of LAMMPS is also important. LAMMPSlib is suitable for
166versions after approximately 2011. Prior to this the python interface
167is slightly different from that used by LAMMPSlib. It is not difficult
168to change to the earlier format.
170**LAMMPS and LAMMPSlib**
172The LAMMPS calculator is another calculator that uses LAMMPS (the
173program) to calculate the energy by generating input files and running
174a separate LAMMPS job to perform the analysis. The output data is then
175read back into python. LAMMPSlib makes direct use of the LAMMPS (the
176program) python interface. As well as directly running any LAMMPS
177command line it allows the values of any of LAMMPS variables to be
178extracted and returned to python.
180**Example**
182Provided that the respective potential file is in the working directory, one
183can simply run (note that LAMMPS needs to be compiled to work with EAM
184potentials)
186::
188 from ase import Atom, Atoms
189 from ase.build import bulk
190 from ase.calculators.lammpslib import LAMMPSlib
192 cmds = ["pair_style eam/alloy",
193 "pair_coeff * * NiAlH_jea.eam.alloy Ni H"]
195 Ni = bulk('Ni', cubic=True)
196 H = Atom('H', position=Ni.cell.diagonal()/2)
197 NiH = Ni + H
199 lammps = LAMMPSlib(lmpcmds=cmds, log_file='test.log')
201 NiH.calc = lammps
202 print("Energy ", NiH.get_potential_energy())
205**Implementation**
207LAMMPS provides a set of python functions to allow execution of the
208underlying C++ LAMMPS code. The functions used by the LAMMPSlib
209interface are::
211 from lammps import lammps
213 lmp = lammps(cmd_args) # initiate LAMMPS object with command line args
215 lmp.scatter_atoms('x',1,3,positions) # atom coords to LAMMPS C array
216 lmp.command(cmd) # executes a one line cmd string
217 lmp.extract_variable(...) # extracts a per atom variable
218 lmp.extract_global(...) # extracts a global variable
219 lmp.close() # close the lammps object
221For a single Ni atom model the following lammps file commands would be run
222by invoking the get_potential_energy() method::
224 units metal
225 atom_style atomic
226 atom_modify map array sort 0 0
228 region cell prism 0 xhi 0 yhi 0 zhi xy xz yz units box
229 create_box 1 cell
230 create_atoms 1 single 0 0 0 units box
231 mass * 1.0
233 ## user lmpcmds get executed here
234 pair_style eam/alloy
235 pair_coeff * * NiAlH_jea.eam.alloy Ni
236 ## end of user lmmpcmds
238 run 0
240where xhi, yhi and zhi are the lattice vector lengths and xy,
241xz and yz are the tilt of the lattice vectors, all to be edited.
244**Notes**
246.. _LAMMPS: http://lammps.sandia.gov/
248* Units: The default lammps_header sets the units to Angstrom and eV
249 and for compatibility with ASE Stress is in GPa.
251* The global energy is currently extracted from LAMMPS using
252 extract_variable since lammps.lammps currently extract_global only
253 accepts the following ['dt', 'boxxlo', 'boxxhi', 'boxylo', 'boxyhi',
254 'boxzlo', 'boxzhi', 'natoms', 'nlocal'].
256* If an error occurs while lammps is in control it will crash
257 Python. Check the output of the log file to find the lammps error.
259* If the are commands directly sent to the LAMMPS object this may
260 change the energy value of the model. However the calculator will not
261 know of it and still return the original energy value.
263"""
265 implemented_properties = ['energy', 'free_energy', 'forces', 'stress',
266 'energies']
268 started = False
269 initialized = False
271 default_parameters = dict(
272 atom_types=None,
273 atom_type_masses=None,
274 log_file=None,
275 lammps_name='',
276 keep_alive=True,
277 lammps_header=['units metal',
278 'atom_style atomic',
279 'atom_modify map array sort 0 0'],
280 amendments=None,
281 post_changebox_cmds=None,
282 boundary=True,
283 create_box=True,
284 create_atoms=True,
285 read_molecular_info=False,
286 comm=None)
288 def __init__(self, *args, **kwargs):
289 Calculator.__init__(self, *args, **kwargs)
290 self.lmp = None
292 def __enter__(self):
293 return self
295 def __exit__(self, *args):
296 self.clean()
298 def clean(self):
299 if self.started:
300 self.lmp.close()
301 self.started = False
302 self.initialized = False
303 self.lmp = None
305 def set_cell(self, atoms: Atoms, change: bool = False):
306 self.prism = Prism(atoms.cell, atoms.pbc)
307 _ = self.prism.get_lammps_prism()
308 xhi, yhi, zhi, xy, xz, yz = convert(_, "distance", "ASE", self.units)
309 box_hi = [xhi, yhi, zhi]
311 if change:
312 cell_cmd = ('change_box all '
313 'x final 0 {} y final 0 {} z final 0 {} '
314 'xy final {} xz final {} yz final {} units box'
315 ''.format(xhi, yhi, zhi, xy, xz, yz))
316 if self.parameters.post_changebox_cmds is not None:
317 for cmd in self.parameters.post_changebox_cmds:
318 self.lmp.command(cmd)
319 else:
320 # just in case we'll want to run with a funny shape box,
321 # and here command will only happen once, and before
322 # any calculation
323 if self.parameters.create_box:
324 self.lmp.command('box tilt large')
326 # Check if there are any indefinite boundaries. If so,
327 # shrink-wrapping will end up being used, but we want to
328 # define the LAMMPS region and box fairly tight around the
329 # atoms to avoid losing any
330 lammps_boundary_conditions = self.lammpsbc(atoms).split()
331 if 's' in lammps_boundary_conditions:
332 pos = self.prism.vector_to_lammps(atoms.positions)
333 posmin = np.amin(pos, axis=0)
334 posmax = np.amax(pos, axis=0)
336 for i in range(3):
337 if lammps_boundary_conditions[i] == 's':
338 box_hi[i] = 1.05 * abs(posmax[i] - posmin[i])
340 cell_cmd = ('region cell prism '
341 '0 {} 0 {} 0 {} '
342 '{} {} {} units box'
343 ''.format(*box_hi, xy, xz, yz))
345 self.lmp.command(cell_cmd)
347 def set_lammps_pos(self, atoms: Atoms):
348 # Create local copy of positions that are wrapped along any periodic
349 # directions
350 pos = convert(atoms.positions, "distance", "ASE", self.units)
352 # wrap only after scaling and rotating to reduce chances of
353 # lammps neighbor list bugs.
354 pos = self.prism.vector_to_lammps(pos, wrap=True)
356 # Convert ase position matrix to lammps-style position array
357 # contiguous in memory
358 lmp_positions = list(pos.ravel())
360 # Convert that lammps-style array into a C object
361 c_double_array = (ctypes.c_double * len(lmp_positions))
362 lmp_c_positions = c_double_array(*lmp_positions)
363 # self.lmp.put_coosrds(lmp_c_positions)
364 self.lmp.scatter_atoms('x', 1, 3, lmp_c_positions)
366 def calculate(self, atoms, properties, system_changes):
367 self.propagate(atoms, properties, system_changes, 0)
369 def propagate(self, atoms, properties, system_changes, n_steps, dt=None,
370 dt_not_real_time=False, velocity_field=None):
371 """"atoms: Atoms object
372 Contains positions, unit-cell, ...
373 properties: list of str
374 List of what needs to be calculated. Can be any combination
375 of 'energy', 'forces', 'stress', 'dipole', 'charges', 'magmom'
376 and 'magmoms'.
377 system_changes: list of str
378 List of what has changed since last calculation. Can be
379 any combination of these five: 'positions', 'numbers', 'cell',
380 'pbc', 'initial_charges' and 'initial_magmoms'.
381 """
382 if len(system_changes) == 0:
383 return
385 if not self.started:
386 self.start_lammps()
387 if not self.initialized:
388 self.initialise_lammps(atoms)
389 else: # still need to reset cell
390 # NOTE: The whole point of ``post_changebox_cmds`` is that they're
391 # executed after any call to LAMMPS' change_box command. Here, we
392 # rely on the fact that self.set_cell(), where we have currently
393 # placed the execution of ``post_changebox_cmds``, gets called
394 # after this initial change_box call.
396 # Apply only requested boundary condition changes. Note this needs
397 # to happen before the call to set_cell since 'change_box' will
398 # apply any shrink-wrapping *after* it's updated the cell
399 # dimensions
400 if 'pbc' in system_changes:
401 change_box_str = 'change_box all boundary {}'
402 change_box_cmd = change_box_str.format(self.lammpsbc(atoms))
403 self.lmp.command(change_box_cmd)
405 # Reset positions so that if they are crazy from last
406 # propagation, change_box (in set_cell()) won't hang.
407 # Could do this only after testing for crazy positions?
408 # Could also use scatter_atoms() to set values (requires
409 # MPI comm), or extra_atoms() to get pointers to local
410 # data structures to zero, but then we would have to be
411 # careful with parallelism.
412 self.lmp.command("set atom * x 0.0 y 0.0 z 0.0")
413 self.set_cell(atoms, change=True)
415 if self.parameters.atom_types is None:
416 raise NameError("atom_types are mandatory.")
418 do_rebuild = (
419 not np.array_equal(atoms.numbers, self.previous_atoms_numbers)
420 or any(_ in system_changes for _ in ('numbers', 'initial_charges'))
421 )
422 if not do_rebuild:
423 do_redo_atom_types = not np.array_equal(
424 atoms.numbers, self.previous_atoms_numbers)
425 else:
426 do_redo_atom_types = False
428 self.lmp.command('echo none') # don't echo the atom positions
429 if do_rebuild:
430 self.rebuild(atoms)
431 elif do_redo_atom_types:
432 self.redo_atom_types(atoms)
433 self.lmp.command('echo log') # switch back log
435 self.set_lammps_pos(atoms)
437 if self.parameters.amendments is not None:
438 for cmd in self.parameters.amendments:
439 self.lmp.command(cmd)
441 if n_steps > 0:
442 if velocity_field is None:
443 vel = convert(
444 atoms.get_velocities(),
445 "velocity",
446 "ASE",
447 self.units)
448 else:
449 # FIXME: Do we need to worry about converting to lammps units
450 # here?
451 vel = atoms.arrays[velocity_field]
453 # Transform the velocities to new coordinate system
454 vel = self.prism.vector_to_lammps(vel)
456 # Convert ase velocities matrix to lammps-style velocities array
457 lmp_velocities = list(vel.ravel())
459 # Convert that lammps-style array into a C object
460 c_double_array = (ctypes.c_double * len(lmp_velocities))
461 lmp_c_velocities = c_double_array(*lmp_velocities)
462 self.lmp.scatter_atoms('v', 1, 3, lmp_c_velocities)
464 # Run for 0 time to calculate
465 if dt is not None:
466 if dt_not_real_time:
467 self.lmp.command('timestep %.30f' % dt)
468 else:
469 self.lmp.command('timestep %.30f' %
470 convert(dt, "time", "ASE", self.units))
471 self.lmp.command('run %d' % n_steps)
473 if n_steps > 0:
474 # TODO this must be slower than native copy, but why is it broken?
475 pos = np.array(
476 [x for x in self.lmp.gather_atoms("x", 1, 3)]).reshape(-1, 3)
477 pos = self.prism.vector_to_ase(pos)
479 # Convert from LAMMPS units to ASE units
480 pos = convert(pos, "distance", self.units, "ASE")
482 atoms.set_positions(pos)
484 vel = np.array(
485 [v for v in self.lmp.gather_atoms("v", 1, 3)]).reshape(-1, 3)
486 vel = self.prism.vector_to_lammps(vel)
487 if velocity_field is None:
488 atoms.set_velocities(convert(vel, 'velocity', self.units,
489 'ASE'))
491 # Extract the forces and energy
492 self.results['energy'] = convert(
493 self.lmp.extract_variable('pe', None, 0),
494 "energy", self.units, "ASE"
495 )
496 self.results['free_energy'] = self.results['energy']
498 ids = self.lmp.numpy.extract_atom("id")
499 # if ids doesn't match atoms then data is MPI distributed, which
500 # we can't handle
501 assert len(ids) == len(atoms)
502 self.results["energies"] = convert(
503 self.lmp.numpy.extract_compute('pe_peratom',
504 self.LMP_STYLE_ATOM,
505 self.LMP_TYPE_VECTOR),
506 "energy", self.units, "ASE"
507 )
508 self.results["energies"][ids - 1] = self.results["energies"]
510 stress = np.empty(6)
511 stress_vars = ['pxx', 'pyy', 'pzz', 'pyz', 'pxz', 'pxy']
513 for i, var in enumerate(stress_vars):
514 stress[i] = self.lmp.extract_variable(var, None, 0)
516 stress_mat = np.zeros((3, 3))
517 stress_mat[0, 0] = stress[0]
518 stress_mat[1, 1] = stress[1]
519 stress_mat[2, 2] = stress[2]
520 stress_mat[1, 2] = stress[3]
521 stress_mat[2, 1] = stress[3]
522 stress_mat[0, 2] = stress[4]
523 stress_mat[2, 0] = stress[4]
524 stress_mat[0, 1] = stress[5]
525 stress_mat[1, 0] = stress[5]
527 stress_mat = self.prism.tensor2_to_ase(stress_mat)
529 stress[0] = stress_mat[0, 0]
530 stress[1] = stress_mat[1, 1]
531 stress[2] = stress_mat[2, 2]
532 stress[3] = stress_mat[1, 2]
533 stress[4] = stress_mat[0, 2]
534 stress[5] = stress_mat[0, 1]
536 self.results['stress'] = convert(-stress, "pressure", self.units, "ASE")
538 # definitely yields atom-id ordered force array
539 f = convert(np.array(self.lmp.gather_atoms("f", 1, 3)).reshape(-1, 3),
540 "force", self.units, "ASE")
541 self.results['forces'] = self.prism.vector_to_ase(f)
543 # otherwise check_state will always trigger a new calculation
544 self.atoms = atoms.copy()
546 if not self.parameters["keep_alive"]:
547 self.clean()
549 def lammpsbc(self, atoms):
550 """Determine LAMMPS boundary types based on ASE pbc settings. For
551 non-periodic dimensions, if the cell length is finite then
552 fixed BCs ('f') are used; if the cell length is approximately
553 zero, shrink-wrapped BCs ('s') are used."""
555 retval = ''
556 pbc = atoms.get_pbc()
557 if np.all(pbc):
558 retval = 'p p p'
559 else:
560 cell = atoms.get_cell()
561 for i in range(3):
562 if pbc[i]:
563 retval += 'p '
564 else:
565 # See if we're using indefinite ASE boundaries along this
566 # direction
567 if np.linalg.norm(cell[i]) < np.finfo(cell[i][0]).tiny:
568 retval += 's '
569 else:
570 retval += 'f '
572 return retval.strip()
574 def rebuild(self, atoms):
575 try:
576 n_diff = len(atoms.numbers) - len(self.previous_atoms_numbers)
577 except Exception: # XXX Which kind of exception?
578 n_diff = len(atoms.numbers)
580 if n_diff > 0:
581 if any(("reax/c" in cmd) for cmd in self.parameters.lmpcmds):
582 self.lmp.command("pair_style lj/cut 2.5")
583 self.lmp.command("pair_coeff * * 1 1")
585 for cmd in self.parameters.lmpcmds:
586 if (("pair_style" in cmd) or ("pair_coeff" in cmd) or
587 ("qeq/reax" in cmd)):
588 self.lmp.command(cmd)
590 cmd = f"create_atoms 1 random {n_diff} 1 NULL"
591 self.lmp.command(cmd)
592 elif n_diff < 0:
593 cmd = "group delatoms id {}:{}".format(
594 len(atoms.numbers) + 1, len(self.previous_atoms_numbers))
595 self.lmp.command(cmd)
596 cmd = "delete_atoms group delatoms"
597 self.lmp.command(cmd)
599 self.redo_atom_types(atoms)
601 def redo_atom_types(self, atoms):
602 current_types = {
603 (i + 1, self.parameters.atom_types[sym]) for i, sym
604 in enumerate(atoms.get_chemical_symbols())}
606 try:
607 previous_types = {
608 (i + 1, self.parameters.atom_types[ase_chemical_symbols[Z]])
609 for i, Z in enumerate(self.previous_atoms_numbers)}
610 except Exception: # XXX which kind of exception?
611 previous_types = set()
613 for (i, i_type) in current_types - previous_types:
614 cmd = f"set atom {i} type {i_type}"
615 self.lmp.command(cmd)
617 # set charges only when LAMMPS `atom_style` permits charges
618 # https://docs.lammps.org/Library_properties.html#extract-atom-flags
619 if self.lmp.extract_setting('q_flag') == 1:
620 charges = atoms.get_initial_charges()
621 if np.any(charges != 0.0):
622 for i, q in enumerate(charges):
623 self.lmp.command(f'set atom {i + 1} charge {q}')
625 self.previous_atoms_numbers = atoms.numbers.copy()
627 def restart_lammps(self, atoms):
628 if self.started:
629 self.lmp.command("clear")
630 # hope there's no other state to be reset
631 self.started = False
632 self.initialized = False
633 self.previous_atoms_numbers = []
634 self.start_lammps()
635 self.initialise_lammps(atoms)
637 def start_lammps(self):
638 # Only import lammps when running a calculation
639 # so it is not required to use other parts of the
640 # module
641 from lammps import LMP_STYLE_ATOM, LMP_TYPE_VECTOR, lammps
643 self.LMP_STYLE_ATOM = LMP_STYLE_ATOM
644 self.LMP_TYPE_VECTOR = LMP_TYPE_VECTOR
646 # start lammps process
647 if self.parameters.log_file is None:
648 cmd_args = ['-echo', 'log', '-log', 'none', '-screen', 'none',
649 '-nocite']
650 else:
651 cmd_args = ['-echo', 'log', '-log', self.parameters.log_file,
652 '-screen', 'none', '-nocite']
654 self.cmd_args = cmd_args
656 if self.lmp is None:
657 self.lmp = lammps(self.parameters.lammps_name, self.cmd_args,
658 comm=self.parameters.comm)
660 # Run header commands to set up lammps (units, etc.)
661 for cmd in self.parameters.lammps_header:
662 self.lmp.command(cmd)
664 for cmd in self.parameters.lammps_header:
665 if "units" in cmd:
666 self.units = cmd.split()[1]
668 if 'lammps_header_extra' in self.parameters:
669 if self.parameters.lammps_header_extra is not None:
670 for cmd in self.parameters.lammps_header_extra:
671 self.lmp.command(cmd)
673 self.started = True
675 def initialise_lammps(self, atoms):
676 # Initialising commands
677 if self.parameters.boundary:
678 # if the boundary command is in the supplied commands use that
679 # otherwise use atoms pbc
680 for cmd in self.parameters.lmpcmds:
681 if 'boundary' in cmd:
682 break
683 else:
684 self.lmp.command('boundary ' + self.lammpsbc(atoms))
686 # Initialize cell
687 self.set_cell(atoms, change=not self.parameters.create_box)
689 if self.parameters.atom_types is None:
690 # if None is given, create from atoms object in order of appearance
691 s = atoms.get_chemical_symbols()
692 _, idx = np.unique(s, return_index=True)
693 s_red = np.array(s)[np.sort(idx)].tolist()
694 self.parameters.atom_types = {j: i + 1 for i, j in enumerate(s_red)}
696 # Initialize box
697 if self.parameters.create_box:
698 # count number of known types
699 n_types = len(self.parameters.atom_types)
700 create_box_command = f'create_box {n_types} cell'
701 self.lmp.command(create_box_command)
703 # Initialize the atoms with their types
704 # positions do not matter here
705 if self.parameters.create_atoms:
706 self.lmp.command('echo none') # don't echo the atom positions
707 self.rebuild(atoms)
708 self.lmp.command('echo log') # turn back on
709 else:
710 self.previous_atoms_numbers = atoms.numbers.copy()
712 # execute the user commands
713 for cmd in self.parameters.lmpcmds + ["compute pe_peratom all pe/atom"]:
714 self.lmp.command(cmd)
716 # Set masses after user commands, e.g. to override
717 # EAM-provided masses
718 for sym in self.parameters.atom_types:
719 if self.parameters.atom_type_masses is None:
720 mass = ase_atomic_masses[ase_atomic_numbers[sym]]
721 else:
722 mass = self.parameters.atom_type_masses[sym]
723 self.lmp.command('mass %d %.30f' % (
724 self.parameters.atom_types[sym],
725 convert(mass, "mass", "ASE", self.units)))
727 # Define force & energy variables for extraction
728 self.lmp.command('variable pxx equal pxx')
729 self.lmp.command('variable pyy equal pyy')
730 self.lmp.command('variable pzz equal pzz')
731 self.lmp.command('variable pxy equal pxy')
732 self.lmp.command('variable pxz equal pxz')
733 self.lmp.command('variable pyz equal pyz')
735 # I am not sure why we need this next line but LAMMPS will
736 # raise an error if it is not there. Perhaps it is needed to
737 # ensure the cell stresses are calculated
738 self.lmp.command('thermo_style custom pe pxx emol ecoul')
740 self.lmp.command('variable fx atom fx')
741 self.lmp.command('variable fy atom fy')
742 self.lmp.command('variable fz atom fz')
744 # do we need this if we extract from a global ?
745 self.lmp.command('variable pe equal pe')
747 self.lmp.command("neigh_modify delay 0 every 1 check yes")
749 self.initialized = True