diff --git a/python/lammps.py b/python/lammps.py index a512efdcd..2db657fba 100644 --- a/python/lammps.py +++ b/python/lammps.py @@ -1,875 +1,932 @@ # ---------------------------------------------------------------------- # LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator # http://lammps.sandia.gov, Sandia National Laboratories # Steve Plimpton, sjplimp@sandia.gov # # Copyright (2003) Sandia Corporation. Under the terms of Contract # DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains # certain rights in this software. This software is distributed under # the GNU General Public License. # # See the README file in the top-level LAMMPS directory. # ------------------------------------------------------------------------- # Python wrappers on LAMMPS library via ctypes # for python3 compatibility from __future__ import print_function # imports for simple LAMMPS python wrapper module "lammps" import sys,traceback,types from ctypes import * from os.path import dirname,abspath,join from inspect import getsourcefile # imports for advanced LAMMPS python wrapper modules "PyLammps" and "IPyLammps" from collections import namedtuple import os import select import re import sys +def get_ctypes_int(size): + if size == 4: + return c_int32 + elif size == 8: + return c_int64 + return c_int + class MPIAbortException(Exception): def __init__(self, message): self.message = message def __str__(self): return repr(self.message) class lammps(object): # detect if Python is using version of mpi4py that can pass a communicator has_mpi4py_v2 = False try: from mpi4py import MPI from mpi4py import __version__ as mpi4py_version if mpi4py_version.split('.')[0] == '2': has_mpi4py_v2 = True except: pass # create instance of LAMMPS def __init__(self,name="",cmdargs=None,ptr=None,comm=None): self.comm = comm self.opened = 0 # determine module location modpath = dirname(abspath(getsourcefile(lambda:0))) self.lib = None # if a pointer to a LAMMPS object is handed in, # all symbols should already be available try: if ptr: self.lib = CDLL("",RTLD_GLOBAL) except: self.lib = None # load liblammps.so unless name is given # if name = "g++", load liblammps_g++.so # try loading the LAMMPS shared object from the location # of lammps.py with an absolute path, # so that LD_LIBRARY_PATH does not need to be set for regular install # fall back to loading with a relative path, # typically requires LD_LIBRARY_PATH to be set appropriately if not self.lib: try: if not name: self.lib = CDLL(join(modpath,"liblammps.so"),RTLD_GLOBAL) else: self.lib = CDLL(join(modpath,"liblammps_%s.so" % name), RTLD_GLOBAL) except: if not name: self.lib = CDLL("liblammps.so",RTLD_GLOBAL) else: self.lib = CDLL("liblammps_%s.so" % name,RTLD_GLOBAL) # if no ptr provided, create an instance of LAMMPS # don't know how to pass an MPI communicator from PyPar # but we can pass an MPI communicator from mpi4py v2.0.0 and later # no_mpi call lets LAMMPS use MPI_COMM_WORLD # cargs = array of C strings from args # if ptr, then are embedding Python in LAMMPS input script # ptr is the desired instance of LAMMPS # just convert it to ctypes ptr and store in self.lmp if not ptr: # with mpi4py v2, can pass MPI communicator to LAMMPS # need to adjust for type of MPI communicator object # allow for int (like MPICH) or void* (like OpenMPI) if lammps.has_mpi4py_v2 and comm != None: if lammps.MPI._sizeof(lammps.MPI.Comm) == sizeof(c_int): MPI_Comm = c_int else: MPI_Comm = c_void_p narg = 0 cargs = 0 if cmdargs: cmdargs.insert(0,"lammps.py") narg = len(cmdargs) for i in range(narg): if type(cmdargs[i]) is str: cmdargs[i] = cmdargs[i].encode() cargs = (c_char_p*narg)(*cmdargs) self.lib.lammps_open.argtypes = [c_int, c_char_p*narg, \ MPI_Comm, c_void_p()] else: self.lib.lammps_open.argtypes = [c_int, c_int, \ MPI_Comm, c_void_p()] self.lib.lammps_open.restype = None self.opened = 1 self.lmp = c_void_p() comm_ptr = lammps.MPI._addressof(comm) comm_val = MPI_Comm.from_address(comm_ptr) self.lib.lammps_open(narg,cargs,comm_val,byref(self.lmp)) else: self.opened = 1 if cmdargs: cmdargs.insert(0,"lammps.py") narg = len(cmdargs) for i in range(narg): if type(cmdargs[i]) is str: cmdargs[i] = cmdargs[i].encode() cargs = (c_char_p*narg)(*cmdargs) self.lmp = c_void_p() self.lib.lammps_open_no_mpi(narg,cargs,byref(self.lmp)) else: self.lmp = c_void_p() self.lib.lammps_open_no_mpi(0,None,byref(self.lmp)) # could use just this if LAMMPS lib interface supported it # self.lmp = self.lib.lammps_open_no_mpi(0,None) else: # magic to convert ptr to ctypes ptr if sys.version_info >= (3, 0): # Python 3 (uses PyCapsule API) pythonapi.PyCapsule_GetPointer.restype = c_void_p pythonapi.PyCapsule_GetPointer.argtypes = [py_object, c_char_p] self.lmp = c_void_p(pythonapi.PyCapsule_GetPointer(ptr, None)) else: # Python 2 (uses PyCObject API) pythonapi.PyCObject_AsVoidPtr.restype = c_void_p pythonapi.PyCObject_AsVoidPtr.argtypes = [py_object] self.lmp = c_void_p(pythonapi.PyCObject_AsVoidPtr(ptr)) + # optional numpy support (lazy loading) + self._numpy = None + + # set default types + self.c_bigint = get_ctypes_int(self.extract_setting("bigint")) + self.c_tagint = get_ctypes_int(self.extract_setting("tagint")) + self.c_imageint = get_ctypes_int(self.extract_setting("imageint")) + def __del__(self): if self.lmp and self.opened: self.lib.lammps_close(self.lmp) self.opened = 0 def close(self): if self.opened: self.lib.lammps_close(self.lmp) self.lmp = None self.opened = 0 def version(self): return self.lib.lammps_version(self.lmp) def file(self,file): if file: file = file.encode() self.lib.lammps_file(self.lmp,file) # send a single command def command(self,cmd): if cmd: cmd = cmd.encode() self.lib.lammps_command(self.lmp,cmd) if self.uses_exceptions and self.lib.lammps_has_error(self.lmp): sb = create_string_buffer(100) error_type = self.lib.lammps_get_last_error_message(self.lmp, sb, 100) error_msg = sb.value.decode().strip() if error_type == 2: raise MPIAbortException(error_msg) raise Exception(error_msg) # send a list of commands def commands_list(self,cmdlist): cmds = [x.encode() for x in cmdlist if type(x) is str] args = (c_char_p * len(cmdlist))(*cmds) self.lib.lammps_commands_list(self.lmp,len(cmdlist),args) # send a string of commands def commands_string(self,multicmd): if type(multicmd) is str: multicmd = multicmd.encode() self.lib.lammps_commands_string(self.lmp,c_char_p(multicmd)) # extract global info def extract_global(self,name,type): if name: name = name.encode() if type == 0: self.lib.lammps_extract_global.restype = POINTER(c_int) elif type == 1: self.lib.lammps_extract_global.restype = POINTER(c_double) else: return None ptr = self.lib.lammps_extract_global(self.lmp,name) return ptr[0] # extract per-atom info def extract_atom(self,name,type): if name: name = name.encode() if type == 0: self.lib.lammps_extract_atom.restype = POINTER(c_int) elif type == 1: self.lib.lammps_extract_atom.restype = POINTER(POINTER(c_int)) elif type == 2: self.lib.lammps_extract_atom.restype = POINTER(c_double) elif type == 3: self.lib.lammps_extract_atom.restype = POINTER(POINTER(c_double)) else: return None ptr = self.lib.lammps_extract_atom(self.lmp,name) return ptr + # extract lammps type byte sizes + + def extract_setting(self, name): + if name: name = name.encode() + self.lib.lammps_extract_atom.restype = c_int + return int(self.lib.lammps_extract_setting(self.lmp,name)) + + @property + def numpy(self): + if not self._numpy: + import numpy as np + class LammpsNumpyWrapper: + def __init__(self, lmp): + self.lmp = lmp + + def extract_atom_iarray(self, name, nelem, dim=1): + if dim == 1: + tmp = self.lmp.extract_atom(name, 0) + ptr = cast(tmp, POINTER(c_int * nelem)) + else: + tmp = self.lmp.extract_atom(name, 1) + ptr = cast(tmp[0], POINTER(c_int * nelem * dim)) + + a = np.frombuffer(ptr.contents, dtype=np.intc) + a.shape = (nelem, dim) + return a + + def extract_atom_darray(self, name, nelem, dim=1): + if dim == 1: + tmp = self.lmp.extract_atom(name, 2) + ptr = cast(tmp, POINTER(c_double * nelem)) + else: + tmp = self.lmp.extract_atom(name, 3) + ptr = cast(tmp[0], POINTER(c_double * nelem * dim)) + + a = np.frombuffer(ptr.contents) + a.shape = (nelem, dim) + return a + + self._numpy = LammpsNumpyWrapper(self) + return self._numpy + # extract compute info def extract_compute(self,id,style,type): if id: id = id.encode() if type == 0: if style > 0: return None self.lib.lammps_extract_compute.restype = POINTER(c_double) ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type) return ptr[0] if type == 1: self.lib.lammps_extract_compute.restype = POINTER(c_double) ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type) return ptr if type == 2: self.lib.lammps_extract_compute.restype = POINTER(POINTER(c_double)) ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type) return ptr return None # extract fix info # in case of global datum, free memory for 1 double via lammps_free() # double was allocated by library interface function def extract_fix(self,id,style,type,i=0,j=0): if id: id = id.encode() if style == 0: self.lib.lammps_extract_fix.restype = POINTER(c_double) ptr = self.lib.lammps_extract_fix(self.lmp,id,style,type,i,j) result = ptr[0] self.lib.lammps_free(ptr) return result elif (style == 1) or (style == 2): if type == 1: self.lib.lammps_extract_fix.restype = POINTER(c_double) elif type == 2: self.lib.lammps_extract_fix.restype = POINTER(POINTER(c_double)) else: return None ptr = self.lib.lammps_extract_fix(self.lmp,id,style,type,i,j) return ptr else: return None # extract variable info # free memory for 1 double or 1 vector of doubles via lammps_free() # for vector, must copy nlocal returned values to local c_double vector # memory was allocated by library interface function def extract_variable(self,name,group,type): if name: name = name.encode() if group: group = group.encode() if type == 0: self.lib.lammps_extract_variable.restype = POINTER(c_double) ptr = self.lib.lammps_extract_variable(self.lmp,name,group) result = ptr[0] self.lib.lammps_free(ptr) return result if type == 1: self.lib.lammps_extract_global.restype = POINTER(c_int) nlocalptr = self.lib.lammps_extract_global(self.lmp,"nlocal".encode()) nlocal = nlocalptr[0] result = (c_double*nlocal)() self.lib.lammps_extract_variable.restype = POINTER(c_double) ptr = self.lib.lammps_extract_variable(self.lmp,name,group) for i in range(nlocal): result[i] = ptr[i] self.lib.lammps_free(ptr) return result return None # set variable value # value is converted to string # returns 0 for success, -1 if failed def set_variable(self,name,value): if name: name = name.encode() if value: value = str(value).encode() return self.lib.lammps_set_variable(self.lmp,name,value) # return current value of thermo keyword def get_thermo(self,name): if name: name = name.encode() self.lib.lammps_get_thermo.restype = c_double return self.lib.lammps_get_thermo(self.lmp,name) # return total number of atoms in system def get_natoms(self): return self.lib.lammps_get_natoms(self.lmp) # return vector of atom properties gathered across procs, ordered by atom ID # name = atom property recognized by LAMMPS in atom->extract() # type = 0 for integer values, 1 for double values # count = number of per-atom valus, 1 for type or charge, 3 for x or f # returned data is a 1d vector - doc how it is ordered? # NOTE: how could we insure are converting to correct Python type # e.g. for Python list or NumPy, etc # ditto for extact_atom() above def gather_atoms(self,name,type,count): if name: name = name.encode() natoms = self.lib.lammps_get_natoms(self.lmp) if type == 0: data = ((count*natoms)*c_int)() self.lib.lammps_gather_atoms(self.lmp,name,type,count,data) elif type == 1: data = ((count*natoms)*c_double)() self.lib.lammps_gather_atoms(self.lmp,name,type,count,data) else: return None return data # scatter vector of atom properties across procs, ordered by atom ID # name = atom property recognized by LAMMPS in atom->extract() # type = 0 for integer values, 1 for double values # count = number of per-atom valus, 1 for type or charge, 3 for x or f # assume data is of correct type and length, as created by gather_atoms() # NOTE: how could we insure are passing correct type to LAMMPS # e.g. for Python list or NumPy, etc def scatter_atoms(self,name,type,count,data): if name: name = name.encode() self.lib.lammps_scatter_atoms(self.lmp,name,type,count,data) # create N atoms on all procs # N = global number of atoms # id = ID of each atom (optional, can be None) # type = type of each atom (1 to Ntypes) (required) # x = coords of each atom as (N,3) array (required) # v = velocity of each atom as (N,3) array (optional, can be None) # NOTE: how could we insure are passing correct type to LAMMPS # e.g. for Python list or NumPy, etc # ditto for gather_atoms() above def create_atoms(self,n,id,type,x,v,image=None,shrinkexceed=False): if id: id_lmp = (c_int * n)() id_lmp[:] = id else: id_lmp = id if image: image_lmp = (c_int * n)() image_lmp[:] = image else: image_lmp = image type_lmp = (c_int * n)() type_lmp[:] = type self.lib.lammps_create_atoms(self.lmp,n,id_lmp,type_lmp,x,v,image_lmp,shrinkexceed) @property def uses_exceptions(self): """ Return whether the LAMMPS shared library was compiled with C++ exceptions handling enabled """ try: if self.lib.lammps_has_error: return True except(AttributeError): return False # ------------------------------------------------------------------------- # ------------------------------------------------------------------------- # ------------------------------------------------------------------------- ################################################################################ # Alternative Python Wrapper # Written by Richard Berger ################################################################################ class OutputCapture(object): """ Utility class to capture LAMMPS library output """ def __init__(self): self.stdout_pipe_read, self.stdout_pipe_write = os.pipe() self.stdout_fd = 1 def __enter__(self): self.stdout = os.dup(self.stdout_fd) os.dup2(self.stdout_pipe_write, self.stdout_fd) return self def __exit__(self, type, value, tracebac): os.dup2(self.stdout, self.stdout_fd) os.close(self.stdout) os.close(self.stdout_pipe_read) os.close(self.stdout_pipe_write) # check if we have more to read from the pipe def more_data(self, pipe): r, _, _ = select.select([pipe], [], [], 0) return bool(r) # read the whole pipe def read_pipe(self, pipe): out = "" while self.more_data(pipe): out += os.read(pipe, 1024).decode() return out @property def output(self): return self.read_pipe(self.stdout_pipe_read) class Variable(object): def __init__(self, lammps_wrapper_instance, name, style, definition): self.wrapper = lammps_wrapper_instance self.name = name self.style = style self.definition = definition.split() @property def value(self): if self.style == 'atom': return list(self.wrapper.lmp.extract_variable(self.name, "all", 1)) else: value = self.wrapper.lmp_print('"${%s}"' % self.name).strip() try: return float(value) except ValueError: return value class AtomList(object): def __init__(self, lammps_wrapper_instance): self.lmp = lammps_wrapper_instance self.natoms = self.lmp.system.natoms self.dimensions = self.lmp.system.dimensions def __getitem__(self, index): if self.dimensions == 2: return Atom2D(self.lmp, index + 1) return Atom(self.lmp, index + 1) class Atom(object): def __init__(self, lammps_wrapper_instance, index): self.lmp = lammps_wrapper_instance self.index = index @property def id(self): return int(self.lmp.eval("id[%d]" % self.index)) @property def type(self): return int(self.lmp.eval("type[%d]" % self.index)) @property def mol(self): return self.lmp.eval("mol[%d]" % self.index) @property def mass(self): return self.lmp.eval("mass[%d]" % self.index) @property def position(self): return (self.lmp.eval("x[%d]" % self.index), self.lmp.eval("y[%d]" % self.index), self.lmp.eval("z[%d]" % self.index)) @position.setter def position(self, value): self.lmp.set("atom", self.index, "x", value[0]) self.lmp.set("atom", self.index, "y", value[1]) self.lmp.set("atom", self.index, "z", value[2]) @property def velocity(self): return (self.lmp.eval("vx[%d]" % self.index), self.lmp.eval("vy[%d]" % self.index), self.lmp.eval("vz[%d]" % self.index)) @property def force(self): return (self.lmp.eval("fx[%d]" % self.index), self.lmp.eval("fy[%d]" % self.index), self.lmp.eval("fz[%d]" % self.index)) @property def charge(self): return self.lmp.eval("q[%d]" % self.index) class Atom2D(Atom): def __init__(self, lammps_wrapper_instance, index): super(Atom2D, self).__init__(lammps_wrapper_instance, index) @property def position(self): return (self.lmp.eval("x[%d]" % self.index), self.lmp.eval("y[%d]" % self.index)) @position.setter def position(self, value): self.lmp.set("atom", self.index, "x", value[0]) self.lmp.set("atom", self.index, "y", value[1]) @property def velocity(self): return (self.lmp.eval("vx[%d]" % self.index), self.lmp.eval("vy[%d]" % self.index)) @property def force(self): return (self.lmp.eval("fx[%d]" % self.index), self.lmp.eval("fy[%d]" % self.index)) def get_thermo_data(output): """ traverse output of runs and extract thermo data columns """ if isinstance(output, str): lines = output.splitlines() else: lines = output runs = [] columns = [] in_run = False current_run = {} for line in lines: if line.startswith("Per MPI rank memory allocation"): in_run = True elif in_run and len(columns) == 0: # first line after memory usage are column names columns = line.split() current_run = {} for col in columns: current_run[col] = [] elif line.startswith("Loop time of "): in_run = False columns = None thermo_data = namedtuple('ThermoData', list(current_run.keys()))(*list(current_run.values())) r = {'thermo' : thermo_data } runs.append(namedtuple('Run', list(r.keys()))(*list(r.values()))) elif in_run and len(columns) > 0: values = [float(x) for x in line.split()] for i, col in enumerate(columns): current_run[col].append(values[i]) return runs class PyLammps(object): """ More Python-like wrapper for LAMMPS (e.g., for iPython) See examples/ipython for usage """ def __init__(self,name="",cmdargs=None,ptr=None,comm=None): if ptr: if isinstance(ptr,PyLammps): self.lmp = ptr.lmp elif isinstance(ptr,lammps): self.lmp = ptr else: self.lmp = lammps(name=name,cmdargs=cmdargs,ptr=ptr,comm=comm) else: self.lmp = lammps(name=name,cmdargs=cmdargs,ptr=None,comm=comm) print("LAMMPS output is captured by PyLammps wrapper") self._cmd_history = [] self.runs = [] def __del__(self): if self.lmp: self.lmp.close() self.lmp = None def close(self): if self.lmp: self.lmp.close() self.lmp = None def version(self): return self.lmp.version() def file(self,file): self.lmp.file(file) def write_script(self,filename): """ Write LAMMPS script file containing all commands executed up until now """ with open(filename, "w") as f: for cmd in self._cmd_history: f.write("%s\n" % cmd) def command(self,cmd): self.lmp.command(cmd) self._cmd_history.append(cmd) def run(self, *args, **kwargs): output = self.__getattr__('run')(*args, **kwargs) self.runs += get_thermo_data(output) return output @property def last_run(self): if len(self.runs) > 0: return self.runs[-1] return None @property def atoms(self): return AtomList(self) @property def system(self): output = self.info("system") d = self._parse_info_system(output) return namedtuple('System', d.keys())(*d.values()) @property def communication(self): output = self.info("communication") d = self._parse_info_communication(output) return namedtuple('Communication', d.keys())(*d.values()) @property def computes(self): output = self.info("computes") return self._parse_element_list(output) @property def dumps(self): output = self.info("dumps") return self._parse_element_list(output) @property def fixes(self): output = self.info("fixes") return self._parse_element_list(output) @property def groups(self): output = self.info("groups") return self._parse_groups(output) @property def variables(self): output = self.info("variables") vars = {} for v in self._parse_element_list(output): vars[v['name']] = Variable(self, v['name'], v['style'], v['def']) return vars def eval(self, expr): value = self.lmp_print('"$(%s)"' % expr).strip() try: return float(value) except ValueError: return value def _split_values(self, line): return [x.strip() for x in line.split(',')] def _get_pair(self, value): return [x.strip() for x in value.split('=')] def _parse_info_system(self, output): lines = output[6:-2] system = {} for line in lines: if line.startswith("Units"): system['units'] = self._get_pair(line)[1] elif line.startswith("Atom style"): system['atom_style'] = self._get_pair(line)[1] elif line.startswith("Atom map"): system['atom_map'] = self._get_pair(line)[1] elif line.startswith("Atoms"): parts = self._split_values(line) system['natoms'] = int(self._get_pair(parts[0])[1]) system['ntypes'] = int(self._get_pair(parts[1])[1]) system['style'] = self._get_pair(parts[2])[1] elif line.startswith("Kspace style"): system['kspace_style'] = self._get_pair(line)[1] elif line.startswith("Dimensions"): system['dimensions'] = int(self._get_pair(line)[1]) elif line.startswith("Orthogonal box"): system['orthogonal_box'] = [float(x) for x in self._get_pair(line)[1].split('x')] elif line.startswith("Boundaries"): system['boundaries'] = self._get_pair(line)[1] elif line.startswith("xlo"): keys, values = [self._split_values(x) for x in self._get_pair(line)] for key, value in zip(keys, values): system[key] = float(value) elif line.startswith("ylo"): keys, values = [self._split_values(x) for x in self._get_pair(line)] for key, value in zip(keys, values): system[key] = float(value) elif line.startswith("zlo"): keys, values = [self._split_values(x) for x in self._get_pair(line)] for key, value in zip(keys, values): system[key] = float(value) elif line.startswith("Molecule type"): system['molecule_type'] = self._get_pair(line)[1] elif line.startswith("Bonds"): parts = self._split_values(line) system['nbonds'] = int(self._get_pair(parts[0])[1]) system['nbondtypes'] = int(self._get_pair(parts[1])[1]) system['bond_style'] = self._get_pair(parts[2])[1] elif line.startswith("Angles"): parts = self._split_values(line) system['nangles'] = int(self._get_pair(parts[0])[1]) system['nangletypes'] = int(self._get_pair(parts[1])[1]) system['angle_style'] = self._get_pair(parts[2])[1] elif line.startswith("Dihedrals"): parts = self._split_values(line) system['ndihedrals'] = int(self._get_pair(parts[0])[1]) system['nangletypes'] = int(self._get_pair(parts[1])[1]) system['dihedral_style'] = self._get_pair(parts[2])[1] elif line.startswith("Impropers"): parts = self._split_values(line) system['nimpropers'] = int(self._get_pair(parts[0])[1]) system['nimpropertypes'] = int(self._get_pair(parts[1])[1]) system['improper_style'] = self._get_pair(parts[2])[1] return system def _parse_info_communication(self, output): lines = output[6:-3] comm = {} for line in lines: if line.startswith("MPI library"): comm['mpi_version'] = line.split(':')[1].strip() elif line.startswith("Comm style"): parts = self._split_values(line) comm['comm_style'] = self._get_pair(parts[0])[1] comm['comm_layout'] = self._get_pair(parts[1])[1] elif line.startswith("Processor grid"): comm['proc_grid'] = [int(x) for x in self._get_pair(line)[1].split('x')] elif line.startswith("Communicate velocities for ghost atoms"): comm['ghost_velocity'] = (self._get_pair(line)[1] == "yes") elif line.startswith("Nprocs"): parts = self._split_values(line) comm['nprocs'] = int(self._get_pair(parts[0])[1]) comm['nthreads'] = int(self._get_pair(parts[1])[1]) return comm def _parse_element_list(self, output): lines = output[6:-3] elements = [] for line in lines: element_info = self._split_values(line.split(':')[1].strip()) element = {'name': element_info[0]} for key, value in [self._get_pair(x) for x in element_info[1:]]: element[key] = value elements.append(element) return elements def _parse_groups(self, output): lines = output[6:-3] groups = [] group_pattern = re.compile(r"(?P.+) \((?P.+)\)") for line in lines: m = group_pattern.match(line.split(':')[1].strip()) group = {'name': m.group('name'), 'type': m.group('type')} groups.append(group) return groups def lmp_print(self, s): """ needed for Python2 compatibility, since print is a reserved keyword """ return self.__getattr__("print")(s) def __getattr__(self, name): def handler(*args, **kwargs): cmd_args = [name] + [str(x) for x in args] with OutputCapture() as capture: self.command(' '.join(cmd_args)) output = capture.output if 'verbose' in kwargs and kwargs['verbose']: print(output) lines = output.splitlines() if len(lines) > 1: return lines elif len(lines) == 1: return lines[0] return None return handler class IPyLammps(PyLammps): """ iPython wrapper for LAMMPS which adds embedded graphics capabilities """ def __init__(self,name="",cmdargs=None,ptr=None,comm=None): super(IPyLammps, self).__init__(name=name,cmdargs=cmdargs,ptr=ptr,comm=comm) def image(self, filename="snapshot.png", group="all", color="type", diameter="type", size=None, view=None, center=None, up=None, zoom=1.0): cmd_args = [group, "image", filename, color, diameter] if size: width = size[0] height = size[1] cmd_args += ["size", width, height] if view: theta = view[0] phi = view[1] cmd_args += ["view", theta, phi] if center: flag = center[0] Cx = center[1] Cy = center[2] Cz = center[3] cmd_args += ["center", flag, Cx, Cy, Cz] if up: Ux = up[0] Uy = up[1] Uz = up[2] cmd_args += ["up", Ux, Uy, Uz] if zoom: cmd_args += ["zoom", zoom] cmd_args.append("modify backcolor white") self.write_dump(*cmd_args) from IPython.core.display import Image return Image('snapshot.png') def video(self, filename): from IPython.display import HTML return HTML("") diff --git a/src/library.cpp b/src/library.cpp index 22b54f73a..8f87d6fab 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -1,1088 +1,1112 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ // C or Fortran style library interface to LAMMPS // customize by adding new LAMMPS-specific functions #include #include #include #include "library.h" #include "lmptype.h" #include "lammps.h" #include "universe.h" #include "input.h" #include "atom_vec.h" #include "atom.h" #include "domain.h" #include "update.h" #include "group.h" #include "input.h" #include "variable.h" #include "modify.h" #include "output.h" #include "thermo.h" #include "compute.h" #include "fix.h" #include "comm.h" #include "memory.h" #include "error.h" +#include "force.h" using namespace LAMMPS_NS; // ---------------------------------------------------------------------- // utility macros // ---------------------------------------------------------------------- /* ---------------------------------------------------------------------- macros for optional code path which captures all exceptions and stores the last error message. These assume there is a variable lmp which is a pointer to the current LAMMPS instance. Usage: BEGIN_CAPTURE { // code paths which might throw exception ... } END_CAPTURE ------------------------------------------------------------------------- */ #ifdef LAMMPS_EXCEPTIONS #define BEGIN_CAPTURE \ Error * error = lmp->error; \ try #define END_CAPTURE \ catch(LAMMPSAbortException & ae) { \ int nprocs = 0; \ MPI_Comm_size(ae.universe, &nprocs ); \ \ if (nprocs > 1) { \ error->set_last_error(ae.message.c_str(), ERROR_ABORT); \ } else { \ error->set_last_error(ae.message.c_str(), ERROR_NORMAL); \ } \ } catch(LAMMPSException & e) { \ error->set_last_error(e.message.c_str(), ERROR_NORMAL); \ } #else #define BEGIN_CAPTURE #define END_CAPTURE #endif // ---------------------------------------------------------------------- // helper functions, not in library API // ---------------------------------------------------------------------- /* ---------------------------------------------------------------------- concatenate one or more LAMMPS input lines starting at ptr removes NULL terminator when last printable char of line = '&' by replacing both NULL and '&' with space character repeat as many times as needed on return, ptr now points to longer line ------------------------------------------------------------------------- */ void concatenate_lines(char *ptr) { int nend = strlen(ptr); int n = nend-1; while (n && isspace(ptr[n])) n--; while (ptr[n] == '&') { ptr[nend] = ' '; ptr[n] = ' '; strtok(ptr,"\n"); nend = strlen(ptr); n = nend-1; while (n && isspace(ptr[n])) n--; } } // ---------------------------------------------------------------------- // library API functions to create/destroy an instance of LAMMPS // and communicate commands to it // ---------------------------------------------------------------------- /* ---------------------------------------------------------------------- create an instance of LAMMPS and return pointer to it pass in command-line args and MPI communicator to run on ------------------------------------------------------------------------- */ void lammps_open(int argc, char **argv, MPI_Comm communicator, void **ptr) { #ifdef LAMMPS_EXCEPTIONS try { LAMMPS *lmp = new LAMMPS(argc,argv,communicator); *ptr = (void *) lmp; } catch(LAMMPSException & e) { fprintf(stderr, "LAMMPS Exception: %s", e.message.c_str()); *ptr = (void *) NULL; } #else LAMMPS *lmp = new LAMMPS(argc,argv,communicator); *ptr = (void *) lmp; #endif } /* ---------------------------------------------------------------------- create an instance of LAMMPS and return pointer to it caller doesn't know MPI communicator, so use MPI_COMM_WORLD initialize MPI if needed ------------------------------------------------------------------------- */ void lammps_open_no_mpi(int argc, char **argv, void **ptr) { int flag; MPI_Initialized(&flag); if (!flag) { int argc = 0; char **argv = NULL; MPI_Init(&argc,&argv); } MPI_Comm communicator = MPI_COMM_WORLD; #ifdef LAMMPS_EXCEPTIONS try { LAMMPS *lmp = new LAMMPS(argc,argv,communicator); *ptr = (void *) lmp; } catch(LAMMPSException & e) { fprintf(stderr, "LAMMPS Exception: %s", e.message.c_str()); *ptr = (void*) NULL; } #else LAMMPS *lmp = new LAMMPS(argc,argv,communicator); *ptr = (void *) lmp; #endif } /* ---------------------------------------------------------------------- destruct an instance of LAMMPS ------------------------------------------------------------------------- */ void lammps_close(void *ptr) { LAMMPS *lmp = (LAMMPS *) ptr; delete lmp; } /* ---------------------------------------------------------------------- get the numerical representation of the current LAMMPS version ------------------------------------------------------------------------- */ int lammps_version(void *ptr) { LAMMPS *lmp = (LAMMPS *) ptr; return atoi(lmp->universe->num_ver); } /* ---------------------------------------------------------------------- process an input script in filename str ------------------------------------------------------------------------- */ void lammps_file(void *ptr, char *str) { LAMMPS *lmp = (LAMMPS *) ptr; BEGIN_CAPTURE { lmp->input->file(str); } END_CAPTURE } /* ---------------------------------------------------------------------- process a single input command in str does not matter if str ends in newline return command name to caller ------------------------------------------------------------------------- */ char *lammps_command(void *ptr, char *str) { LAMMPS *lmp = (LAMMPS *) ptr; char *result = NULL; BEGIN_CAPTURE { result = lmp->input->one(str); } END_CAPTURE return result; } /* ---------------------------------------------------------------------- process multiple input commands in cmds = list of strings does not matter if each string ends in newline create long contatentated string for processing by commands_string() insert newlines in concatenated string as needed ------------------------------------------------------------------------- */ void lammps_commands_list(void *ptr, int ncmd, char **cmds) { LAMMPS *lmp = (LAMMPS *) ptr; int n = ncmd+1; for (int i = 0; i < ncmd; i++) n += strlen(cmds[i]); char *str = (char *) lmp->memory->smalloc(n,"lib/commands/list:str"); str[0] = '\0'; n = 0; for (int i = 0; i < ncmd; i++) { strcpy(&str[n],cmds[i]); n += strlen(cmds[i]); if (str[n-1] != '\n') { str[n] = '\n'; str[n+1] = '\0'; n++; } } lammps_commands_string(ptr,str); lmp->memory->sfree(str); } /* ---------------------------------------------------------------------- process multiple input commands in single long str, separated by newlines single command can span multiple lines via continuation characters multi-line commands enabled by triple quotes will not work ------------------------------------------------------------------------- */ void lammps_commands_string(void *ptr, char *str) { LAMMPS *lmp = (LAMMPS *) ptr; // make copy of str so can strtok() it int n = strlen(str) + 1; char *copy = new char[n]; strcpy(copy,str); BEGIN_CAPTURE { char *ptr = copy; for (int i=0; i < n-1; ++i) { // handle continuation character as last character in line or string if ((copy[i] == '&') && (copy[i+1] == '\n')) copy[i+1] = copy[i] = ' '; else if ((copy[i] == '&') && (copy[i+1] == '\0')) copy[i] = ' '; if (copy[i] == '\n') { copy[i] = '\0'; lmp->input->one(ptr); ptr = copy + i+1; } else if (copy[i+1] == '\0') lmp->input->one(ptr); } } END_CAPTURE delete [] copy; } /* ---------------------------------------------------------------------- clean-up function to free memory allocated by lib and returned to caller ------------------------------------------------------------------------- */ void lammps_free(void *ptr) { free(ptr); } // ---------------------------------------------------------------------- // library API functions to extract info from LAMMPS or set info in LAMMPS // ---------------------------------------------------------------------- /* ---------------------------------------------------------------------- add LAMMPS-specific library functions all must receive LAMMPS pointer as argument customize by adding a function here and in library.h header file ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- extract a LAMMPS setting as an integer only use for settings that require return of an int customize by adding names ------------------------------------------------------------------------- */ int lammps_extract_setting(void *ptr, char *name) { if (strcmp(name,"bigint") == 0) return sizeof(bigint); if (strcmp(name,"tagint") == 0) return sizeof(tagint); if (strcmp(name,"imageint") == 0) return sizeof(imageint); return -1; } /* ---------------------------------------------------------------------- extract a pointer to an internal LAMMPS global entity name = desired quantity, e.g. dt or boxyhi or natoms returns a void pointer to the entity which the caller can cast to the proper data type returns a NULL if name not listed below this function need only be invoked once the returned pointer is a permanent valid reference to the quantity customize by adding names ------------------------------------------------------------------------- */ void *lammps_extract_global(void *ptr, char *name) { LAMMPS *lmp = (LAMMPS *) ptr; if (strcmp(name,"dt") == 0) return (void *) &lmp->update->dt; if (strcmp(name,"boxlo") == 0) return (void *) lmp->domain->boxlo; if (strcmp(name,"boxhi") == 0) return (void *) lmp->domain->boxhi; if (strcmp(name,"boxxlo") == 0) return (void *) &lmp->domain->boxlo[0]; if (strcmp(name,"boxxhi") == 0) return (void *) &lmp->domain->boxhi[0]; if (strcmp(name,"boxylo") == 0) return (void *) &lmp->domain->boxlo[1]; if (strcmp(name,"boxyhi") == 0) return (void *) &lmp->domain->boxhi[1]; if (strcmp(name,"boxzlo") == 0) return (void *) &lmp->domain->boxlo[2]; if (strcmp(name,"boxzhi") == 0) return (void *) &lmp->domain->boxhi[2]; if (strcmp(name,"periodicity") == 0) return (void *) lmp->domain->periodicity; if (strcmp(name,"xy") == 0) return (void *) &lmp->domain->xy; if (strcmp(name,"xz") == 0) return (void *) &lmp->domain->xz; if (strcmp(name,"yz") == 0) return (void *) &lmp->domain->yz; if (strcmp(name,"natoms") == 0) return (void *) &lmp->atom->natoms; if (strcmp(name,"nbonds") == 0) return (void *) &lmp->atom->nbonds; if (strcmp(name,"nangles") == 0) return (void *) &lmp->atom->nangles; if (strcmp(name,"ndihedrals") == 0) return (void *) &lmp->atom->ndihedrals; if (strcmp(name,"nimpropers") == 0) return (void *) &lmp->atom->nimpropers; if (strcmp(name,"nlocal") == 0) return (void *) &lmp->atom->nlocal; if (strcmp(name,"nghost") == 0) return (void *) &lmp->atom->nghost; if (strcmp(name,"nmax") == 0) return (void *) &lmp->atom->nmax; + if (strcmp(name,"ntypes") == 0) return (void *) &lmp->atom->ntypes; if (strcmp(name,"ntimestep") == 0) return (void *) &lmp->update->ntimestep; if (strcmp(name,"units") == 0) return (void *) lmp->update->unit_style; if (strcmp(name,"triclinic") == 0) return (void *) &lmp->domain->triclinic; if (strcmp(name,"q_flag") == 0) return (void *) &lmp->atom->q_flag; // update->atime can be referenced as a pointer // thermo "timer" data cannot be, since it is computed on request // lammps_get_thermo() can access all thermo keywords by value if (strcmp(name,"atime") == 0) return (void *) &lmp->update->atime; if (strcmp(name,"atimestep") == 0) return (void *) &lmp->update->atimestep; + // global constants defined by units + + if (strcmp(name,"boltz") == 0) return (void *) &lmp->force->boltz; + if (strcmp(name,"hplanck") == 0) return (void *) &lmp->force->hplanck; + if (strcmp(name,"mvv2e") == 0) return (void *) &lmp->force->mvv2e; + if (strcmp(name,"ftm2v") == 0) return (void *) &lmp->force->ftm2v; + if (strcmp(name,"mv2d") == 0) return (void *) &lmp->force->mv2d; + if (strcmp(name,"nktv2p") == 0) return (void *) &lmp->force->nktv2p; + if (strcmp(name,"qqr2e") == 0) return (void *) &lmp->force->qqr2e; + if (strcmp(name,"qe2f") == 0) return (void *) &lmp->force->qe2f; + if (strcmp(name,"vxmu2f") == 0) return (void *) &lmp->force->vxmu2f; + if (strcmp(name,"xxt2kmu") == 0) return (void *) &lmp->force->xxt2kmu; + if (strcmp(name,"dielectric") == 0) return (void *) &lmp->force->dielectric; + if (strcmp(name,"qqrd2e") == 0) return (void *) &lmp->force->qqrd2e; + if (strcmp(name,"e_mass") == 0) return (void *) &lmp->force->e_mass; + if (strcmp(name,"hhmrr2e") == 0) return (void *) &lmp->force->hhmrr2e; + if (strcmp(name,"mvh2r") == 0) return (void *) &lmp->force->mvh2r; + + if (strcmp(name,"angstrom") == 0) return (void *) &lmp->force->angstrom; + if (strcmp(name,"femtosecond") == 0) return (void *) &lmp->force->femtosecond; + if (strcmp(name,"qelectron") == 0) return (void *) &lmp->force->qelectron; + return NULL; } /* ---------------------------------------------------------------------- extract simulation box parameters see domain.h for definition of these arguments domain->init() call needed to set box_change ------------------------------------------------------------------------- */ void lammps_extract_box(void *ptr, double *boxlo, double *boxhi, double *xy, double *yz, double *xz, int *periodicity, int *box_change) { LAMMPS *lmp = (LAMMPS *) ptr; Domain *domain = lmp->domain; domain->init(); boxlo[0] = domain->boxlo[0]; boxlo[1] = domain->boxlo[1]; boxlo[2] = domain->boxlo[2]; boxhi[0] = domain->boxhi[0]; boxhi[1] = domain->boxhi[1]; boxhi[2] = domain->boxhi[2]; *xy = domain->xy; *yz = domain->yz; *xz = domain->xz; periodicity[0] = domain->periodicity[0]; periodicity[1] = domain->periodicity[1]; periodicity[2] = domain->periodicity[2]; *box_change = domain->box_change; } /* ---------------------------------------------------------------------- extract a pointer to an internal LAMMPS atom-based entity name = desired quantity, e.g. x or mass returns a void pointer to the entity which the caller can cast to the proper data type returns a NULL if Atom::extract() does not recognize the name the returned pointer is not a permanent valid reference to the per-atom quantity, since LAMMPS may reallocate per-atom data customize by adding names to Atom::extract() ------------------------------------------------------------------------- */ void *lammps_extract_atom(void *ptr, char *name) { LAMMPS *lmp = (LAMMPS *) ptr; return lmp->atom->extract(name); } /* ---------------------------------------------------------------------- extract a pointer to an internal LAMMPS compute-based entity the compute is invoked if its value(s) is not current id = compute ID style = 0 for global data, 1 for per-atom data, 2 for local data type = 0 for scalar, 1 for vector, 2 for array for global data, returns a pointer to the compute's internal data structure for the entity caller should cast it to (double *) for a scalar or vector caller should cast it to (double **) for an array for per-atom or local data, returns a pointer to the compute's internal data structure for the entity caller should cast it to (double *) for a vector caller should cast it to (double **) for an array returns a void pointer to the compute's internal data structure for the entity which the caller can cast to the proper data type returns a NULL if id is not recognized or style/type not supported the returned pointer is not a permanent valid reference to the compute data, this function should be re-invoked IMPORTANT: if the compute is not current it will be invoked LAMMPS cannot easily check here if it is valid to invoke the compute, so caller must insure that it is OK ------------------------------------------------------------------------- */ void *lammps_extract_compute(void *ptr, char *id, int style, int type) { LAMMPS *lmp = (LAMMPS *) ptr; BEGIN_CAPTURE { int icompute = lmp->modify->find_compute(id); if (icompute < 0) return NULL; Compute *compute = lmp->modify->compute[icompute]; if (style == 0) { if (type == 0) { if (!compute->scalar_flag) return NULL; if (compute->invoked_scalar != lmp->update->ntimestep) compute->compute_scalar(); return (void *) &compute->scalar; } if (type == 1) { if (!compute->vector_flag) return NULL; if (compute->invoked_vector != lmp->update->ntimestep) compute->compute_vector(); return (void *) compute->vector; } if (type == 2) { if (!compute->array_flag) return NULL; if (compute->invoked_array != lmp->update->ntimestep) compute->compute_array(); return (void *) compute->array; } } if (style == 1) { if (!compute->peratom_flag) return NULL; if (type == 1) { if (compute->invoked_peratom != lmp->update->ntimestep) compute->compute_peratom(); return (void *) compute->vector_atom; } if (type == 2) { if (compute->invoked_peratom != lmp->update->ntimestep) compute->compute_peratom(); return (void *) compute->array_atom; } } if (style == 2) { if (!compute->local_flag) return NULL; if (type == 1) { if (compute->invoked_local != lmp->update->ntimestep) compute->compute_local(); return (void *) compute->vector_local; } if (type == 2) { if (compute->invoked_local != lmp->update->ntimestep) compute->compute_local(); return (void *) compute->array_local; } } } END_CAPTURE return NULL; } /* ---------------------------------------------------------------------- extract a pointer to an internal LAMMPS fix-based entity id = fix ID style = 0 for global data, 1 for per-atom data, 2 for local data type = 0 for scalar, 1 for vector, 2 for array i,j = indices needed only to specify which global vector or array value for global data, returns a pointer to a memory location which is allocated by this function which the caller can cast to a (double *) which points to the value for per-atom or local data, returns a pointer to the fix's internal data structure for the entity caller should cast it to (double *) for a vector caller should cast it to (double **) for an array returns a NULL if id is not recognized or style/type not supported IMPORTANT: for global data, this function allocates a double to store the value in, so the caller must free this memory to avoid a leak, e.g. double *dptr = (double *) lammps_extract_fix(); double value = *dptr; lammps_free(dptr); IMPORTANT: LAMMPS cannot easily check here when info extracted from the fix is valid, so caller must insure that it is OK ------------------------------------------------------------------------- */ void *lammps_extract_fix(void *ptr, char *id, int style, int type, int i, int j) { LAMMPS *lmp = (LAMMPS *) ptr; BEGIN_CAPTURE { int ifix = lmp->modify->find_fix(id); if (ifix < 0) return NULL; Fix *fix = lmp->modify->fix[ifix]; if (style == 0) { double *dptr = (double *) malloc(sizeof(double)); if (type == 0) { if (!fix->scalar_flag) return NULL; *dptr = fix->compute_scalar(); return (void *) dptr; } if (type == 1) { if (!fix->vector_flag) return NULL; *dptr = fix->compute_vector(i); return (void *) dptr; } if (type == 2) { if (!fix->array_flag) return NULL; *dptr = fix->compute_array(i,j); return (void *) dptr; } } if (style == 1) { if (!fix->peratom_flag) return NULL; if (type == 1) return (void *) fix->vector_atom; if (type == 2) return (void *) fix->array_atom; } if (style == 2) { if (!fix->local_flag) return NULL; if (type == 1) return (void *) fix->vector_local; if (type == 2) return (void *) fix->array_local; } } END_CAPTURE return NULL; } /* ---------------------------------------------------------------------- extract a pointer to an internal LAMMPS evaluated variable name = variable name, must be equal-style or atom-style variable group = group ID for evaluating an atom-style variable, else NULL for equal-style variable, returns a pointer to a memory location which is allocated by this function which the caller can cast to a (double *) which points to the value for atom-style variable, returns a pointer to the vector of per-atom values on each processor, which the caller can cast to a (double *) which points to the values returns a NULL if name is not recognized or not equal-style or atom-style IMPORTANT: for both equal-style and atom-style variables, this function allocates memory to store the variable data in so the caller must free this memory to avoid a leak e.g. for equal-style variables double *dptr = (double *) lammps_extract_variable(); double value = *dptr; lammps_free(dptr); e.g. for atom-style variables double *vector = (double *) lammps_extract_variable(); use the vector values lammps_free(vector); IMPORTANT: LAMMPS cannot easily check here when it is valid to evaluate the variable or any fixes or computes or thermodynamic info it references, so caller must insure that it is OK ------------------------------------------------------------------------- */ void *lammps_extract_variable(void *ptr, char *name, char *group) { LAMMPS *lmp = (LAMMPS *) ptr; BEGIN_CAPTURE { int ivar = lmp->input->variable->find(name); if (ivar < 0) return NULL; if (lmp->input->variable->equalstyle(ivar)) { double *dptr = (double *) malloc(sizeof(double)); *dptr = lmp->input->variable->compute_equal(ivar); return (void *) dptr; } if (lmp->input->variable->atomstyle(ivar)) { int igroup = lmp->group->find(group); if (igroup < 0) return NULL; int nlocal = lmp->atom->nlocal; double *vector = (double *) malloc(nlocal*sizeof(double)); lmp->input->variable->compute_atom(ivar,igroup,vector,1,0); return (void *) vector; } } END_CAPTURE return NULL; } /* ---------------------------------------------------------------------- reset simulation box parameters see domain.h for definition of these arguments assumes domain->set_initial_box() has been invoked previously ------------------------------------------------------------------------- */ void lammps_reset_box(void *ptr, double *boxlo, double *boxhi, double xy, double yz, double xz) { LAMMPS *lmp = (LAMMPS *) ptr; Domain *domain = lmp->domain; domain->boxlo[0] = boxlo[0]; domain->boxlo[1] = boxlo[1]; domain->boxlo[2] = boxlo[2]; domain->boxhi[0] = boxhi[0]; domain->boxhi[1] = boxhi[1]; domain->boxhi[2] = boxhi[2]; domain->xy = xy; domain->yz = yz; domain->xz = xz; domain->set_global_box(); lmp->comm->set_proc_grid(); domain->set_local_box(); } /* ---------------------------------------------------------------------- set the value of a STRING variable to str return -1 if variable doesn't exist or not a STRING variable return 0 for success ------------------------------------------------------------------------- */ int lammps_set_variable(void *ptr, char *name, char *str) { LAMMPS *lmp = (LAMMPS *) ptr; int err = -1; BEGIN_CAPTURE { err = lmp->input->variable->set_string(name,str); } END_CAPTURE return err; } /* ---------------------------------------------------------------------- return the current value of a thermo keyword as a double unlike lammps_extract_global() this does not give access to the storage of the data in question instead it triggers the Thermo class to compute the current value and returns it ------------------------------------------------------------------------- */ double lammps_get_thermo(void *ptr, char *name) { LAMMPS *lmp = (LAMMPS *) ptr; double dval = 0.0; BEGIN_CAPTURE { lmp->output->thermo->evaluate_keyword(name,&dval); } END_CAPTURE return dval; } /* ---------------------------------------------------------------------- return the total number of atoms in the system useful before call to lammps_get_atoms() so can pre-allocate vector ------------------------------------------------------------------------- */ int lammps_get_natoms(void *ptr) { LAMMPS *lmp = (LAMMPS *) ptr; if (lmp->atom->natoms > MAXSMALLINT) return 0; int natoms = static_cast (lmp->atom->natoms); return natoms; } /* ---------------------------------------------------------------------- gather the named atom-based entity across all processors atom IDs must be consecutive from 1 to N name = desired quantity, e.g. x or charge type = 0 for integer values, 1 for double values count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f return atom-based values in 1d data, ordered by count, then by atom ID e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],... data must be pre-allocated by caller to correct length ------------------------------------------------------------------------- */ void lammps_gather_atoms(void *ptr, char *name, int type, int count, void *data) { LAMMPS *lmp = (LAMMPS *) ptr; BEGIN_CAPTURE { // error if tags are not defined or not consecutive int flag = 0; if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1; if (lmp->atom->natoms > MAXSMALLINT) flag = 1; if (flag) { if (lmp->comm->me == 0) lmp->error->warning(FLERR,"Library error in lammps_gather_atoms"); return; } int natoms = static_cast (lmp->atom->natoms); int i,j,offset; void *vptr = lmp->atom->extract(name); if(vptr == NULL) { lmp->error->warning(FLERR,"lammps_gather_atoms: unknown property name"); return; } // copy = Natom length vector of per-atom values // use atom ID to insert each atom's values into copy // MPI_Allreduce with MPI_SUM to merge into data, ordered by atom ID if (type == 0) { int *vector = NULL; int **array = NULL; const int imgunpack = (count == 3) && (strcmp(name,"image") == 0); if ((count == 1) || imgunpack) vector = (int *) vptr; else array = (int **) vptr; int *copy; lmp->memory->create(copy,count*natoms,"lib/gather:copy"); for (i = 0; i < count*natoms; i++) copy[i] = 0; tagint *tag = lmp->atom->tag; int nlocal = lmp->atom->nlocal; if (count == 1) for (i = 0; i < nlocal; i++) copy[tag[i]-1] = vector[i]; else if (imgunpack) { for (i = 0; i < nlocal; i++) { offset = count*(tag[i]-1); const int image = vector[i]; copy[offset++] = (image & IMGMASK) - IMGMAX; copy[offset++] = ((image >> IMGBITS) & IMGMASK) - IMGMAX; copy[offset++] = ((image >> IMG2BITS) & IMGMASK) - IMGMAX; } } else for (i = 0; i < nlocal; i++) { offset = count*(tag[i]-1); for (j = 0; j < count; j++) copy[offset++] = array[i][j]; } MPI_Allreduce(copy,data,count*natoms,MPI_INT,MPI_SUM,lmp->world); lmp->memory->destroy(copy); } else { double *vector = NULL; double **array = NULL; if (count == 1) vector = (double *) vptr; else array = (double **) vptr; double *copy; lmp->memory->create(copy,count*natoms,"lib/gather:copy"); for (i = 0; i < count*natoms; i++) copy[i] = 0.0; tagint *tag = lmp->atom->tag; int nlocal = lmp->atom->nlocal; if (count == 1) { for (i = 0; i < nlocal; i++) copy[tag[i]-1] = vector[i]; } else { for (i = 0; i < nlocal; i++) { offset = count*(tag[i]-1); for (j = 0; j < count; j++) copy[offset++] = array[i][j]; } } MPI_Allreduce(copy,data,count*natoms,MPI_DOUBLE,MPI_SUM,lmp->world); lmp->memory->destroy(copy); } } END_CAPTURE } /* ---------------------------------------------------------------------- scatter the named atom-based entity across all processors atom IDs must be consecutive from 1 to N name = desired quantity, e.g. x or charge type = 0 for integer values, 1 for double values count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f data = atom-based values in 1d data, ordered by count, then by atom ID e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],... ------------------------------------------------------------------------- */ void lammps_scatter_atoms(void *ptr, char *name, int type, int count, void *data) { LAMMPS *lmp = (LAMMPS *) ptr; BEGIN_CAPTURE { // error if tags are not defined or not consecutive or no atom map int flag = 0; if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1; if (lmp->atom->natoms > MAXSMALLINT) flag = 1; if (lmp->atom->map_style == 0) flag = 1; if (flag) { if (lmp->comm->me == 0) lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms"); return; } int natoms = static_cast (lmp->atom->natoms); int i,j,m,offset; void *vptr = lmp->atom->extract(name); if(vptr == NULL) { lmp->error->warning(FLERR, "lammps_scatter_atoms: unknown property name"); return; } // copy = Natom length vector of per-atom values // use atom ID to insert each atom's values into copy // MPI_Allreduce with MPI_SUM to merge into data, ordered by atom ID if (type == 0) { int *vector = NULL; int **array = NULL; const int imgpack = (count == 3) && (strcmp(name,"image") == 0); if ((count == 1) || imgpack) vector = (int *) vptr; else array = (int **) vptr; int *dptr = (int *) data; if (count == 1) { for (i = 0; i < natoms; i++) if ((m = lmp->atom->map(i+1)) >= 0) vector[m] = dptr[i]; } else if (imgpack) { for (i = 0; i < natoms; i++) if ((m = lmp->atom->map(i+1)) >= 0) { offset = count*i; int image = dptr[offset++] + IMGMAX; image += (dptr[offset++] + IMGMAX) << IMGBITS; image += (dptr[offset++] + IMGMAX) << IMG2BITS; vector[m] = image; } } else { for (i = 0; i < natoms; i++) if ((m = lmp->atom->map(i+1)) >= 0) { offset = count*i; for (j = 0; j < count; j++) array[m][j] = dptr[offset++]; } } } else { double *vector = NULL; double **array = NULL; if (count == 1) vector = (double *) vptr; else array = (double **) vptr; double *dptr = (double *) data; if (count == 1) { for (i = 0; i < natoms; i++) if ((m = lmp->atom->map(i+1)) >= 0) vector[m] = dptr[i]; } else { for (i = 0; i < natoms; i++) { if ((m = lmp->atom->map(i+1)) >= 0) { offset = count*i; for (j = 0; j < count; j++) array[m][j] = dptr[offset++]; } } } } } END_CAPTURE } /* ---------------------------------------------------------------------- create N atoms and assign them to procs based on coords id = atom IDs (optional, NULL will generate 1 to N) type = N-length vector of atom types (required) x = 3N-length 1d vector of atom coords (required) v = 3N-length 1d vector of atom velocities (optional, NULL if just 0.0) image flags can be treated in two ways: (a) image = vector of current image flags each atom will be remapped into periodic box by domain->ownatom() image flag will be incremented accordingly and stored with atom (b) image = NULL each atom will be remapped into periodic box by domain->ownatom() image flag will be set to 0 by atom->avec->create_atom() shrinkexceed = 1 allows atoms to be outside a shrinkwrapped boundary passed to ownatom() which will assign them to boundary proc important if atoms may be (slightly) outside non-periodic dim e.g. due to restoring a snapshot from a previous run and previous box id and image must be 32-bit integers x,v = ordered by xyz, then by atom e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],... ------------------------------------------------------------------------- */ void lammps_create_atoms(void *ptr, int n, tagint *id, int *type, double *x, double *v, imageint *image, int shrinkexceed) { LAMMPS *lmp = (LAMMPS *) ptr; BEGIN_CAPTURE { // error if box does not exist or tags not defined int flag = 0; if (lmp->domain->box_exist == 0) flag = 1; if (lmp->atom->tag_enable == 0) flag = 1; if (flag) { if (lmp->comm->me == 0) lmp->error->warning(FLERR,"Library error in lammps_create_atoms"); return; } // loop over N atoms of entire system // if this proc owns it based on coords, invoke create_atom() // optionally set atom tags and velocities Atom *atom = lmp->atom; Domain *domain = lmp->domain; int nlocal = atom->nlocal; bigint natoms_prev = atom->natoms; int nlocal_prev = nlocal; double xdata[3]; for (int i = 0; i < n; i++) { xdata[0] = x[3*i]; xdata[1] = x[3*i+1]; xdata[2] = x[3*i+2]; imageint * img = image ? &image[i] : NULL; tagint tag = id ? id[i] : -1; if (!domain->ownatom(tag, xdata, img, shrinkexceed)) continue; atom->avec->create_atom(type[i],xdata); if (id) atom->tag[nlocal] = id[i]; else atom->tag[nlocal] = i+1; if (v) { atom->v[nlocal][0] = v[3*i]; atom->v[nlocal][1] = v[3*i+1]; atom->v[nlocal][2] = v[3*i+2]; } if (image) atom->image[nlocal] = image[i]; nlocal++; } // need to reset atom->natoms inside LAMMPS bigint ncurrent = nlocal; MPI_Allreduce(&ncurrent,&lmp->atom->natoms,1,MPI_LMP_BIGINT, MPI_SUM,lmp->world); // init per-atom fix/compute/variable values for created atoms atom->data_fix_compute_variable(nlocal_prev,nlocal); // if global map exists, reset it // invoke map_init() b/c atom count has grown if (lmp->atom->map_style) { lmp->atom->map_init(); lmp->atom->map_set(); } // warn if new natoms is not correct if (lmp->atom->natoms != natoms_prev + n) { char str[128]; sprintf(str,"Library warning in lammps_create_atoms, " "invalid total atoms %ld %ld",lmp->atom->natoms,natoms_prev+n); if (lmp->comm->me == 0) lmp->error->warning(FLERR,str); } } END_CAPTURE } // ---------------------------------------------------------------------- // library API functions for error handling // ---------------------------------------------------------------------- #ifdef LAMMPS_EXCEPTIONS /* ---------------------------------------------------------------------- check if a new error message ------------------------------------------------------------------------- */ int lammps_has_error(void *ptr) { LAMMPS * lmp = (LAMMPS *) ptr; Error * error = lmp->error; return error->get_last_error() ? 1 : 0; } /* ---------------------------------------------------------------------- copy the last error message of LAMMPS into a character buffer return value encodes which type of error: 1 = normal error (recoverable) 2 = abort error (non-recoverable) ------------------------------------------------------------------------- */ int lammps_get_last_error_message(void *ptr, char * buffer, int buffer_size) { LAMMPS * lmp = (LAMMPS *) ptr; Error * error = lmp->error; if(error->get_last_error()) { int error_type = error->get_last_error_type(); strncpy(buffer, error->get_last_error(), buffer_size-1); error->set_last_error(NULL, ERROR_NONE); return error_type; } return 0; } #endif