Merge pull request #1674 from rbberger/library_interface_update

Python script, as follows:

.. parsed-literal::

.. code-block:: Python

   from lammps import lammps

can run LAMMPS, grab data, change it, and put it back into LAMMPS.

.. parsed-literal::

.. code-block:: Python

   lmp = lammps()           # create a LAMMPS object using the default liblammps.so library

                            # 4 optional args are allowed: name, cmdargs, ptr, comm

The lines

.. parsed-literal::

.. code-block:: Python

   from lammps import lammps

   lmp = lammps()

If the ptr argument is set like this:

.. parsed-literal::

.. code-block:: Python

   lmp = lammps(ptr=lmpptr)

script, and coordinate and run multiple simulations, e.g.

.. parsed-literal::

.. code-block:: Python

   from lammps import lammps

   lmp1 = lammps()

like this:

.. parsed-literal::

.. code-block:: Python

   from ctypes import \*

   natoms = lmp.get_natoms()

  Python script.

----------

.. autoclass:: lammps.lammps

   :members:

   :no-undoc-members:

.. autoclass:: lammps.NeighList

   :members:

   :no-undoc-members:

.. _lws: http://lammps.sandia.gov

.. _ld: Manual.html

"Higher level section"_Python_head.html - "LAMMPS WWW Site"_lws -

"LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c

:link(lws,http://lammps.sandia.gov)

:link(ld,Manual.html)

:link(lc,Commands_all.html)

:line

Python library interface :h3

As described previously, the Python interface to LAMMPS consists of a

Python "lammps" module, the source code for which is in

python/lammps.py, which creates a "lammps" object, with a set of

methods that can be invoked on that object.  The sample Python code

below assumes you have first imported the "lammps" module in your

Python script, as follows:

from lammps import lammps :pre

These are the methods defined by the lammps module.  If you look at

the files src/library.cpp and src/library.h you will see they

correspond one-to-one with calls you can make to the LAMMPS library

from a C++ or C or Fortran program, and which are described on the

"Howto library"_Howto_library.html doc page.

The python/examples directory has Python scripts which show how Python

can run LAMMPS, grab data, change it, and put it back into LAMMPS.

lmp = lammps()           # create a LAMMPS object using the default liblammps.so library

                         # 4 optional args are allowed: name, cmdargs, ptr, comm

lmp = lammps(ptr=lmpptr) # use lmpptr as previously created LAMMPS object

lmp = lammps(comm=split) # create a LAMMPS object with a custom communicator, requires mpi4py 2.0.0 or later

lmp = lammps(name="g++")   # create a LAMMPS object using the liblammps_g++.so library

lmp = lammps(name="g++",cmdargs=list)    # add LAMMPS command-line args, e.g. list = \["-echo","screen"\] :pre

lmp.close()              # destroy a LAMMPS object :pre

version = lmp.version()  # return the numerical version id, e.g. LAMMPS 2 Sep 2015 -> 20150902 :pre

lmp.file(file)           # run an entire input script, file = "in.lj"

lmp.command(cmd)         # invoke a single LAMMPS command, cmd = "run 100"

lmp.commands_list(cmdlist)     # invoke commands in cmdlist = ["run 10", "run 20"]

lmp.commands_string(multicmd)  # invoke commands in multicmd = "run 10\nrun 20" :pre

size = lmp.extract_setting(name)     # return data type info :pre

xlo = lmp.extract_global(name,type)  # extract a global quantity

                                     # name = "boxxlo", "nlocal", etc

                                     # type = 0 = int

                                     #        1 = double :pre

boxlo,boxhi,xy,yz,xz,periodicity,box_change = lmp.extract_box()  # extract box info :pre

coords = lmp.extract_atom(name,type)      # extract a per-atom quantity

                                          # name = "x", "type", etc

                                          # type = 0 = vector of ints

                                          #        1 = array of ints

                                          #        2 = vector of doubles

                                          #        3 = array of doubles :pre

eng = lmp.extract_compute(id,style,type)  # extract value(s) from a compute

v3 = lmp.extract_fix(id,style,type,i,j)   # extract value(s) from a fix

                                          # id = ID of compute or fix

                                          # style = 0 = global data

                                          #         1 = per-atom data

                                          #         2 = local data

                                          # type = 0 = scalar

                                          #        1 = vector

                                          #        2 = array

                                          # i,j = indices of value in global vector or array :pre

var = lmp.extract_variable(name,group,flag)  # extract value(s) from a variable

                                             # name = name of variable

                                             # group = group ID (ignored for equal-style variables)

                                             # flag = 0 = equal-style variable

                                             #        1 = atom-style variable :pre

value = lmp.get_thermo(name)              # return current value of a thermo keyword

natoms = lmp.get_natoms()                 # total # of atoms as int :pre

flag = lmp.set_variable(name,value)       # set existing named string-style variable to value, flag = 0 if successful

lmp.reset_box(boxlo,boxhi,xy,yz,xz)       # reset the simulation box size :pre

data = lmp.gather_atoms(name,type,count)  # return per-atom property of all atoms gathered into data, ordered by atom ID

                                          # name = "x", "charge", "type", etc

data = lmp.gather_atoms_concat(name,type,count)  # ditto, but concatenated atom values from each proc (unordered)

data = lmp.gather_atoms_subset(name,type,count,ndata,ids)  # ditto, but for subset of Ndata atoms with IDs :pre

lmp.scatter_atoms(name,type,count,data)   # scatter per-atom property to all atoms from data, ordered by atom ID

                                          # name = "x", "charge", "type", etc

                                          # count = # of per-atom values, 1 or 3, etc :pre

lmp.scatter_atoms_subset(name,type,count,ndata,ids,data)  # ditto, but for subset of Ndata atoms with IDs :pre

lmp.create_atoms(n,ids,types,x,v,image,shrinkexceed)   # create N atoms with IDs, types, x, v, and image flags :pre

:line

The lines

from lammps import lammps

lmp = lammps() :pre

create an instance of LAMMPS, wrapped in a Python class by the lammps

Python module, and return an instance of the Python class as lmp.  It

is used to make all subsequent calls to the LAMMPS library.

Additional arguments to lammps() can be used to tell Python the name

of the shared library to load or to pass arguments to the LAMMPS

instance, the same as if LAMMPS were launched from a command-line

prompt.

If the ptr argument is set like this:

lmp = lammps(ptr=lmpptr) :pre

then lmpptr must be an argument passed to Python via the LAMMPS

"python"_python.html command, when it is used to define a Python

function that is invoked by the LAMMPS input script.  This mode of

calling Python from LAMMPS is described in the "Python

call"_Python_call.html doc page.  The variable lmpptr refers to the

instance of LAMMPS that called the embedded Python interpreter.  Using

it as an argument to lammps() allows the returned Python class

instance "lmp" to make calls to that instance of LAMMPS.  See the

"python"_python.html command doc page for examples using this syntax.

Note that you can create multiple LAMMPS objects in your Python

script, and coordinate and run multiple simulations, e.g.

from lammps import lammps

lmp1 = lammps()

lmp2 = lammps()

lmp1.file("in.file1")

lmp2.file("in.file2") :pre

The file(), command(), commands_list(), commands_string() methods

allow an input script, a single command, or multiple commands to be

invoked.

The extract_setting(), extract_global(), extract_box(),

extract_atom(), extract_compute(), extract_fix(), and

extract_variable() methods return values or pointers to data

structures internal to LAMMPS.

For extract_global() see the src/library.cpp file for the list of

valid names.  New names could easily be added.  A double or integer is

returned.  You need to specify the appropriate data type via the type

argument.

For extract_atom(), a pointer to internal LAMMPS atom-based data is

returned, which you can use via normal Python subscripting.  See the

extract() method in the src/atom.cpp file for a list of valid names.

Again, new names could easily be added if the property you want is not

listed.  A pointer to a vector of doubles or integers, or a pointer to

an array of doubles (double **) or integers (int **) is returned.  You

need to specify the appropriate data type via the type argument.

For extract_compute() and extract_fix(), the global, per-atom, or

local data calculated by the compute or fix can be accessed.  What is

returned depends on whether the compute or fix calculates a scalar or

vector or array.  For a scalar, a single double value is returned.  If

the compute or fix calculates a vector or array, a pointer to the

internal LAMMPS data is returned, which you can use via normal Python

subscripting.  The one exception is that for a fix that calculates a

global vector or array, a single double value from the vector or array

is returned, indexed by I (vector) or I and J (array).  I,J are

zero-based indices.  The I,J arguments can be left out if not needed.

See the "Howto output"_Howto_output.html doc page for a discussion of

global, per-atom, and local data, and of scalar, vector, and array

data types.  See the doc pages for individual "computes"_compute.html

and "fixes"_fix.html for a description of what they calculate and

store.

For extract_variable(), an "equal-style or atom-style

variable"_variable.html is evaluated and its result returned.

For equal-style variables a single double value is returned and the

group argument is ignored.  For atom-style variables, a vector of

doubles is returned, one value per atom, which you can use via normal

Python subscripting. The values will be zero for atoms not in the

specified group.

The get_thermo() method returns the current value of a thermo

keyword as a float.

The get_natoms() method returns the total number of atoms in the

simulation, as an int.

The set_variable() method sets an existing string-style variable to a

new string value, so that subsequent LAMMPS commands can access the

variable.

The reset_box() method resets the size and shape of the simulation

box, e.g. as part of restoring a previously extracted and saved state

of a simulation.

The gather methods collect peratom info of the requested type (atom

coords, atom types, forces, etc) from all processors, and returns the

same vector of values to each calling processor.  The scatter

functions do the inverse.  They distribute a vector of peratom values,

passed by all calling processors, to individual atoms, which may be

owned by different processors.

Note that the data returned by the gather methods,

e.g. gather_atoms("x"), is different from the data structure returned

by extract_atom("x") in four ways.  (1) Gather_atoms() returns a

vector which you index as x\[i\]; extract_atom() returns an array

which you index as x\[i\]\[j\].  (2) Gather_atoms() orders the atoms

by atom ID while extract_atom() does not.  (3) Gather_atoms() returns

a list of all atoms in the simulation; extract_atoms() returns just

the atoms local to each processor.  (4) Finally, the gather_atoms()

data structure is a copy of the atom coords stored internally in

LAMMPS, whereas extract_atom() returns an array that effectively

points directly to the internal data.  This means you can change

values inside LAMMPS from Python by assigning a new values to the

extract_atom() array.  To do this with the gather_atoms() vector, you

need to change values in the vector, then invoke the scatter_atoms()

method.

For the scatter methods, the array of coordinates passed to must be a

ctypes vector of ints or doubles, allocated and initialized something

like this:

from ctypes import *

natoms = lmp.get_natoms()

n3 = 3*natoms

x = (n3*c_double)()

x\[0\] = x coord of atom with ID 1

x\[1\] = y coord of atom with ID 1

x\[2\] = z coord of atom with ID 1

x\[3\] = x coord of atom with ID 2

...

x\[n3-1\] = z coord of atom with ID natoms

lmp.scatter_atoms("x",1,3,x) :pre

Alternatively, you can just change values in the vector returned by

the gather methods, since they are also ctypes vectors.

:line

As noted above, these Python class methods correspond one-to-one with

the functions in the LAMMPS library interface in src/library.cpp and

library.h.  This means you can extend the Python wrapper via the

following steps:

Add a new interface function to src/library.cpp and

src/library.h. :ulb,l

Rebuild LAMMPS as a shared library. :l

Add a wrapper method to python/lammps.py for this interface

function. :l

You should now be able to invoke the new interface function from a

Python script. :l

:ule

extensions = [

    'sphinx.ext.mathjax',

    'sphinx.ext.imgmath',

    'sphinx.ext.autodoc',

    'table_from_list',

]

# 2017-12-07: commented out, since this package is broken with Sphinx 16.x

from sphinx.highlighting import lexers

lexers['LAMMPS'] = LAMMPSLexer.LAMMPSLexer(startinline=True)

sys.path.append(os.path.join(os.path.dirname(__file__), '../../../python'))

# 3d Lennard-Jones melt

units		lj

atom_style	atomic

lattice		fcc 0.8442

region		box block 0 2 0 2 0 2

create_box	1 box

create_atoms	1 box

mass		1 1.0

velocity	all create 3.0 87287

pair_style	lj/cut 2.5

pair_coeff	1 1 1.0 1.0 2.5

neighbor	0.1 bin

neigh_modify	every 20 delay 0 check no

python         post_force_callback here """

from __future__ import print_function

from lammps import lammps

def post_force_callback(lmp, v):

  try:

    import os

    pid = os.getpid()

    pid_prefix = "[{}] ".format(pid)

    L = lammps(ptr=lmp)

    t = L.extract_global("ntimestep", 0)

    print(pid_prefix, "### POST_FORCE ###", t)

    #mylist = L.get_neighlist(0)

    mylist = L.find_pair_neighlist("lj/cut", request=0)

    print(pid_prefix, mylist)

    nlocal = L.extract_global("nlocal", 0)

    nghost = L.extract_global("nghost", 0)

    ntypes = L.extract_global("ntypes", 0)

    mass = L.numpy.extract_atom_darray("mass", ntypes+1)

    atype = L.numpy.extract_atom_iarray("type", nlocal+nghost)

    x = L.numpy.extract_atom_darray("x", nlocal+nghost, dim=3)

    v = L.numpy.extract_atom_darray("v", nlocal+nghost, dim=3)

    f = L.numpy.extract_atom_darray("f", nlocal+nghost, dim=3)

    for iatom, numneigh, neighs in mylist:

      print(pid_prefix, "- {}".format(iatom), x[iatom], v[iatom], f[iatom], " : ",  numneigh, "Neighbors")

      for jatom in neighs:

        if jatom < nlocal:

            print(pid_prefix, "    *  ", jatom, x[jatom], v[jatom], f[jatom])

        else:

            print(pid_prefix, "    * [GHOST]", jatom, x[jatom], v[jatom], f[jatom])

  except Exception as e:

    print(e)

"""

fix		1 all nve

fix     3 all python/invoke 50 post_force post_force_callback

#dump		id all atom 1 dump.melt

#dump		2 all image 1 image.*.jpg type type &

#		axes yes 0.8 0.02 view 60 -30

#dump_modify	2 pad 3

#dump		3 all movie 1 movie.mpg type type &

#		axes yes 0.8 0.02 view 60 -30

#dump_modify	3 pad 3

thermo		1

run		100

  def __str__(self):

    return repr(self.message)

class NeighList:

    """This is a wrapper class that exposes the contents of a neighbor list

    It can be used like a regular Python list.

    Internally it uses the lower-level LAMMPS C-library interface.

    :param lmp: reference to instance of :class:`lammps`

    :type  lmp: lammps

    :param idx: neighbor list index

    :type  idx: int

    """

    def __init__(self, lmp, idx):

        self.lmp = lmp

        self.idx = idx

    def __str__(self):

        return "Neighbor List ({} atoms)".format(self.size)

    def __repr__(self):

        return self.__str__()

    @property

    def size(self):

        """

        :return: number of elements in neighbor list

        """

        return self.lmp.get_neighlist_size(self.idx)

    def get(self, element):

        """

        :return: tuple with atom local index, number of neighbors and array of neighbor local atom indices

        :rtype:  (int, int, numpy.array)

        """

        iatom, numneigh, neighbors = self.lmp.get_neighlist_element_neighbors(self.idx, element)

        return iatom, numneigh, neighbors

    # the methods below implement the iterator interface, so NeighList can be used like a regular Python list

    def __getitem__(self, element):

        return self.get(element)

    def __len__(self):

        return self.size

    def __iter__(self):

        inum = self.size

        for ii in range(inum):

            yield self.get(ii)

class lammps(object):

  # detect if Python is using version of mpi4py that can pass a communicator

    modpath = dirname(abspath(getsourcefile(lambda:0)))

    self.lib = None

    self.lmp = None

    # if a pointer to a LAMMPS object is handed in,

    # all symbols should already be available

      [c_void_p,c_char_p,c_int,c_int,c_int,POINTER(c_int),c_void_p]

    self.lib.lammps_scatter_atoms_subset.restype = None

    self.lib.lammps_find_pair_neighlist.argtypes = [c_void_p, c_char_p, c_int, c_int, c_int]

    self.lib.lammps_find_pair_neighlist.restype  = c_int

    self.lib.lammps_find_fix_neighlist.argtypes = [c_void_p, c_char_p, c_int]

    self.lib.lammps_find_fix_neighlist.restype  = c_int

    self.lib.lammps_find_compute_neighlist.argtypes = [c_void_p, c_char_p, c_int]

    self.lib.lammps_find_compute_neighlist.restype  = c_int

    self.lib.lammps_neighlist_num_elements.argtypes = [c_void_p, c_int]

    self.lib.lammps_neighlist_num_elements.restype  = c_int

    self.lib.lammps_neighlist_element_neighbors.argtypes = [c_void_p, c_int, c_int, POINTER(c_int), POINTER(c_int), POINTER(POINTER(c_int))]

    self.lib.lammps_neighlist_element_neighbors.restype  = None

    # 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

    self.lib.lammps_set_fix_external_callback(self.lmp, fix_name.encode(), cFunc, cCaller)

  def get_neighlist(self, idx):

    """Returns an instance of :class:`NeighList` which wraps access to the neighbor list with the given index

    :param idx: index of neighbor list

    :type  idx: int

    :return: an instance of :class:`NeighList` wrapping access to neighbor list data

    :rtype:  NeighList

    """

    if idx < 0:

        return None

    return NeighList(self, idx)

  def find_pair_neighlist(self, style, exact=True, nsub=0, request=0):

    """Find neighbor list index of pair style neighbor list

    Try finding pair instance that matches style. If exact is set, the pair must

    match style exactly. If exact is 0, style must only be contained. If pair is

    of style pair/hybrid, style is instead matched the nsub-th hybrid sub-style.

    Once the pair instance has been identified, multiple neighbor list requests

    may be found. Every neighbor list is uniquely identified by its request

    index. Thus, providing this request index ensures that the correct neighbor

    list index is returned.

    :param style: name of pair style that should be searched for

    :type  style: string

    :param exact: controls whether style should match exactly or only must be contained in pair style name, defaults to True

    :type  exact: bool, optional

    :param nsub:  match nsub-th hybrid sub-style, defaults to 0

    :type  nsub:  int, optional

    :param request:   index of neighbor list request, in case there are more than one, defaults to 0

    :type  request:   int, optional

    :return: neighbor list index if found, otherwise -1

    :rtype:  int

     """

    style = style.encode()

    exact = int(exact)

    idx = self.lib.lammps_find_pair_neighlist(self.lmp, style, exact, nsub, request)

    return self.get_neighlist(idx)

  def find_fix_neighlist(self, fixid, request=0):

    """Find neighbor list index of fix neighbor list

    :param fixid: name of fix

    :type  fixid: string

    :param request:   index of neighbor list request, in case there are more than one, defaults to 0

    :type  request:   int, optional

    :return: neighbor list index if found, otherwise -1

    :rtype:  int

     """

    fixid = fixid.encode()

    idx = self.lib.lammps_find_fix_neighlist(self.lmp, fixid, request)

    return self.get_neighlist(idx)

  def find_compute_neighlist(self, computeid, request=0):

    """Find neighbor list index of compute neighbor list

    :param computeid: name of compute

    :type  computeid: string

    :param request:   index of neighbor list request, in case there are more than one, defaults to 0

    :type  request:   int, optional

    :return: neighbor list index if found, otherwise -1

    :rtype:  int

     """

    computeid = computeid.encode()

    idx = self.lib.lammps_find_compute_neighlist(self.lmp, computeid, request)

    return self.get_neighlist(idx)

  def get_neighlist_size(self, idx):

    """Return the number of elements in neighbor list with the given index

    :param idx: neighbor list index

    :type  idx: int

    :return: number of elements in neighbor list with index idx

    :rtype:  int

     """

    return self.lib.lammps_neighlist_num_elements(self.lmp, idx)

  def get_neighlist_element_neighbors(self, idx, element):

    """Return data of neighbor list entry

    :param element: neighbor list index

    :type  element: int

    :param element: neighbor list element index

    :type  element: int

    :return: tuple with atom local index, number of neighbors and array of neighbor local atom indices

    :rtype:  (int, int, numpy.array)

    """

    c_iatom = c_int()

    c_numneigh = c_int()

    c_neighbors = POINTER(c_int)()

    self.lib.lammps_neighlist_element_neighbors(self.lmp, idx, element, byref(c_iatom), byref(c_numneigh), byref(c_neighbors))

    neighbors = self.numpy.iarray(c_int, c_neighbors, c_numneigh.value, 1)

    return c_iatom.value, c_numneigh.value, neighbors

# -------------------------------------------------------------------------

# -------------------------------------------------------------------------

# -------------------------------------------------------------------------