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