Coverage for /builds/kinetik161/ase/ase/calculators/lammpslib.py: 73.50%
283 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +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(0, 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', 'charges' and '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 = (not np.array_equal(atoms.numbers,
419 self.previous_atoms_numbers)
420 or ("numbers" in system_changes))
421 if not do_rebuild:
422 do_redo_atom_types = not np.array_equal(
423 atoms.numbers, self.previous_atoms_numbers)
424 else:
425 do_redo_atom_types = False
427 self.lmp.command('echo none') # don't echo the atom positions
428 if do_rebuild:
429 self.rebuild(atoms)
430 elif do_redo_atom_types:
431 self.redo_atom_types(atoms)
432 self.lmp.command('echo log') # switch back log
434 self.set_lammps_pos(atoms)
436 if self.parameters.amendments is not None:
437 for cmd in self.parameters.amendments:
438 self.lmp.command(cmd)
440 if n_steps > 0:
441 if velocity_field is None:
442 vel = convert(
443 atoms.get_velocities(),
444 "velocity",
445 "ASE",
446 self.units)
447 else:
448 # FIXME: Do we need to worry about converting to lammps units
449 # here?
450 vel = atoms.arrays[velocity_field]
452 # Transform the velocities to new coordinate system
453 vel = self.prism.vector_to_lammps(vel)
455 # Convert ase velocities matrix to lammps-style velocities array
456 lmp_velocities = list(vel.ravel())
458 # Convert that lammps-style array into a C object
459 c_double_array = (ctypes.c_double * len(lmp_velocities))
460 lmp_c_velocities = c_double_array(*lmp_velocities)
461 self.lmp.scatter_atoms('v', 1, 3, lmp_c_velocities)
463 # Run for 0 time to calculate
464 if dt is not None:
465 if dt_not_real_time:
466 self.lmp.command('timestep %.30f' % dt)
467 else:
468 self.lmp.command('timestep %.30f' %
469 convert(dt, "time", "ASE", self.units))
470 self.lmp.command('run %d' % n_steps)
472 if n_steps > 0:
473 # TODO this must be slower than native copy, but why is it broken?
474 pos = np.array(
475 [x for x in self.lmp.gather_atoms("x", 1, 3)]).reshape(-1, 3)
476 pos = self.prism.vector_to_ase(pos)
478 # Convert from LAMMPS units to ASE units
479 pos = convert(pos, "distance", self.units, "ASE")
481 atoms.set_positions(pos)
483 vel = np.array(
484 [v for v in self.lmp.gather_atoms("v", 1, 3)]).reshape(-1, 3)
485 vel = self.prism.vector_to_lammps(vel)
486 if velocity_field is None:
487 atoms.set_velocities(convert(vel, 'velocity', self.units,
488 'ASE'))
490 # Extract the forces and energy
491 self.results['energy'] = convert(
492 self.lmp.extract_variable('pe', None, 0),
493 "energy", self.units, "ASE"
494 )
495 self.results['free_energy'] = self.results['energy']
497 ids = self.lmp.numpy.extract_atom("id")
498 # if ids doesn't match atoms then data is MPI distributed, which
499 # we can't handle
500 assert len(ids) == len(atoms)
501 self.results["energies"] = convert(
502 self.lmp.numpy.extract_compute('pe_peratom',
503 self.LMP_STYLE_ATOM,
504 self.LMP_TYPE_VECTOR),
505 "energy", self.units, "ASE"
506 )
507 self.results["energies"][ids - 1] = self.results["energies"]
509 stress = np.empty(6)
510 stress_vars = ['pxx', 'pyy', 'pzz', 'pyz', 'pxz', 'pxy']
512 for i, var in enumerate(stress_vars):
513 stress[i] = self.lmp.extract_variable(var, None, 0)
515 stress_mat = np.zeros((3, 3))
516 stress_mat[0, 0] = stress[0]
517 stress_mat[1, 1] = stress[1]
518 stress_mat[2, 2] = stress[2]
519 stress_mat[1, 2] = stress[3]
520 stress_mat[2, 1] = stress[3]
521 stress_mat[0, 2] = stress[4]
522 stress_mat[2, 0] = stress[4]
523 stress_mat[0, 1] = stress[5]
524 stress_mat[1, 0] = stress[5]
526 stress_mat = self.prism.tensor2_to_ase(stress_mat)
528 stress[0] = stress_mat[0, 0]
529 stress[1] = stress_mat[1, 1]
530 stress[2] = stress_mat[2, 2]
531 stress[3] = stress_mat[1, 2]
532 stress[4] = stress_mat[0, 2]
533 stress[5] = stress_mat[0, 1]
535 self.results['stress'] = convert(-stress, "pressure", self.units, "ASE")
537 # definitely yields atom-id ordered force array
538 f = convert(np.array(self.lmp.gather_atoms("f", 1, 3)).reshape(-1, 3),
539 "force", self.units, "ASE")
540 self.results['forces'] = self.prism.vector_to_ase(f)
542 # otherwise check_state will always trigger a new calculation
543 self.atoms = atoms.copy()
545 if not self.parameters["keep_alive"]:
546 self.clean()
548 def lammpsbc(self, atoms):
549 """Determine LAMMPS boundary types based on ASE pbc settings. For
550 non-periodic dimensions, if the cell length is finite then
551 fixed BCs ('f') are used; if the cell length is approximately
552 zero, shrink-wrapped BCs ('s') are used."""
554 retval = ''
555 pbc = atoms.get_pbc()
556 if np.all(pbc):
557 retval = 'p p p'
558 else:
559 cell = atoms.get_cell()
560 for i in range(0, 3):
561 if pbc[i]:
562 retval += 'p '
563 else:
564 # See if we're using indefinite ASE boundaries along this
565 # direction
566 if np.linalg.norm(cell[i]) < np.finfo(cell[i][0]).tiny:
567 retval += 's '
568 else:
569 retval += 'f '
571 return retval.strip()
573 def rebuild(self, atoms):
574 try:
575 n_diff = len(atoms.numbers) - len(self.previous_atoms_numbers)
576 except Exception: # XXX Which kind of exception?
577 n_diff = len(atoms.numbers)
579 if n_diff > 0:
580 if any(("reax/c" in cmd) for cmd in self.parameters.lmpcmds):
581 self.lmp.command("pair_style lj/cut 2.5")
582 self.lmp.command("pair_coeff * * 1 1")
584 for cmd in self.parameters.lmpcmds:
585 if (("pair_style" in cmd) or ("pair_coeff" in cmd) or
586 ("qeq/reax" in cmd)):
587 self.lmp.command(cmd)
589 cmd = f"create_atoms 1 random {n_diff} 1 NULL"
590 self.lmp.command(cmd)
591 elif n_diff < 0:
592 cmd = "group delatoms id {}:{}".format(
593 len(atoms.numbers) + 1, len(self.previous_atoms_numbers))
594 self.lmp.command(cmd)
595 cmd = "delete_atoms group delatoms"
596 self.lmp.command(cmd)
598 self.redo_atom_types(atoms)
600 def redo_atom_types(self, atoms):
601 current_types = {
602 (i + 1, self.parameters.atom_types[sym]) for i, sym
603 in enumerate(atoms.get_chemical_symbols())}
605 try:
606 previous_types = {
607 (i + 1, self.parameters.atom_types[ase_chemical_symbols[Z]])
608 for i, Z in enumerate(self.previous_atoms_numbers)}
609 except Exception: # XXX which kind of exception?
610 previous_types = set()
612 for (i, i_type) in current_types - previous_types:
613 cmd = f"set atom {i} type {i_type}"
614 self.lmp.command(cmd)
616 self.previous_atoms_numbers = atoms.numbers.copy()
618 def restart_lammps(self, atoms):
619 if self.started:
620 self.lmp.command("clear")
621 # hope there's no other state to be reset
622 self.started = False
623 self.initialized = False
624 self.previous_atoms_numbers = []
625 self.start_lammps()
626 self.initialise_lammps(atoms)
628 def start_lammps(self):
629 # Only import lammps when running a calculation
630 # so it is not required to use other parts of the
631 # module
632 from lammps import LMP_STYLE_ATOM, LMP_TYPE_VECTOR, lammps
634 self.LMP_STYLE_ATOM = LMP_STYLE_ATOM
635 self.LMP_TYPE_VECTOR = LMP_TYPE_VECTOR
637 # start lammps process
638 if self.parameters.log_file is None:
639 cmd_args = ['-echo', 'log', '-log', 'none', '-screen', 'none',
640 '-nocite']
641 else:
642 cmd_args = ['-echo', 'log', '-log', self.parameters.log_file,
643 '-screen', 'none', '-nocite']
645 self.cmd_args = cmd_args
647 if self.lmp is None:
648 self.lmp = lammps(self.parameters.lammps_name, self.cmd_args,
649 comm=self.parameters.comm)
651 # Run header commands to set up lammps (units, etc.)
652 for cmd in self.parameters.lammps_header:
653 self.lmp.command(cmd)
655 for cmd in self.parameters.lammps_header:
656 if "units" in cmd:
657 self.units = cmd.split()[1]
659 if 'lammps_header_extra' in self.parameters:
660 if self.parameters.lammps_header_extra is not None:
661 for cmd in self.parameters.lammps_header_extra:
662 self.lmp.command(cmd)
664 self.started = True
666 def initialise_lammps(self, atoms):
667 # Initialising commands
668 if self.parameters.boundary:
669 # if the boundary command is in the supplied commands use that
670 # otherwise use atoms pbc
671 for cmd in self.parameters.lmpcmds:
672 if 'boundary' in cmd:
673 break
674 else:
675 self.lmp.command('boundary ' + self.lammpsbc(atoms))
677 # Initialize cell
678 self.set_cell(atoms, change=not self.parameters.create_box)
680 if self.parameters.atom_types is None:
681 # if None is given, create from atoms object in order of appearance
682 s = atoms.get_chemical_symbols()
683 _, idx = np.unique(s, return_index=True)
684 s_red = np.array(s)[np.sort(idx)].tolist()
685 self.parameters.atom_types = {j: i + 1 for i, j in enumerate(s_red)}
687 # Initialize box
688 if self.parameters.create_box:
689 # count number of known types
690 n_types = len(self.parameters.atom_types)
691 create_box_command = f'create_box {n_types} cell'
692 self.lmp.command(create_box_command)
694 # Initialize the atoms with their types
695 # positions do not matter here
696 if self.parameters.create_atoms:
697 self.lmp.command('echo none') # don't echo the atom positions
698 self.rebuild(atoms)
699 self.lmp.command('echo log') # turn back on
700 else:
701 self.previous_atoms_numbers = atoms.numbers.copy()
703 # execute the user commands
704 for cmd in self.parameters.lmpcmds + ["compute pe_peratom all pe/atom"]:
705 self.lmp.command(cmd)
707 # Set masses after user commands, e.g. to override
708 # EAM-provided masses
709 for sym in self.parameters.atom_types:
710 if self.parameters.atom_type_masses is None:
711 mass = ase_atomic_masses[ase_atomic_numbers[sym]]
712 else:
713 mass = self.parameters.atom_type_masses[sym]
714 self.lmp.command('mass %d %.30f' % (
715 self.parameters.atom_types[sym],
716 convert(mass, "mass", "ASE", self.units)))
718 # Define force & energy variables for extraction
719 self.lmp.command('variable pxx equal pxx')
720 self.lmp.command('variable pyy equal pyy')
721 self.lmp.command('variable pzz equal pzz')
722 self.lmp.command('variable pxy equal pxy')
723 self.lmp.command('variable pxz equal pxz')
724 self.lmp.command('variable pyz equal pyz')
726 # I am not sure why we need this next line but LAMMPS will
727 # raise an error if it is not there. Perhaps it is needed to
728 # ensure the cell stresses are calculated
729 self.lmp.command('thermo_style custom pe pxx emol ecoul')
731 self.lmp.command('variable fx atom fx')
732 self.lmp.command('variable fy atom fy')
733 self.lmp.command('variable fz atom fz')
735 # do we need this if we extract from a global ?
736 self.lmp.command('variable pe equal pe')
738 self.lmp.command("neigh_modify delay 0 every 1 check yes")
740 self.initialized = True