Merge pull request #742 from rbberger/fix_python_move

Fix python/move

doc/src/Section_commands.txt

@@ -620,7 +620,8 @@ USER-INTEL, k = KOKKOS, o = USER-OMP, t = OPT.
 "press/berendsen"_fix_press_berendsen.html,
 "print"_fix_print.html,
 "property/atom"_fix_property_atom.html,
-"python"_fix_python.html,
+"python/invoke"_fix_python_invoke.html,
+"python/move"_fix_python_move.html,
 "qeq/comb (o)"_fix_qeq_comb.html,
 "qeq/dynamic"_fix_qeq.html,
 "qeq/fire"_fix_qeq.html,

doc/src/Section_python.txt

@@ -123,7 +123,7 @@ code directly from an input script:
 
 "python"_python.html
 "variable python"_variable.html
-"fix python"_fix_python.html
+"fix python/invoke"_fix_python_invoke.html
 "pair_style python"_pair_python.html :ul
 
 The "python"_python.html command which can be used to define and
@@ -165,7 +165,7 @@ doc page for its python-style variables for more info, including
 examples of Python code you can write for both pure Python operations
 and callbacks to LAMMPS.
 
-The "fix python"_fix_python.html command can execute
+The "fix python/invoke"_fix_python_invoke.html command can execute
 Python code at selected timesteps during a simulation run.
 
 The "pair_style python"_pair_python command allows you to define

doc/src/fix_python_invoke.txt

@@ -6,14 +6,14 @@
 
 :line
 
-fix python command :h3
+fix python/invoke command :h3
 
 [Syntax:]
 
-fix ID group-ID python N callback function_name :pre
+fix ID group-ID python/invoke N callback function_name :pre
 
 ID, group-ID are ignored by this fix :ulb,l
-python = style name of this fix command :l
+python/invoke = style name of this fix command :l
 N = execute every N steps :l
 callback = {post_force} or {end_of_step} :l
   {post_force} = callback after force computations on atoms every N time steps
@@ -36,8 +36,8 @@ def end_of_step_callback(lammps_ptr):
     # access LAMMPS state using Python interface
 """ :pre
 
-fix pf  all python 50 post_force post_force_callback
-fix eos all python 50 end_of_step end_of_step_callback :pre
+fix pf  all python/invoke 50 post_force post_force_callback
+fix eos all python/invoke 50 end_of_step end_of_step_callback :pre
 
 [Description:]
 

doc/src/fix_python_move.txt

@@ -0,0 +1,102 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+fix python/move command :h3
+
+[Syntax:]
+
+fix python/move pymodule.CLASS :pre
+
+pymodule.CLASS = use class [CLASS] in module/file [pymodule] to compute how to move atoms
+
+[Examples:]
+
+fix  1 all python/move py_nve.NVE
+fix  1 all python/move py_nve.NVE_OPT :pre
+
+[Description:]
+
+The {python/move} fix style provides a way to define ways how particles
+are moved during an MD run from python script code, that is loaded from
+a file into LAMMPS and executed at the various steps where other fixes
+can be executed. This python script must contain specific python class
+definitions.
+
+This allows to implement complex position updates and also modified
+time integration methods. Due to python being an interpreted language,
+however, the performance of this fix can be moderately to significantly
+slower than the corresponding C++ code. For specific cases, this
+performance penalty can be limited through effective use of NumPy.
+
+:line
+
+The python module file has to start with the following code:
+
+from __future__ import print_function
+import lammps
+import ctypes
+import traceback
+import numpy as np
+#
+class LAMMPSFix(object):
+    def __init__(self, ptr, group_name="all"):
+        self.lmp = lammps.lammps(ptr=ptr)
+        self.group_name = group_name
+#
+class LAMMPSFixMove(LAMMPSFix):
+    def __init__(self, ptr, group_name="all"):
+        super(LAMMPSFixMove, self).__init__(ptr, group_name)
+#
+    def init(self):
+        pass
+#
+    def initial_integrate(self, vflag):
+        pass
+#
+    def final_integrate(self):
+        pass
+#
+    def initial_integrate_respa(self, vflag, ilevel, iloop):
+        pass
+#
+    def final_integrate_respa(self, ilevel, iloop):
+        pass
+#
+    def reset_dt(self):
+        pass :pre
+
+Any classes implementing new atom motion functionality have to be
+derived from the [LAMMPSFixMove] class, overriding the available
+methods as needed.
+
+Examples for how to do this are in the {examples/python} folder.
+
+:line
+
+[Restart, fix_modify, output, run start/stop, minimize info:]
+
+No information about this fix is written to "binary restart
+files"_restart.html.  None of the "fix_modify"_fix_modify.html options
+are relevant to this fix.  No global or per-atom quantities are stored
+by this fix for access by various "output
+commands"_Section_howto.html#howto_15.  No parameter of this fix can
+be used with the {start/stop} keywords of the "run"_run.html command.
+This fix is not invoked during "energy minimization"_minimize.html.
+
+[Restrictions:]
+
+This pair style is part of the PYTHON package.  It is only enabled if
+LAMMPS was built with that package.  See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info.
+
+[Related commands:]
+
+"fix nve"_fix_nve.html, "fix python/invoke"_fix_python_invoke.html
+
+[Default:] none
+

doc/src/fixes.txt

@@ -115,7 +115,8 @@ Fixes :h1
    fix_press_berendsen
    fix_print
    fix_property_atom
-   fix_python
+   fix_python_invoke
+   fix_python_move
    fix_qbmsst
    fix_qeq
    fix_qeq_comb

doc/src/lammps.book

@@ -244,7 +244,8 @@ fix_pour.html
 fix_press_berendsen.html
 fix_print.html
 fix_property_atom.html
-fix_python.html
+fix_python_invoke.html
+fix_python_move.html
 fix_qbmsst.html
 fix_qeq.html
 fix_qeq_comb.html

doc/src/python.txt

@@ -406,7 +406,7 @@ cases, LAMMPS has no simple way to check that something illogical is
 being attempted.
 
 The same applies to Python functions called during a simulation run at
-each time step using "fix python"_fix_python.html.
+each time step using "fix python/invoke"_fix_python_invoke.html.
 
 :line
 
@@ -493,6 +493,6 @@ different source files, problems may occur.
 
 [Related commands:]
 
-"shell"_shell.html, "variable"_variable.html, "fix python"_fix_python.html
+"shell"_shell.html, "variable"_variable.html, "fix python/invoke"_fix_python_invoke.html
 
 [Default:] none

examples/python/README.fix_python_move

@@ -0,0 +1,16 @@
+This folder contains several LAMMPS input scripts and a python module
+file py_nve.py to demonstrate the use of the fix style python/move
+to reimplement NVE using Python.
+
+in.fix_python_move_nve_melt:
+This is a version of the melt example which replaces the default NVE integrator
+with a Python implementation. Fix python/move is used to create an
+instance of the py_nve.NVE class which implements the required interface.
+It demonstrates how to access LAMMPS data as numpy arrays. This gives direct
+access to memory owned by the C++ code, allows easy manipulation through numpy
+operations and avoids unnecessary copies.
+
+in.fix_python_move_nve_melt_opt:
+This version of melt example uses NVE_Opt instead of NVE. While this Python
+implementation is still much slower than the native version, it shows that
+simple code transformations can lead to speedups.

examples/python/in.fix_python_invoke

@@ -33,8 +33,8 @@ def post_force_callback(lmp, v):
 """
 
 fix		1 all nve
-fix     2 all python 50 end_of_step end_of_step_callback
-fix     3 all python 50 post_force post_force_callback
+fix     2 all python/invoke 50 end_of_step end_of_step_callback
+fix     3 all python/invoke 50 post_force post_force_callback
 
 #dump		id all atom 50 dump.melt
 

examples/python/in.fix_python_move_nve_melt

@@ -0,0 +1,23 @@
+# 3d Lennard-Jones melt
+
+units		lj
+atom_style	atomic
+
+lattice		fcc 0.8442
+region		box block 0 10 0 10 0 10
+create_box	1 box
+create_atoms	1 box
+mass		* 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.3 bin
+neigh_modify	every 20 delay 0 check no
+
+fix		1 all python/move py_nve.NVE
+
+thermo		50
+run		250

examples/python/in.fix_python_move_nve_melt_opt

@@ -0,0 +1,23 @@
+# 3d Lennard-Jones melt
+
+units		lj
+atom_style	atomic
+
+lattice		fcc 0.8442
+region		box block 0 10 0 10 0 10
+create_box	1 box
+create_atoms	1 box
+mass		* 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.3 bin
+neigh_modify	every 20 delay 0 check no
+
+fix		1 all python/move py_nve.NVE_Opt
+
+thermo		50
+run		250

examples/python/py_nve.py

@@ -0,0 +1,118 @@
+from __future__ import print_function
+import lammps
+import ctypes
+import traceback
+import numpy as np
+
+class LAMMPSFix(object):
+    def __init__(self, ptr, group_name="all"):
+        self.lmp = lammps.lammps(ptr=ptr)
+        self.group_name = group_name
+
+class LAMMPSFixMove(LAMMPSFix):
+    def __init__(self, ptr, group_name="all"):
+        super(LAMMPSFixMove, self).__init__(ptr, group_name)
+
+    def init(self):
+        pass
+
+    def initial_integrate(self, vflag):
+        pass
+
+    def final_integrate(self):
+        pass
+
+    def initial_integrate_respa(self, vflag, ilevel, iloop):
+        pass
+
+    def final_integrate_respa(self, ilevel, iloop):
+        pass
+
+    def reset_dt(self):
+        pass
+
+
+class NVE(LAMMPSFixMove):
+    """ Python implementation of fix/nve """
+    def __init__(self, ptr, group_name="all"):
+        super(NVE, self).__init__(ptr)
+        assert(self.group_name == "all")
+
+    def init(self):
+        dt = self.lmp.extract_global("dt", 1)
+        ftm2v = self.lmp.extract_global("ftm2v", 1)
+        self.ntypes = self.lmp.extract_global("ntypes", 0)
+        self.dtv = dt
+        self.dtf = 0.5 * dt * ftm2v
+
+    def initial_integrate(self, vflag):
+        nlocal = self.lmp.extract_global("nlocal", 0)
+        mass = self.lmp.numpy.extract_atom_darray("mass", self.ntypes+1)
+        atype = self.lmp.numpy.extract_atom_iarray("type", nlocal)
+        x = self.lmp.numpy.extract_atom_darray("x", nlocal, dim=3)
+        v = self.lmp.numpy.extract_atom_darray("v", nlocal, dim=3)
+        f = self.lmp.numpy.extract_atom_darray("f", nlocal, dim=3)
+
+        for i in range(x.shape[0]):
+            dtfm = self.dtf / mass[int(atype[i])]
+            v[i,:]+= dtfm * f[i,:]
+            x[i,:] += self.dtv * v[i,:]
+
+    def final_integrate(self):
+        nlocal = self.lmp.extract_global("nlocal", 0)
+        mass = self.lmp.numpy.extract_atom_darray("mass", self.ntypes+1)
+        atype = self.lmp.numpy.extract_atom_iarray("type", nlocal)
+        v = self.lmp.numpy.extract_atom_darray("v", nlocal, dim=3)
+        f = self.lmp.numpy.extract_atom_darray("f", nlocal, dim=3)
+
+        for i in range(v.shape[0]):
+            dtfm = self.dtf / mass[int(atype[i])]
+            v[i,:] += dtfm * f[i,:]
+
+
+class NVE_Opt(LAMMPSFixMove):
+    """ Performance-optimized Python implementation of fix/nve """
+    def __init__(self, ptr, group_name="all"):
+        super(NVE_Opt, self).__init__(ptr)
+        assert(self.group_name == "all")
+
+    def init(self):
+        dt = self.lmp.extract_global("dt", 1)
+        ftm2v = self.lmp.extract_global("ftm2v", 1)
+        self.ntypes = self.lmp.extract_global("ntypes", 0)
+        self.dtv = dt
+        self.dtf = 0.5 * dt * ftm2v
+        self.mass = self.lmp.numpy.extract_atom_darray("mass", self.ntypes+1)
+
+    def initial_integrate(self, vflag):        
+        nlocal = self.lmp.extract_global("nlocal", 0)
+        atype = self.lmp.numpy.extract_atom_iarray("type", nlocal)
+        x = self.lmp.numpy.extract_atom_darray("x", nlocal, dim=3)
+        v = self.lmp.numpy.extract_atom_darray("v", nlocal, dim=3)
+        f = self.lmp.numpy.extract_atom_darray("f", nlocal, dim=3)
+        dtf = self.dtf
+        dtv = self.dtv
+        mass = self.mass
+
+        dtfm = dtf / np.take(mass, atype)
+        dtfm.reshape((nlocal, 1))
+
+        for d in range(x.shape[1]):
+            v[:,d] += dtfm[:,0] * f[:,d]
+            x[:,d] += dtv * v[:,d]
+
+    def final_integrate(self):
+        nlocal = self.lmp.extract_global("nlocal", 0)
+        mass = self.lmp.numpy.extract_atom_darray("mass", self.ntypes+1)
+        atype = self.lmp.numpy.extract_atom_iarray("type", nlocal)
+        v = self.lmp.numpy.extract_atom_darray("v", nlocal, dim=3)
+        f = self.lmp.numpy.extract_atom_darray("f", nlocal, dim=3)
+        dtf = self.dtf
+        dtv = self.dtv
+        mass = self.mass
+
+        dtfm = dtf / np.take(mass, atype)
+        dtfm.reshape((nlocal, 1))
+
+        for d in range(v.shape[1]):
+            v[:,d] += dtfm[:,0] * f[:,d]

python/lammps.py

@@ -37,7 +37,7 @@ def get_ctypes_int(size):
     return c_int32
   elif size == 8:
     return c_int64
-  return c_int
+  return c_int 
 
 class MPIAbortException(Exception):
   def __init__(self, message):
@@ -266,25 +266,41 @@ class lammps(object):
         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))
+        def _ctype_to_numpy_int(self, ctype_int):
+          if ctype_int == c_int32:
+            return np.int32
+          elif ctype_int == c_int64:
+            return np.int64
+          return np.intc
 
-          a = np.frombuffer(ptr.contents, dtype=np.intc)
+        def extract_atom_iarray(self, name, nelem, dim=1):
+          if name in ['id', 'molecule']:
+            c_int_type = self.lmp.c_tagint
+          elif name in ['image']:
+            c_int_type = self.lmp.c_imageint
+          else:
+            c_int_type = c_int
+
+          np_int_type = self._ctype_to_numpy_int(c_int_type)
+
+          if dim == 1:
+            tmp = self.lmp.extract_atom(name, 0)
+            ptr = cast(tmp, POINTER(c_int_type * nelem))
+          else:
+            tmp = self.lmp.extract_atom(name, 1)
+            ptr = cast(tmp[0], POINTER(c_int_type * nelem * dim))
+
+          a = np.frombuffer(ptr.contents, dtype=np_int_type)
           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))
+            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))
+            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)

src/.gitignore

@@ -910,8 +910,10 @@
 /python_impl.cpp
 /python_impl.h
 /python_compat.h
-/fix_python.cpp
-/fix_python.h
+/fix_python_move.cpp
+/fix_python_move.h
+/fix_python_invoke.cpp
+/fix_python_invoke.h
 /pair_python.cpp
 /pair_python.h
 /reader_molfile.cpp

src/Makefile

@@ -17,12 +17,12 @@ SHLINK = liblammps.so
 OBJDIR =   Obj_$@
 OBJSHDIR = Obj_shared_$@
 
-SRC =	$(filter-out library.cpp,$(wildcard *.cpp))
-INC =	$(filter-out library.h,$(wildcard *.h))
+SRC =	$(wildcard *.cpp)
+INC =	$(wildcard *.h)
 OBJ = 	$(SRC:.cpp=.o)
 
-SRCLIB = $(filter-out main.cpp,$(SRC)) library.cpp
-OBJLIB = $(filter-out main.o,$(OBJ)) library.o
+SRCLIB = $(filter-out main.cpp,$(SRC))
+OBJLIB = $(filter-out main.o,$(OBJ))
 
 # Command-line options for mode: exe (default), shexe, lib, shlib
 
@@ -176,7 +176,7 @@ help:
 	  then cp Makefile.package.settings.empty Makefile.package.settings; fi
 	@cp Makefile.package Makefile.package.settings $(objdir)
 	@cd $(objdir); rm -f .depend; \
-	$(MAKE) $(MFLAGS) "SRC = $(SRC) library.cpp" "INC = $(INC) library.h" depend || :
+	$(MAKE) $(MFLAGS) "SRC = $(SRC)" "INC = $(INC)" depend || :
 ifeq ($(mode),exe)
 	@cd $(objdir); \
 	$(MAKE) $(MFLAGS) "OBJ = $(OBJ)" "INC = $(INC)" "SHFLAGS =" \

src/PYTHON/fix_python_invoke.cpp

@@ -18,7 +18,7 @@
 #include <Python.h>
 #include <stdio.h>
 #include <string.h>
-#include "fix_python.h"
+#include "fix_python_invoke.h"
 #include "atom.h"
 #include "force.h"
 #include "update.h"
@@ -32,13 +32,13 @@ using namespace FixConst;
 
 /* ---------------------------------------------------------------------- */
 
-FixPython::FixPython(LAMMPS *lmp, int narg, char **arg) :
+FixPythonInvoke::FixPythonInvoke(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
-  if (narg != 6) error->all(FLERR,"Illegal fix python command");
+  if (narg != 6) error->all(FLERR,"Illegal fix python/invoke command");
 
   nevery = force->inumeric(FLERR,arg[3]);
-  if (nevery <= 0) error->all(FLERR,"Illegal fix python command");
+  if (nevery <= 0) error->all(FLERR,"Illegal fix python/invoke command");
 
   // ensure Python interpreter is initialized
   python->init();
@@ -48,7 +48,7 @@ FixPython::FixPython(LAMMPS *lmp, int narg, char **arg) :
   } else if (strcmp(arg[4],"end_of_step") == 0) {
     selected_callback = END_OF_STEP;
   } else {
-    error->all(FLERR,"Unsupported callback name for fix/python");
+    error->all(FLERR,"Unsupported callback name for fix python/invoke");
   }
 
   // get Python function
@@ -74,14 +74,14 @@ FixPython::FixPython(LAMMPS *lmp, int narg, char **arg) :
 
 /* ---------------------------------------------------------------------- */
 
-int FixPython::setmask()
+int FixPythonInvoke::setmask()
 {
   return selected_callback;
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixPython::end_of_step()
+void FixPythonInvoke::end_of_step()
 {
   PyGILState_STATE gstate = PyGILState_Ensure();
 
@@ -96,7 +96,7 @@ void FixPython::end_of_step()
 
 /* ---------------------------------------------------------------------- */
 
-void FixPython::post_force(int vflag)
+void FixPythonInvoke::post_force(int vflag)
 {
   if (update->ntimestep % nevery != 0) return;
 

src/PYTHON/fix_python_invoke.h

@@ -13,21 +13,22 @@
 
 #ifdef FIX_CLASS
 
-FixStyle(python,FixPython)
+FixStyle(python,FixPythonInvoke)
+FixStyle(python/invoke,FixPythonInvoke)
 
 #else
 
-#ifndef LMP_FIX_PYTHON_H
-#define LMP_FIX_PYTHON_H
+#ifndef LMP_FIX_PYTHON_INVOKE_H
+#define LMP_FIX_PYTHON_INVOKE_H
 
 #include "fix.h"
 
 namespace LAMMPS_NS {
 
-class FixPython : public Fix {
+class FixPythonInvoke : public Fix {
  public:
-  FixPython(class LAMMPS *, int, char **);
-  virtual ~FixPython() {}
+  FixPythonInvoke(class LAMMPS *, int, char **);
+  virtual ~FixPythonInvoke() {}
   int setmask();
   virtual void end_of_step();
   virtual void post_force(int);

src/PYTHON/fix_python_move.cpp

@@ -0,0 +1,239 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Richard Berger (Temple U)
+------------------------------------------------------------------------- */
+
+#include <Python.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "fix_python_move.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "memory.h"
+#include "neigh_list.h"
+#include "python.h"
+#include "error.h"
+#include "python_compat.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixPythonMove::FixPythonMove(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  dynamic_group_allow = 1;
+  time_integrate = 1;
+
+  python->init();
+
+  py_move = NULL;
+
+  PyGILState_STATE gstate = PyGILState_Ensure();
+
+  // add current directory to PYTHONPATH
+  PyObject * py_path = PySys_GetObject("path");
+  PyList_Append(py_path, PY_STRING_FROM_STRING("."));
+
+
+  // create integrator instance
+  char * full_cls_name = arg[3];
+  char * lastpos = strrchr(full_cls_name, '.');
+
+  if (lastpos == NULL) {
+    error->all(FLERR,"Fix python/integrate requires fully qualified class name");
+  }
+
+  size_t module_name_length = strlen(full_cls_name) - strlen(lastpos);
+  size_t cls_name_length = strlen(lastpos)-1;
+
+  char * module_name = new char[module_name_length+1];
+  char * cls_name = new char[cls_name_length+1];
+  strncpy(module_name, full_cls_name, module_name_length);
+  module_name[module_name_length] = 0;
+
+  strcpy(cls_name, lastpos+1);
+
+  PyObject * pModule = PyImport_ImportModule(module_name);
+  if (!pModule) {
+    PyErr_Print();
+    PyErr_Clear();
+    PyGILState_Release(gstate);
+    error->all(FLERR,"Loading python integrator module failure");
+  }
+
+  // create LAMMPS atom type to potential file type mapping in python class
+  // by calling 'lammps_pair_style.map_coeff(name,type)'
+
+  PyObject *py_move_type = PyObject_GetAttrString(pModule, cls_name);
+  if (!py_move_type) {
+    PyErr_Print();
+    PyErr_Clear();
+    PyGILState_Release(gstate);
+    error->all(FLERR,"Could not find integrator class in module'");
+  }
+
+  delete [] module_name;
+  delete [] cls_name;
+
+  PyObject * ptr = PY_VOID_POINTER(lmp);
+  PyObject * arglist = Py_BuildValue("(O)", ptr);
+  PyObject * py_move_obj = PyObject_CallObject(py_move_type, arglist);
+  Py_DECREF(arglist);
+
+  if (!py_move_obj) {
+    PyErr_Print();
+    PyErr_Clear();
+    PyGILState_Release(gstate);
+    error->all(FLERR,"Could not instantiate instance of integrator class'");
+  }
+
+  // check object interface
+  py_move = (void *) py_move_obj;
+
+  PyGILState_Release(gstate);
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixPythonMove::~FixPythonMove()
+{
+  PyGILState_STATE gstate = PyGILState_Ensure();
+  if(py_move) Py_DECREF((PyObject*) py_move);
+  PyGILState_Release(gstate);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixPythonMove::setmask()
+{
+  int mask = 0;
+  mask |= INITIAL_INTEGRATE;
+  mask |= FINAL_INTEGRATE;
+  mask |= INITIAL_INTEGRATE_RESPA;
+  mask |= FINAL_INTEGRATE_RESPA;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPythonMove::init()
+{
+  PyGILState_STATE gstate = PyGILState_Ensure();
+  PyObject *py_move_obj = (PyObject *) py_move;
+  PyObject *py_init = PyObject_GetAttrString(py_move_obj,"init");
+  if (!py_init) {
+    PyErr_Print();
+    PyErr_Clear();
+    PyGILState_Release(gstate);
+    error->all(FLERR,"Could not find 'init' method'");
+  }
+  PyObject * result = PyEval_CallObject(py_init, NULL);
+  PyGILState_Release(gstate);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPythonMove::initial_integrate(int vflag)
+{
+  PyGILState_STATE gstate = PyGILState_Ensure();
+  PyObject *py_move_obj = (PyObject *) py_move;
+  PyObject *py_initial_integrate = PyObject_GetAttrString(py_move_obj,"initial_integrate");
+  if (!py_initial_integrate) {
+    PyErr_Print();
+    PyErr_Clear();
+    PyGILState_Release(gstate);
+    error->all(FLERR,"Could not find 'initial_integrate' method'");
+  }
+  PyObject * arglist = Py_BuildValue("(i)", vflag);
+  PyObject * result = PyEval_CallObject(py_initial_integrate, arglist);
+  Py_DECREF(arglist);
+  PyGILState_Release(gstate);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPythonMove::final_integrate()
+{
+  PyGILState_STATE gstate = PyGILState_Ensure();
+  PyObject *py_move_obj = (PyObject *) py_move;
+  PyObject *py_final_integrate = PyObject_GetAttrString(py_move_obj,"final_integrate");
+  if (!py_final_integrate) {
+    PyErr_Print();
+    PyErr_Clear();
+    PyGILState_Release(gstate);
+    error->all(FLERR,"Could not find 'final_integrate' method'");
+  }
+  PyObject * result = PyEval_CallObject(py_final_integrate, NULL);
+  PyGILState_Release(gstate);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPythonMove::initial_integrate_respa(int vflag, int ilevel, int iloop)
+{
+  PyGILState_STATE gstate = PyGILState_Ensure();
+  PyObject *py_move_obj = (PyObject *) py_move;
+  PyObject *py_initial_integrate_respa = PyObject_GetAttrString(py_move_obj,"initial_integrate_respa");
+  if (!py_initial_integrate_respa) {
+    PyErr_Print();
+    PyErr_Clear();
+    PyGILState_Release(gstate);
+    error->all(FLERR,"Could not find 'initial_integrate_respa' method'");
+  }
+  PyObject * arglist = Py_BuildValue("(iii)", vflag, ilevel, iloop);
+  PyObject * result = PyEval_CallObject(py_initial_integrate_respa, arglist);
+  Py_DECREF(arglist);
+  PyGILState_Release(gstate);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPythonMove::final_integrate_respa(int ilevel, int iloop)
+{
+  PyGILState_STATE gstate = PyGILState_Ensure();
+  PyObject *py_move_obj = (PyObject *) py_move;
+  PyObject *py_final_integrate_respa = PyObject_GetAttrString(py_move_obj,"final_integrate_respa");
+  if (!py_final_integrate_respa) {
+    PyErr_Print();
+    PyErr_Clear();
+    PyGILState_Release(gstate);
+    error->all(FLERR,"Could not find 'final_integrate_respa' method'");
+  }
+  PyObject * arglist = Py_BuildValue("(ii)", ilevel, iloop);
+  PyObject * result = PyEval_CallObject(py_final_integrate_respa, arglist);
+  Py_DECREF(arglist);
+  PyGILState_Release(gstate);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPythonMove::reset_dt()
+{
+  PyGILState_STATE gstate = PyGILState_Ensure();
+  PyObject *py_move_obj = (PyObject *) py_move;
+  PyObject *py_reset_dt = PyObject_GetAttrString(py_move_obj,"reset_dt");
+  if (!py_reset_dt) {
+    PyErr_Print();
+    PyErr_Clear();
+    PyGILState_Release(gstate);
+    error->all(FLERR,"Could not find 'reset_dt' method'");
+  }
+  PyObject * result = PyEval_CallObject(py_reset_dt, NULL);
+  PyGILState_Release(gstate);
+}

src/PYTHON/fix_python_move.h

@@ -0,0 +1,67 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   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.
+
+   Pair zero is a dummy pair interaction useful for requiring a
+   force cutoff distance in the absense of pair-interactions or
+   with hybrid/overlay if a larger force cutoff distance is required.
+
+   This can be used in conjunction with bond/create to create bonds
+   that are longer than the cutoff of a given force field, or to
+   calculate radial distribution functions for models without
+   pair interactions.
+
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(python/move,FixPythonMove)
+
+#else
+
+#ifndef LMP_FIX_PYTHON_MOVE_H
+#define LMP_FIX_PYTHON_MOVE_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixPythonMove : public Fix {
+ public:
+  FixPythonMove(LAMMPS *lmp, int narg, char **arg);
+  virtual ~FixPythonMove();
+
+  int setmask();
+  virtual void init();
+  virtual void initial_integrate(int);
+  virtual void final_integrate();
+  virtual void initial_integrate_respa(int, int, int);
+  virtual void final_integrate_respa(int, int);
+  virtual void reset_dt();
+
+ protected:
+  void * py_move;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+*/

src/PYTHON/pair_python.cpp

@@ -181,6 +181,7 @@ void PairPython::compute(int eflag, int vflag)
           error->all(FLERR,"Calling 'compute_force' function failed");
         }
         fpair = factor_lj*PyFloat_AsDouble(py_value);
+        Py_DECREF(py_value);
 
         f[i][0] += delx*fpair;
         f[i][1] += dely*fpair;
@@ -194,6 +195,7 @@ void PairPython::compute(int eflag, int vflag)
         if (eflag) {
           py_value = PyObject_CallObject(py_compute_energy,py_compute_args);
           evdwl = factor_lj*PyFloat_AsDouble(py_value);
+          Py_DECREF(py_value);
         } else evdwl = 0.0;
 
         if (evflag) ev_tally(i,j,nlocal,newton_pair,
@@ -467,6 +469,7 @@ double PairPython::single(int i, int j, int itype, int jtype, double rsq,
     error->all(FLERR,"Calling 'compute_force' function failed");
   }
   fforce = factor_lj*PyFloat_AsDouble(py_value);
+  Py_DECREF(py_value);
 
   py_value = PyObject_CallObject(py_compute_energy,py_compute_args);
   if (!py_value) {
@@ -476,6 +479,7 @@ double PairPython::single(int i, int j, int itype, int jtype, double rsq,
     error->all(FLERR,"Calling 'compute_energy' function failed");
   }
   double evdwl = factor_lj*PyFloat_AsDouble(py_value);
+  Py_DECREF(py_value);
 
   Py_DECREF(py_compute_args);
   PyGILState_Release(gstate);

src/PYTHON/python_impl.cpp

@@ -326,6 +326,7 @@ void PythonImpl::invoke_function(int ifunc, char *result)
   pValue = PyObject_CallObject(pFunc,pArgs);
 
   if (!pValue) {
+    PyErr_Print();
     PyGILState_Release(gstate);
     error->one(FLERR,"Python function evaluation failed");
   }

src/Purge.list

@@ -16,6 +16,9 @@ style_region.h
 style_neigh_bin.h
 style_neigh_pair.h
 style_neigh_stencil.h
+# deleted on 15 December 2017
+fix_python.cpp
+fix_python.h
 # deleted on 1 December 2017
 npair_half_bin_newtoff_intel.cpp
 npair_half_bin_newtoff_intel.h