Coverage for /builds/debichem-team/python-ase/ase/calculators/lammpslib.py : 77.50%
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', 'forces', 'stress']
254 started = False
255 initialized = False
257 default_parameters = dict(
258 atom_types=None,
259 atom_type_masses=None,
260 log_file=None,
261 lammps_name='',
262 keep_alive=False,
263 lammps_header=['units metal',
264 'atom_style atomic',
265 'atom_modify map array sort 0 0'],
266 amendments=None,
267 post_changebox_cmds=None,
268 boundary=True,
269 create_box=True,
270 create_atoms=True,
271 read_molecular_info=False,
272 comm=None)
274 def __init__(self, *args, **kwargs):
275 Calculator.__init__(self, *args, **kwargs)
276 self.lmp = None
278 def __del__(self):
279 if self.started:
280 self.lmp.close()
281 self.started = False
282 self.lmp = None
284 def set_cell(self, atoms, change=False):
285 lammps_cell, self.coord_transform = convert_cell(atoms.get_cell())
287 xhi, xy, xz, _, yhi, yz, _, _, zhi = convert(
288 lammps_cell.flatten(order='C'), "distance", "ASE", self.units)
289 box_hi = [xhi, yhi, zhi]
291 if change:
292 cell_cmd = ('change_box all '
293 'x final 0 {} y final 0 {} z final 0 {} '
294 'xy final {} xz final {} yz final {} units box'
295 ''.format(xhi, yhi, zhi, xy, xz, yz))
296 if self.parameters.post_changebox_cmds is not None:
297 for cmd in self.parameters.post_changebox_cmds:
298 self.lmp.command(cmd)
299 else:
300 # just in case we'll want to run with a funny shape box,
301 # and here command will only happen once, and before
302 # any calculation
303 if self.parameters.create_box:
304 self.lmp.command('box tilt large')
306 # Check if there are any indefinite boundaries. If so,
307 # shrink-wrapping will end up being used, but we want to
308 # define the LAMMPS region and box fairly tight around the
309 # atoms to avoid losing any
310 lammps_boundary_conditions = self.lammpsbc(atoms).split()
311 if 's' in lammps_boundary_conditions:
312 pos = atoms.get_positions()
313 if self.coord_transform is not None:
314 pos = np.dot(self.coord_transform, pos.transpose())
315 pos = pos.transpose()
316 posmin = np.amin(pos, axis=0)
317 posmax = np.amax(pos, axis=0)
319 for i in range(0, 3):
320 if lammps_boundary_conditions[i] == 's':
321 box_hi[i] = 1.05 * abs(posmax[i] - posmin[i])
323 cell_cmd = ('region cell prism '
324 '0 {} 0 {} 0 {} '
325 '{} {} {} units box'
326 ''.format(*box_hi, xy, xz, yz))
328 self.lmp.command(cell_cmd)
330 def set_lammps_pos(self, atoms):
331 # Create local copy of positions that are wrapped along any periodic
332 # directions
333 cell = convert(atoms.cell, "distance", "ASE", self.units)
334 pos = convert(atoms.positions, "distance", "ASE", self.units)
336 # If necessary, transform the positions to new coordinate system
337 if self.coord_transform is not None:
338 pos = np.dot(pos, self.coord_transform.T)
339 cell = np.dot(cell, self.coord_transform.T)
341 # wrap only after scaling and rotating to reduce chances of
342 # lammps neighbor list bugs.
343 pos = wrap_positions(pos, cell, atoms.get_pbc())
345 # Convert ase position matrix to lammps-style position array
346 # contiguous in memory
347 lmp_positions = list(pos.ravel())
349 # Convert that lammps-style array into a C object
350 c_double_array = (ctypes.c_double * len(lmp_positions))
351 lmp_c_positions = c_double_array(*lmp_positions)
352 # self.lmp.put_coosrds(lmp_c_positions)
353 self.lmp.scatter_atoms('x', 1, 3, lmp_c_positions)
355 def calculate(self, atoms, properties, system_changes):
356 self.propagate(atoms, properties, system_changes, 0)
358 def propagate(self, atoms, properties, system_changes, n_steps, dt=None,
359 dt_not_real_time=False, velocity_field=None):
360 """"atoms: Atoms object
361 Contains positions, unit-cell, ...
362 properties: list of str
363 List of what needs to be calculated. Can be any combination
364 of 'energy', 'forces', 'stress', 'dipole', 'charges', 'magmom'
365 and 'magmoms'.
366 system_changes: list of str
367 List of what has changed since last calculation. Can be
368 any combination of these five: 'positions', 'numbers', 'cell',
369 'pbc', 'charges' and 'magmoms'.
370 """
371 if len(system_changes) == 0:
372 return
374 self.coord_transform = None
376 if not self.started:
377 self.start_lammps()
378 if not self.initialized:
379 self.initialise_lammps(atoms)
380 else: # still need to reset cell
381 # NOTE: The whole point of ``post_changebox_cmds`` is that they're
382 # executed after any call to LAMMPS' change_box command. Here, we
383 # rely on the fact that self.set_cell(), where we have currently
384 # placed the execution of ``post_changebox_cmds``, gets called
385 # after this initial change_box call.
387 # Apply only requested boundary condition changes. Note this needs
388 # to happen before the call to set_cell since 'change_box' will
389 # apply any shrink-wrapping *after* it's updated the cell
390 # dimensions
391 if 'pbc' in system_changes:
392 change_box_str = 'change_box all boundary {}'
393 change_box_cmd = change_box_str.format(self.lammpsbc(atoms))
394 self.lmp.command(change_box_cmd)
396 # Reset positions so that if they are crazy from last
397 # propagation, change_box (in set_cell()) won't hang.
398 # Could do this only after testing for crazy positions?
399 # Could also use scatter_atoms() to set values (requires
400 # MPI comm), or extra_atoms() to get pointers to local
401 # data structures to zero, but then we would have to be
402 # careful with parallelism.
403 self.lmp.command("set atom * x 0.0 y 0.0 z 0.0")
404 self.set_cell(atoms, change=True)
406 if self.parameters.atom_types is None:
407 raise NameError("atom_types are mandatory.")
409 do_rebuild = (not np.array_equal(atoms.numbers,
410 self.previous_atoms_numbers)
411 or ("numbers" in system_changes))
412 if not do_rebuild:
413 do_redo_atom_types = not np.array_equal(
414 atoms.numbers, self.previous_atoms_numbers)
415 else:
416 do_redo_atom_types = False
418 self.lmp.command('echo none') # don't echo the atom positions
419 if do_rebuild:
420 self.rebuild(atoms)
421 elif do_redo_atom_types:
422 self.redo_atom_types(atoms)
423 self.lmp.command('echo log') # switch back log
425 self.set_lammps_pos(atoms)
427 if self.parameters.amendments is not None:
428 for cmd in self.parameters.amendments:
429 self.lmp.command(cmd)
431 if n_steps > 0:
432 if velocity_field is None:
433 vel = convert(
434 atoms.get_velocities(),
435 "velocity",
436 "ASE",
437 self.units)
438 else:
439 # FIXME: Do we need to worry about converting to lammps units
440 # here?
441 vel = atoms.arrays[velocity_field]
443 # If necessary, transform the velocities to new coordinate system
444 if self.coord_transform is not None:
445 vel = np.dot(self.coord_transform, vel.T).T
447 # Convert ase velocities matrix to lammps-style velocities array
448 lmp_velocities = list(vel.ravel())
450 # Convert that lammps-style array into a C object
451 c_double_array = (ctypes.c_double * len(lmp_velocities))
452 lmp_c_velocities = c_double_array(*lmp_velocities)
453 self.lmp.scatter_atoms('v', 1, 3, lmp_c_velocities)
455 # Run for 0 time to calculate
456 if dt is not None:
457 if dt_not_real_time:
458 self.lmp.command('timestep %.30f' % dt)
459 else:
460 self.lmp.command('timestep %.30f' %
461 convert(dt, "time", "ASE", self.units))
462 self.lmp.command('run %d' % n_steps)
464 if n_steps > 0:
465 # TODO this must be slower than native copy, but why is it broken?
466 pos = np.array(
467 [x for x in self.lmp.gather_atoms("x", 1, 3)]).reshape(-1, 3)
468 if self.coord_transform is not None:
469 pos = np.dot(pos, self.coord_transform)
471 # Convert from LAMMPS units to ASE units
472 pos = convert(pos, "distance", self.units, "ASE")
474 atoms.set_positions(pos)
476 vel = np.array(
477 [v for v in self.lmp.gather_atoms("v", 1, 3)]).reshape(-1, 3)
478 if self.coord_transform is not None:
479 vel = np.dot(vel, self.coord_transform)
480 if velocity_field is None:
481 atoms.set_velocities(convert(vel, 'velocity', self.units,
482 'ASE'))
484 # Extract the forces and energy
485 self.results['energy'] = convert(
486 self.lmp.extract_variable('pe', None, 0),
487 "energy", self.units, "ASE"
488 )
489 self.results['free_energy'] = self.results['energy']
491 stress = np.empty(6)
492 stress_vars = ['pxx', 'pyy', 'pzz', 'pyz', 'pxz', 'pxy']
494 for i, var in enumerate(stress_vars):
495 stress[i] = self.lmp.extract_variable(var, None, 0)
497 stress_mat = np.zeros((3, 3))
498 stress_mat[0, 0] = stress[0]
499 stress_mat[1, 1] = stress[1]
500 stress_mat[2, 2] = stress[2]
501 stress_mat[1, 2] = stress[3]
502 stress_mat[2, 1] = stress[3]
503 stress_mat[0, 2] = stress[4]
504 stress_mat[2, 0] = stress[4]
505 stress_mat[0, 1] = stress[5]
506 stress_mat[1, 0] = stress[5]
507 if self.coord_transform is not None:
508 stress_mat = np.dot(self.coord_transform.T,
509 np.dot(stress_mat, self.coord_transform))
510 stress[0] = stress_mat[0, 0]
511 stress[1] = stress_mat[1, 1]
512 stress[2] = stress_mat[2, 2]
513 stress[3] = stress_mat[1, 2]
514 stress[4] = stress_mat[0, 2]
515 stress[5] = stress_mat[0, 1]
517 self.results['stress'] = convert(-stress, "pressure", self.units, "ASE")
519 # definitely yields atom-id ordered force array
520 f = convert(np.array(self.lmp.gather_atoms("f", 1, 3)).reshape(-1, 3),
521 "force", self.units, "ASE")
523 if self.coord_transform is not None:
524 self.results['forces'] = np.dot(f, self.coord_transform)
525 else:
526 self.results['forces'] = f.copy()
528 # otherwise check_state will always trigger a new calculation
529 self.atoms = atoms.copy()
531 if not self.parameters.keep_alive:
532 self.lmp.close()
534 def lammpsbc(self, atoms):
535 """Determine LAMMPS boundary types based on ASE pbc settings. For
536 non-periodic dimensions, if the cell length is finite then
537 fixed BCs ('f') are used; if the cell length is approximately
538 zero, shrink-wrapped BCs ('s') are used."""
540 retval = ''
541 pbc = atoms.get_pbc()
542 if np.all(pbc):
543 retval = 'p p p'
544 else:
545 cell = atoms.get_cell()
546 for i in range(0, 3):
547 if pbc[i]:
548 retval += 'p '
549 else:
550 # See if we're using indefinite ASE boundaries along this
551 # direction
552 if np.linalg.norm(cell[i]) < np.finfo(cell[i][0]).tiny:
553 retval += 's '
554 else:
555 retval += 'f '
557 return retval.strip()
559 def rebuild(self, atoms):
560 try:
561 n_diff = len(atoms.numbers) - len(self.previous_atoms_numbers)
562 except Exception: # XXX Which kind of exception?
563 n_diff = len(atoms.numbers)
565 if n_diff > 0:
566 if any([("reax/c" in cmd) for cmd in self.parameters.lmpcmds]):
567 self.lmp.command("pair_style lj/cut 2.5")
568 self.lmp.command("pair_coeff * * 1 1")
570 for cmd in self.parameters.lmpcmds:
571 if (("pair_style" in cmd) or ("pair_coeff" in cmd) or
572 ("qeq/reax" in cmd)):
573 self.lmp.command(cmd)
575 cmd = "create_atoms 1 random {} 1 NULL".format(n_diff)
576 self.lmp.command(cmd)
577 elif n_diff < 0:
578 cmd = "group delatoms id {}:{}".format(
579 len(atoms.numbers) + 1, len(self.previous_atoms_numbers))
580 self.lmp.command(cmd)
581 cmd = "delete_atoms group delatoms"
582 self.lmp.command(cmd)
584 self.redo_atom_types(atoms)
586 def redo_atom_types(self, atoms):
587 current_types = set(
588 (i + 1, self.parameters.atom_types[sym]) for i, sym
589 in enumerate(atoms.get_chemical_symbols()))
591 try:
592 previous_types = set(
593 (i + 1, self.parameters.atom_types[ase_chemical_symbols[Z]])
594 for i, Z in enumerate(self.previous_atoms_numbers))
595 except Exception: # XXX which kind of exception?
596 previous_types = set()
598 for (i, i_type) in current_types - previous_types:
599 cmd = "set atom {} type {}".format(i, i_type)
600 self.lmp.command(cmd)
602 self.previous_atoms_numbers = atoms.numbers.copy()
604 def restart_lammps(self, atoms):
605 if self.started:
606 self.lmp.command("clear")
607 # hope there's no other state to be reset
608 self.started = False
609 self.initialized = False
610 self.previous_atoms_numbers = []
611 self.start_lammps()
612 self.initialise_lammps(atoms)
614 def start_lammps(self):
615 # Only import lammps when running a calculation
616 # so it is not required to use other parts of the
617 # module
618 from lammps import lammps
619 # start lammps process
620 if self.parameters.log_file is None:
621 cmd_args = ['-echo', 'log', '-log', 'none', '-screen', 'none',
622 '-nocite']
623 else:
624 cmd_args = ['-echo', 'log', '-log', self.parameters.log_file,
625 '-screen', 'none', '-nocite']
627 self.cmd_args = cmd_args
629 if self.lmp is None:
630 self.lmp = lammps(self.parameters.lammps_name, self.cmd_args,
631 comm=self.parameters.comm)
633 # Run header commands to set up lammps (units, etc.)
634 for cmd in self.parameters.lammps_header:
635 self.lmp.command(cmd)
637 for cmd in self.parameters.lammps_header:
638 if "units" in cmd:
639 self.units = cmd.split()[1]
641 if 'lammps_header_extra' in self.parameters:
642 if self.parameters.lammps_header_extra is not None:
643 for cmd in self.parameters.lammps_header_extra:
644 self.lmp.command(cmd)
646 self.started = True
648 def initialise_lammps(self, atoms):
649 # Initialising commands
650 if self.parameters.boundary:
651 # if the boundary command is in the supplied commands use that
652 # otherwise use atoms pbc
653 for cmd in self.parameters.lmpcmds:
654 if 'boundary' in cmd:
655 break
656 else:
657 self.lmp.command('boundary ' + self.lammpsbc(atoms))
659 # Initialize cell
660 self.set_cell(atoms, change=not self.parameters.create_box)
662 if self.parameters.atom_types is None:
663 # if None is given, create from atoms object in order of appearance
664 s = atoms.get_chemical_symbols()
665 _, idx = np.unique(s, return_index=True)
666 s_red = np.array(s)[np.sort(idx)].tolist()
667 self.parameters.atom_types = {j: i + 1 for i, j in enumerate(s_red)}
669 # Initialize box
670 if self.parameters.create_box:
671 # count number of known types
672 n_types = len(self.parameters.atom_types)
673 create_box_command = 'create_box {} cell'.format(n_types)
674 self.lmp.command(create_box_command)
676 # Initialize the atoms with their types
677 # positions do not matter here
678 if self.parameters.create_atoms:
679 self.lmp.command('echo none') # don't echo the atom positions
680 self.rebuild(atoms)
681 self.lmp.command('echo log') # turn back on
682 else:
683 self.previous_atoms_numbers = atoms.numbers.copy()
685 # execute the user commands
686 for cmd in self.parameters.lmpcmds:
687 self.lmp.command(cmd)
689 # Set masses after user commands, e.g. to override
690 # EAM-provided masses
691 for sym in self.parameters.atom_types:
692 if self.parameters.atom_type_masses is None:
693 mass = ase_atomic_masses[ase_atomic_numbers[sym]]
694 else:
695 mass = self.parameters.atom_type_masses[sym]
696 self.lmp.command('mass %d %.30f' % (
697 self.parameters.atom_types[sym],
698 convert(mass, "mass", "ASE", self.units)))
700 # Define force & energy variables for extraction
701 self.lmp.command('variable pxx equal pxx')
702 self.lmp.command('variable pyy equal pyy')
703 self.lmp.command('variable pzz equal pzz')
704 self.lmp.command('variable pxy equal pxy')
705 self.lmp.command('variable pxz equal pxz')
706 self.lmp.command('variable pyz equal pyz')
708 # I am not sure why we need this next line but LAMMPS will
709 # raise an error if it is not there. Perhaps it is needed to
710 # ensure the cell stresses are calculated
711 self.lmp.command('thermo_style custom pe pxx emol ecoul')
713 self.lmp.command('variable fx atom fx')
714 self.lmp.command('variable fy atom fy')
715 self.lmp.command('variable fz atom fz')
717 # do we need this if we extract from a global ?
718 self.lmp.command('variable pe equal pe')
720 self.lmp.command("neigh_modify delay 0 every 1 check yes")
722 self.initialized = True