diff --git a/src/Makefile b/src/Makefile
index 07bf2443bb91620b7dd01006a6641ca81a0a5132..5cf76f305c6645daf9c828d839f85637f9b6b00f 100755
--- a/src/Makefile
+++ b/src/Makefile
@@ -43,14 +43,14 @@ endif
 
 PACKAGE = asphere body class2 colloid coreshell dipole fld gpu granular kim \
 	  kokkos kspace manybody mc meam misc molecule mpiio opt peri poems \
-	  qeq reax replica rigid shock snap srd voronoi xtc
+	  python qeq reax replica rigid shock snap srd voronoi xtc
 
 PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \
 	   user-cuda user-eff user-fep user-intel user-lb user-misc \
 	   user-molfile user-omp user-phonon user-qmmm user-quip \
 	   user-reaxc user-sph
 
-PACKLIB = gpu kim kokkos meam poems reax voronoi \
+PACKLIB = gpu kim kokkos meam poems python reax voronoi \
 	  user-atc user-awpmd user-colvars user-cuda user-molfile \
 	  user-qmmm user-quip 
 
diff --git a/src/PYTHON/Install.sh b/src/PYTHON/Install.sh
new file mode 100755
index 0000000000000000000000000000000000000000..02c7acc81276126ad59c79fc7afe5499e5b69729
--- /dev/null
+++ b/src/PYTHON/Install.sh
@@ -0,0 +1,69 @@
+# Install/unInstall package files in LAMMPS
+# mode = 0/1/2 for uninstall/install/update
+
+mode=$1
+
+# arg1 = file, arg2 = file it depends on
+
+action () {
+  if (test $mode = 0) then
+    rm -f ../$1
+  elif (! cmp -s $1 ../$1) then
+    if (test -z "$2" || test -e ../$2) then
+      cp $1 ..
+      if (test $mode = 2) then
+        echo "  updating src/$1"
+      fi
+    fi
+  elif (test -n "$2") then
+    if (test ! -e ../$2) then
+      rm -f ../$1
+    fi
+  fi
+}
+
+# force rebuild of files with LMP_KOKKOS switch
+# also variable so its *.d dependence on changed python_wrapper.h is rebuilt
+
+touch ../python_wrapper.h
+touch ../variable.cpp
+
+# all package files with no dependencies
+
+for file in *.cpp *.h; do
+  action $file
+done
+
+# edit 2 Makefile.package files to include/exclude package info
+
+if (test $1 = 1) then
+
+  if (test -e ../Makefile.package) then
+    sed -i -e 's/[^ \t]*python[^ \t]* //' ../Makefile.package
+    sed -i -e 's/[^ \t]*PYTHON[^ \t]* //g' ../Makefile.package
+    sed -i -e 's|^PKG_INC =[ \t]*|&-DLMP_PYTHON |' ../Makefile.package
+    sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(python_SYSINC) |' ../Makefile.package
+    sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(python_SYSLIB) |' ../Makefile.package
+    sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(python_SYSPATH) |' ../Makefile.package
+  fi
+
+  if (test -e ../Makefile.package.settings) then
+    sed -i -e '/^include.*python.*$/d' ../Makefile.package.settings
+    # multiline form needed for BSD sed on Macs
+    sed -i -e '4 i \
+include ..\/..\/lib\/python\/Makefile.lammps
+' ../Makefile.package.settings
+  fi
+
+elif (test $1 = 0) then
+
+  if (test -e ../Makefile.package) then
+    sed -i -e 's/[^ \t]*python[^ \t]* //' ../Makefile.package
+    sed -i -e 's/[^ \t]*PYTHON[^ \t]* //g' ../Makefile.package
+  fi
+
+  if (test -e ../Makefile.package.settings) then
+    sed -i -e '/^include.*python.*$/d' ../Makefile.package.settings
+  fi
+
+fi
diff --git a/src/PYTHON/python.cpp b/src/PYTHON/python.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ba9373a1de87cf429eb27cf9c2eb0878a7ca57d3
--- /dev/null
+++ b/src/PYTHON/python.cpp
@@ -0,0 +1,400 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#include "Python.h"
+#include "python.h"
+#include "force.h"
+#include "input.h"
+#include "variable.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{NONE,INT,DOUBLE,STRING,PTR};
+
+#define VALUELENGTH 64               // also in variable.cpp
+
+/* ---------------------------------------------------------------------- */
+
+Python::Python(LAMMPS *lmp) : Pointers(lmp)
+{
+  python_exists = 1;
+
+  pyMain = NULL;
+
+  // pfuncs stores interface info for each Python function
+
+  nfunc = 0;
+  pfuncs = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+Python::~Python()
+{
+  // clean up
+
+  for (int i = 0; i < nfunc; i++) {
+    delete [] pfuncs[i].name;
+    deallocate(i);
+    PyObject *pFunc = (PyObject *) pfuncs[i].pFunc;
+    Py_XDECREF(pFunc);
+  }
+
+  // shutdown Python interpreter
+
+  if (pyMain) Py_Finalize();
+
+  memory->sfree(pfuncs);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Python::command(int narg, char **arg)
+{
+  if (narg < 2) error->all(FLERR,"Invalid python command");
+
+  // if invoke is only keyword, invoke the previously defined function
+
+  if (narg == 2 && strcmp(arg[1],"invoke") == 0) {
+    int ifunc = find(arg[0]);
+    if (ifunc < 0) error->all(FLERR,"Python invoke of undefined function");
+
+    char *str = NULL;
+    if (noutput) {
+      str = input->variable->pythonstyle(pfuncs[ifunc].ovarname,
+                                         pfuncs[ifunc].name);
+      if (!str) 
+        error->all(FLERR,"Python variable does not match Python function");
+    }
+
+    invoke_function(ifunc,str);
+    return;
+  }
+
+  // parse optional args, invoke is not allowed in this mode
+
+  ninput = noutput = 0;
+  istr = NULL;
+  ostr = NULL;
+  format = NULL;
+  char *pyfile = NULL;
+  char *herestr = NULL;
+  int existflag = 0;
+
+  int iarg = 1;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"input") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Invalid python command");
+      ninput = force->inumeric(FLERR,arg[iarg+1]);
+      if (ninput < 0) error->all(FLERR,"Invalid python command");
+      iarg += 2;
+      istr = new char*[ninput];
+      if (iarg+ninput > narg) error->all(FLERR,"Invalid python command");
+      for (int i = 0; i < ninput; i++) istr[i] = arg[iarg+i];
+      iarg += ninput;
+    } else if (strcmp(arg[iarg],"return") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Invalid python command");
+      noutput = 1;
+      ostr = arg[iarg+1];
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"format") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Invalid python command");
+      int n = strlen(arg[iarg+1]) + 1;
+      format = new char[n];
+      strcpy(format,arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"file") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Invalid python command");
+      int n = strlen(arg[iarg+1]) + 1;
+      pyfile = new char[n];
+      strcpy(pyfile,arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"here") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Invalid python command");
+      herestr = arg[iarg+1];
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"exists") == 0) {
+      existflag = 1;
+      iarg++;
+    } else error->all(FLERR,"Invalid python command");
+  }
+
+  if (pyfile && herestr) error->all(FLERR,"Invalid python command");
+  if (pyfile && existflag) error->all(FLERR,"Invalid python command");
+  if (herestr && existflag) error->all(FLERR,"Invalid python command");
+
+  // create or overwrite entry in pfuncs vector with name = arg[0]
+
+  int ifunc = create_entry(arg[0]);
+
+  // one-time intitialization of Python interpreter
+  // Py_SetArgv() enables finding of *.py module files in current dir
+  //   only needed for module load, not for direct file read into __main__
+  // pymain stores pointer to main module
+
+  if (pyMain == NULL) {
+    if (Py_IsInitialized())
+      error->all(FLERR,"Cannot embed Python when also "
+                 "extending Python with LAMMPS");
+    Py_Initialize();
+
+    //char *arg = (char *) "./lmp";
+    //PySys_SetArgv(1,&arg);
+
+    //PyObject *pName = PyString_FromString("__main__");
+    //if (!pName) error->all(FLERR,"Bad pName");
+    //PyObject *pModule = PyImport_Import(pName);
+    //Py_DECREF(pName);
+
+    PyObject *pModule = PyImport_AddModule("__main__");
+    if (!pModule) error->all(FLERR,"Could not initialize embedded Python");
+    //if (!pModule) error->one(FLERR,"Could not initialize embedded Python");
+    pyMain = (void *) pModule;
+  }
+
+  // send Python code to Python interpreter
+  // file: read the file via PyRun_SimpleFile()
+  // here: process the here string directly
+  // exist: do nothing, assume code has already been run
+
+  if (pyfile) {
+    FILE *fp = fopen(pyfile,"r");
+    if (fp == NULL) error->all(FLERR,"Could not open Python file");
+    int err = PyRun_SimpleFile(fp,pyfile);
+    if (err) error->all(FLERR,"Could not process Python file");
+  } else if (herestr) {
+    int err = PyRun_SimpleString(herestr);
+    if (err) error->all(FLERR,"Could not process Python string");
+  }
+
+  // pFunc = function object for requested function
+
+  PyObject *pModule = (PyObject *) pyMain;
+  PyObject *pFunc = PyObject_GetAttrString(pModule,pfuncs[ifunc].name);
+  if (!pFunc) error->all(FLERR,"Could not find Python function");
+  if (!PyCallable_Check(pFunc)) 
+    error->all(FLERR,"Python function is not callable");
+  pfuncs[ifunc].pFunc = (void *) pFunc;
+
+  // clean-up input storage
+
+  delete [] istr;
+  delete [] format;
+  delete [] pyfile;
+}
+
+/* ------------------------------------------------------------------ */
+
+void Python::invoke_function(int ifunc, char *result)
+{
+  PyObject *pValue;
+  char *str;
+
+  PyObject *pFunc = (PyObject *) pfuncs[ifunc].pFunc;
+
+  // create Python tuple of input arguments
+
+  int ninput = pfuncs[ifunc].ninput;
+  PyObject *pArgs = PyTuple_New(ninput);
+  if (!pArgs) error->all(FLERR,"Could not create Python function arguments");
+
+  for (int i = 0; i < ninput; i++) {
+    int itype = pfuncs[ifunc].itype[i];
+    if (itype == INT) {
+      if (pfuncs[ifunc].ivarflag[i]) {
+        str = input->variable->retrieve(pfuncs[ifunc].svalue[i]);
+        if (!str) 
+          error->all(FLERR,"Could not evaluate Python function input variable");
+        pValue = PyInt_FromLong(atoi(str));
+      } else pValue = PyInt_FromLong(pfuncs[ifunc].ivalue[i]);
+    } else if (itype == DOUBLE) {
+      if (pfuncs[ifunc].ivarflag[i]) {
+        str = input->variable->retrieve(pfuncs[ifunc].svalue[i]);
+        if (!str) 
+          error->all(FLERR,"Could not evaluate Python function input variable");
+        pValue = PyFloat_FromDouble(atof(str));
+      } else pValue = PyFloat_FromDouble(pfuncs[ifunc].dvalue[i]);
+    } else if (itype == STRING) {
+      if (pfuncs[ifunc].ivarflag[i]) {
+        str = input->variable->retrieve(pfuncs[ifunc].svalue[i]);
+        if (!str) 
+          error->all(FLERR,"Could not evaluate Python function input variable");
+        pValue = PyString_FromString(str);
+      } else pValue = PyString_FromString(pfuncs[ifunc].svalue[i]);
+    } else if (itype == PTR) {
+      pValue = PyCObject_FromVoidPtr((void *) lmp,NULL);
+    }
+    PyTuple_SetItem(pArgs,i,pValue); 
+  }
+
+  // call the Python function
+  // error check with one() since only some procs may fail
+
+  pValue = PyObject_CallObject(pFunc,pArgs);
+  if (!pValue) error->one(FLERR,"Python function evaluation failed");
+  Py_DECREF(pArgs);
+
+  // function returned a value
+  // assign it to result string stored by python-style variable
+
+  if (pfuncs[ifunc].noutput) {
+    int otype = pfuncs[ifunc].otype;
+    if (otype == INT) {
+      sprintf(result,"%ld",PyInt_AsLong(pValue));
+    } else if (otype == DOUBLE) {
+      sprintf(result,"%.15g",PyFloat_AsDouble(pValue));
+    } else if (otype == STRING) {
+      char *pystr = PyString_AsString(pValue);
+      strncpy(result,pystr,VALUELENGTH-1);
+    }
+    Py_DECREF(pValue);
+  }
+}
+
+/* ------------------------------------------------------------------ */
+
+int Python::find(char *name)
+{
+  for (int i = 0; i < nfunc; i++)
+    if (strcmp(name,pfuncs[i].name) == 0) return i;
+  return -1;
+}
+
+/* ------------------------------------------------------------------ */
+
+int Python::variable_match(char *name, char *varname, int numeric)
+{
+  int ifunc = find(name);
+  if (ifunc < 0) return -1;
+  if (pfuncs[ifunc].noutput == 0) return -1;
+  if (strcmp(pfuncs[ifunc].ovarname,varname) != 0) return -1;
+  if (numeric && pfuncs[ifunc].otype == STRING) return -1;
+  return ifunc;
+}
+
+/* ------------------------------------------------------------------ */
+
+int Python::create_entry(char *name)
+{
+  // ifunc = index to entry by name in pfuncs vector, can be old or new
+  // free old vectors if overwriting old pfunc
+
+  int ifunc = find(name);
+
+  if (ifunc < 0) {
+    ifunc = nfunc;
+    nfunc++;
+    pfuncs = (PyFunc *) 
+      memory->srealloc(pfuncs,nfunc*sizeof(struct PyFunc),"python:pfuncs");
+    int n = strlen(name) + 1;
+    pfuncs[ifunc].name = new char[n];
+    strcpy(pfuncs[ifunc].name,name);
+  } else deallocate(ifunc);
+
+  pfuncs[ifunc].ninput = ninput;
+  pfuncs[ifunc].noutput = noutput;
+  
+  if (!format && ninput+noutput)
+    error->all(FLERR,"Invalid python command");
+  else if (format && strlen(format) != ninput+noutput) 
+    error->all(FLERR,"Invalid python command");
+
+  // process inputs as values or variables
+
+  pfuncs[ifunc].itype = new int[ninput];
+  pfuncs[ifunc].ivarflag = new int[ninput];
+  pfuncs[ifunc].ivalue = new int[ninput];
+  pfuncs[ifunc].dvalue = new double[ninput];
+  pfuncs[ifunc].svalue = new char*[ninput];
+
+  for (int i = 0; i < ninput; i++) {
+    pfuncs[ifunc].svalue[i] = NULL;
+    char type = format[i];
+    if (type == 'i') {
+      pfuncs[ifunc].itype[i] = INT;
+      if (strstr(istr[i],"v_") == istr[i]) {
+        pfuncs[ifunc].ivarflag[i] = 1;
+        int n = strlen(&istr[i][2]) + 1;
+        pfuncs[ifunc].svalue[i] = new char[n];
+        strcpy(pfuncs[ifunc].svalue[i],&istr[i][2]);
+      } else {
+        pfuncs[ifunc].ivarflag[i] = 0;
+        pfuncs[ifunc].ivalue[i] = force->inumeric(FLERR,istr[i]);
+      }
+    } else if (type == 'f') {
+      pfuncs[ifunc].itype[i] = DOUBLE;
+      if (strstr(istr[i],"v_") == istr[i]) {
+        pfuncs[ifunc].ivarflag[i] = 1;
+        int n = strlen(&istr[i][2]) + 1;
+        pfuncs[ifunc].svalue[i] = new char[n];
+        strcpy(pfuncs[ifunc].svalue[i],&istr[i][2]);
+      } else {
+        pfuncs[ifunc].ivarflag[i] = 0;
+        pfuncs[ifunc].dvalue[i] = force->numeric(FLERR,istr[i]);
+      }
+    } else if (type == 's') {
+      pfuncs[ifunc].itype[i] = STRING;
+      if (strstr(istr[i],"v_") == istr[i]) {
+        pfuncs[ifunc].ivarflag[i] = 1;
+        int n = strlen(&istr[i][2]) + 1;
+        pfuncs[ifunc].svalue[i] = new char[n];
+        strcpy(pfuncs[ifunc].svalue[i],&istr[i][2]);
+      } else {
+        pfuncs[ifunc].ivarflag[i] = 0;
+        int n = strlen(istr[i]) + 1;
+        pfuncs[ifunc].svalue[i] = new char[n];
+        strcpy(pfuncs[ifunc].svalue[i],istr[i]);
+      }
+    } else if (type == 'p') {
+      pfuncs[ifunc].ivarflag[i] = 0;
+      pfuncs[ifunc].itype[i] = PTR;
+      if (strcmp(istr[i],"SELF") != 0) 
+        error->all(FLERR,"Invalid python command");
+
+    } else error->all(FLERR,"Invalid python command");
+  }
+  
+  // process output as value or variable
+
+  pfuncs[ifunc].ovarname = NULL;
+  if (!noutput) return ifunc;
+
+  char type = format[ninput];
+  if (type == 'i') pfuncs[ifunc].otype = INT;
+  else if (type == 'f') pfuncs[ifunc].otype = DOUBLE;
+  else if (type == 's') pfuncs[ifunc].otype = STRING;
+  else error->all(FLERR,"Invalid python command");
+
+  if (strstr(ostr,"v_") != ostr) error->all(FLERR,"Invalid python command");
+  int n = strlen(&ostr[2]) + 1;
+  pfuncs[ifunc].ovarname = new char[n];
+  strcpy(pfuncs[ifunc].ovarname,&ostr[2]);
+
+  return ifunc;
+}
+
+/* ------------------------------------------------------------------ */
+
+void Python::deallocate(int i)
+{
+  delete [] pfuncs[i].itype;
+  delete [] pfuncs[i].ivarflag;
+  delete [] pfuncs[i].ivalue;
+  delete [] pfuncs[i].dvalue;
+  for (int j = 0; j < pfuncs[i].ninput; j++)
+    delete [] pfuncs[i].svalue[j];
+  delete [] pfuncs[i].svalue;
+  delete [] pfuncs[i].ovarname;
+}
diff --git a/src/PYTHON/python.h b/src/PYTHON/python.h
new file mode 100644
index 0000000000000000000000000000000000000000..2d30cb8880e0286bfaa9191875603f192171521b
--- /dev/null
+++ b/src/PYTHON/python.h
@@ -0,0 +1,59 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+
+#ifndef LMP_PYTHON_H
+#define LMP_PYTHON_H
+
+#include "pointers.h"
+
+namespace LAMMPS_NS {
+
+class Python : protected Pointers {
+ public:
+  int python_exists;
+
+  Python(class LAMMPS *);
+  ~Python();
+  void command(int, char **);
+  void invoke_function(int, char *);
+  int find(char *);
+  int variable_match(char *, char *, int);
+
+ private:
+  int ninput,noutput;
+  char **istr;
+  char *ostr,*format;
+  void *pyMain;
+  
+  struct PyFunc {
+    char *name;
+    int ninput,noutput;
+    int *itype,*ivarflag;
+    int *ivalue;
+    double *dvalue;
+    char **svalue;
+    int otype;
+    char *ovarname;
+    void *pFunc;
+  };
+
+  PyFunc *pfuncs;
+  int nfunc;
+
+  int create_entry(char *);
+  void deallocate(int);
+};
+
+}
+
+#endif
diff --git a/src/input.cpp b/src/input.cpp
index 4b6e27880f2ca694b94bcb49c373a0555feed974..455b8dc8696427b466a067a2c4afe74f618392a8 100644
--- a/src/input.cpp
+++ b/src/input.cpp
@@ -160,17 +160,33 @@ void Input::file()
       m = 0;
       while (1) {
         if (maxline-m < 2) reallocate(line,maxline,0);
+
+	// end of file reached, so break
+	// n == 0 if nothing read, else n = line with str terminator
+	
         if (fgets(&line[m],maxline-m,infile) == NULL) {
           if (m) n = strlen(line) + 1;
           else n = 0;
           break;
         }
+
+	// continue if last char read was not a newline
+	// could happen if line is very long
+
         m = strlen(line);
         if (line[m-1] != '\n') continue;
-        
+
+	// continue reading if final printable char is & char
+	// or if odd number of triple quotes
+	// else break with n = line with str terminator
+
         m--;
         while (m >= 0 && isspace(line[m])) m--;
         if (m < 0 || line[m] != '&') {
+	  if (numtriple(line) % 2) {
+	    m += 2;
+	    continue;
+	  }
           line[m+1] = '\0';
           n = m+2;
           break;
@@ -301,11 +317,11 @@ char *Input::one(const char *single)
 /* ----------------------------------------------------------------------
    parse copy of command line by inserting string terminators
    strip comment = all chars from # on
-   replace all $ via variable substitution
+   replace all $ via variable substitution except within quotes
    command = first word
    narg = # of args
    arg[] = individual args
-   treat text between single/double quotes as one arg
+   treat text between single/double/triple quotes as one arg via nextword()
 ------------------------------------------------------------------------- */
 
 void Input::parse()
@@ -317,17 +333,32 @@ void Input::parse()
   strcpy(copy,line);
   
   // strip any # comment by replacing it with 0
-  // do not strip # inside single/double quotes
-  
-  char quote = '\0';
+  // do not strip from a # inside single/double/triple quotes
+  // quoteflag = 1,2,3 when encounter first single/double,triple quote
+  // quoteflag = 0 when encounter matching single/double,triple quote
+
+  int quoteflag = 0;
   char *ptr = copy;
   while (*ptr) {
-    if (*ptr == '#' && !quote) {
+    if (*ptr == '#' && !quoteflag) {
       *ptr = '\0';
       break;
     }
-    if (*ptr == quote) quote = '\0';
-    else if (*ptr == '"' || *ptr == '\'') quote = *ptr;
+    if (quoteflag == 0) {
+      if (strstr(ptr,"\"\"\"") == ptr) {
+	quoteflag = 3;
+	ptr += 2;
+      }
+      else if (*ptr == '"') quoteflag = 2;
+      else if (*ptr == '\'') quoteflag = 1;
+    } else {
+      if (quoteflag == 3 && strstr(ptr,"\"\"\"") == ptr) {
+	quoteflag = 0;
+	ptr += 2;
+      }
+      else if (quoteflag == 2 && *ptr == '"') quoteflag = 0;
+      else if (quoteflag == 1 && *ptr == '\'') quoteflag = 0;
+    }
     ptr++;
   }
   
@@ -344,7 +375,7 @@ void Input::parse()
   
   // point arg[] at each subsequent arg in copy string
   // nextword() inserts string terminators into copy string to delimit args
-  // nextword() treats text between single/double quotes as one arg
+  // nextword() treats text between single/double/triple quotes as one arg
   
   narg = 0;
   ptr = next;
@@ -364,31 +395,53 @@ void Input::parse()
    find next word in str
    insert 0 at end of word
    ignore leading whitespace
-   treat text between single/double quotes as one arg
+   treat text between single/double/triple quotes as one arg
    matching quote must be followed by whitespace char if not end of string
    strip quotes from returned word
-   return ptr to start of word
-   return next = ptr after word or NULL if word ended with 0
-   return NULL if no word in string
+   return ptr to start of word or NULL if no word in string
+   also return next = ptr after word
 ------------------------------------------------------------------------- */
 
 char *Input::nextword(char *str, char **next)
 {
   char *start,*stop;
   
+  // start = first non-whitespace char
+
   start = &str[strspn(str," \t\n\v\f\r")];
   if (*start == '\0') return NULL;
   
-  if (*start == '"' || *start == '\'') {
-    stop = strchr(&start[1],*start);
+  // if start is single/double/triple quote:
+  //   start = first char beyond quote
+  //   stop = first char of matching quote
+  //   next = first char beyond matching quote
+  //   next must be NULL or whitespace
+  // if start is not single/double/triple quote:
+  //   stop = first whitespace char after start
+  //   next = char after stop, or stop itself if stop is NULL
+
+  if (strstr(start,"\"\"\"") == start) {
+    stop = strstr(&start[3],"\"\"\"");
     if (!stop) error->all(FLERR,"Unbalanced quotes in input line");
-    if (stop[1] && !isspace(stop[1]))
+    start += 3;
+    *next = stop+3;
+    if (**next && !isspace(**next))
       error->all(FLERR,"Input line quote not followed by whitespace");
+  } else if (*start == '"' || *start == '\'') {
+    stop = strchr(&start[1],*start);
+    if (!stop) error->all(FLERR,"Unbalanced quotes in input line");
     start++;
-  } else stop = &start[strcspn(start," \t\n\v\f\r")];
-  
-  if (*stop == '\0') *next = NULL;
-  else *next = stop+1;
+    *next = stop+1;
+    if (**next && !isspace(**next))
+      error->all(FLERR,"Input line quote not followed by whitespace");
+  } else {
+    stop = &start[strcspn(start," \t\n\v\f\r")];
+    if (*stop == '\0') *next = stop;
+    else *next = stop+1;
+  }
+
+  // set stop to NULL to terminate word
+
   *stop = '\0';
   return start;
 }
@@ -404,7 +457,7 @@ void Input::substitute(char *&str, char *&str2, int &max, int &max2, int flag)
 {
   // use str2 as scratch space to expand str, then copy back to str
   // reallocate str and str2 as necessary
-  // do not replace $ inside single/double quotes
+  // do not replace $ inside single/double/triple quotes
   // var = pts at variable name, ended by NULL
   //   if $ is followed by '{', trailing '}' becomes NULL
   //   else $x becomes x followed by NULL
@@ -413,7 +466,7 @@ void Input::substitute(char *&str, char *&str2, int &max, int &max2, int flag)
   int i,n,paren_count;
   char immediate[256];
   char *var,*value,*beyond;
-  char quote = '\0';
+  int quoteflag = 0;
   char *ptr = str;
   
   n = strlen(str) + 1;
@@ -422,10 +475,11 @@ void Input::substitute(char *&str, char *&str2, int &max, int &max2, int flag)
   char *ptr2 = str2;
   
   while (*ptr) {
+
     // variable substitution
     
-    if (*ptr == '$' && !quote) {
-      
+    if (*ptr == '$' && !quoteflag) {
+
       // value = ptr to expanded variable
       // variable name between curly braces, e.g. ${a}
       
@@ -440,7 +494,7 @@ void Input::substitute(char *&str, char *&str2, int &max, int &max2, int flag)
         beyond = ptr + strlen(var) + 3;
         value = variable->retrieve(var);
         
-        // immediate variable between parenthesis, e.g. $(1/2)
+      // immediate variable between parenthesis, e.g. $(1/2)
         
       } else if (*(ptr+1) == '(') {
         var = ptr+2;
@@ -492,10 +546,29 @@ void Input::substitute(char *&str, char *&str2, int &max, int &max2, int flag)
       
       continue;
     }
-    
-    if (*ptr == quote) quote = '\0';
-    else if (*ptr == '"' || *ptr == '\'') quote = *ptr;
-    
+
+    // quoteflag = 1,2,3 when encounter first single/double,triple quote
+    // quoteflag = 0 when encounter matching single/double,triple quote
+    // copy 2 extra triple quote chars into str2
+
+    if (quoteflag == 0) {
+      if (strstr(ptr,"\"\"\"") == ptr) {
+	quoteflag = 3;
+	*ptr2++ = *ptr++;
+	*ptr2++ = *ptr++;
+      } 
+      else if (*ptr == '"') quoteflag = 2;
+      else if (*ptr == '\'') quoteflag = 1;
+    } else {
+      if (quoteflag == 3 && strstr(ptr,"\"\"\"") == ptr) {
+	quoteflag = 0;
+	*ptr2++ = *ptr++;
+	*ptr2++ = *ptr++;
+      }
+      else if (quoteflag == 2 && *ptr == '"') quoteflag = 0;
+      else if (quoteflag == 1 && *ptr == '\'') quoteflag = 0;
+    }
+
     // copy current character into str2
     
     *ptr2++ = *ptr++;
@@ -509,6 +582,21 @@ void Input::substitute(char *&str, char *&str2, int &max, int &max2, int flag)
   strcpy(str,str2);
 }
 
+/* ----------------------------------------------------------------------
+   return number of triple quotes in line
+------------------------------------------------------------------------- */
+
+int Input::numtriple(char *line)
+{
+  int count = 0;
+  char *ptr = line;
+  while ((ptr = strstr(ptr,"\"\"\""))) {
+    ptr += 3;
+    count++;
+  }
+  return count;
+}
+
 /* ----------------------------------------------------------------------
    rellocate a string
    if n > 0: set max >= n in increments of DELTALINE
@@ -543,6 +631,7 @@ int Input::execute_command()
   else if (!strcmp(command,"next")) next_command();
   else if (!strcmp(command,"partition")) partition();
   else if (!strcmp(command,"print")) print();
+  else if (!strcmp(command,"python")) python();
   else if (!strcmp(command,"quit")) quit();
   else if (!strcmp(command,"shell")) shell();
   else if (!strcmp(command,"variable")) variable_command();
@@ -978,6 +1067,13 @@ void Input::print()
 
 /* ---------------------------------------------------------------------- */
 
+void Input::python()
+{
+  variable->python_command(narg,arg);
+}
+
+/* ---------------------------------------------------------------------- */
+
 void Input::quit()
 {
   if (narg) error->all(FLERR,"Illegal quit command");
diff --git a/src/input.h b/src/input.h
index a6df1a4a8c30305d6ec94a27ea303e18b738dddd..d1c3d6c1adf1f19e90983cd45b90ccf15a05fb83 100644
--- a/src/input.h
+++ b/src/input.h
@@ -58,6 +58,7 @@ class Input : protected Pointers {
 
   void parse();                          // parse an input text line
   char *nextword(char *, char **);       // find next word in string with quotes
+  int numtriple(char *);                 // count number of triple quotes
   void reallocate(char *&, int &, int);  // reallocate a char string
   int execute_command();                 // execute a single command
 
@@ -71,6 +72,7 @@ class Input : protected Pointers {
   void next_command();
   void partition();
   void print();
+  void python();
   void quit();
   void shell();
   void variable_command();
diff --git a/src/library.cpp b/src/library.cpp
index 8fa32444d034d79cab0d369a5cde2a0fa195c521..545bef4e9f5f0497191266576afdb30118659001 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -357,6 +357,19 @@ void *lammps_extract_variable(void *ptr, char *name, char *group)
   return NULL;
 }
 
+/* ----------------------------------------------------------------------
+   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 = lmp->input->variable->set_string(name,str);
+  return err;
+}
+
 /* ----------------------------------------------------------------------
    return the total number of atoms in the system
    useful before call to lammps_get_atoms() so can pre-allocate vector
@@ -365,6 +378,7 @@ void *lammps_extract_variable(void *ptr, char *name, char *group)
 int lammps_get_natoms(void *ptr)
 {
   LAMMPS *lmp = (LAMMPS *) ptr;
+
   if (lmp->atom->natoms > MAXSMALLINT) return 0;
   int natoms = static_cast<int> (lmp->atom->natoms);
   return natoms;
diff --git a/src/library.h b/src/library.h
index fc85c1280900bdb6c489296a0df902e064988de8..fddc7312a4c1d9c76587701c3eb46f1c351fb579 100644
--- a/src/library.h
+++ b/src/library.h
@@ -37,6 +37,7 @@ void *lammps_extract_compute(void *, char *, int, int);
 void *lammps_extract_fix(void *, char *, int, int, int, int);
 void *lammps_extract_variable(void *, char *, char *);
 
+int lammps_set_variable(void *, char *, char *);
 int lammps_get_natoms(void *);
 void lammps_gather_atoms(void *, char *, int, int, void *);
 void lammps_scatter_atoms(void *, char *, int, int, void *);
diff --git a/src/python_wrapper.h b/src/python_wrapper.h
new file mode 100644
index 0000000000000000000000000000000000000000..2975d5192f3fd854bc18871683a50f7bd4e16dfe
--- /dev/null
+++ b/src/python_wrapper.h
@@ -0,0 +1,47 @@
+/* -*- 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.
+------------------------------------------------------------------------- */
+
+#ifndef LMP_PYTHON_WRAPPER_H
+#define LMP_PYTHON_WRAPPER_H
+
+// true interface to embedded Python
+// used when PYTHON package is installed
+
+#ifdef LMP_PYTHON
+
+#include "python.h"
+
+#else
+
+// dummy interface to PYTHON
+// needed for compiling when PYTHON is not installed
+
+namespace LAMMPS_NS {
+
+class Python {
+ public:
+  int python_exists;
+
+  Python(class LAMMPS *) {python_exists = 0;}
+  ~Python() {}
+  void command(int, char **) {}
+  void invoke_function(int, char *) {}
+  int find(char *) {return -1;}
+  int variable_match(char *, char *, int) {return -1;}
+
+};
+
+}
+
+#endif
+#endif
diff --git a/src/variable.cpp b/src/variable.cpp
index 79c8ad65e20b1ee16be7584891c63b115ab175a0..f9828a4a12f14648ff248c7d27abb512ce6100b8 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -28,14 +28,15 @@
 #include "compute.h"
 #include "fix.h"
 #include "fix_store.h"
+#include "force.h"
 #include "output.h"
 #include "thermo.h"
 #include "random_mars.h"
 #include "math_const.h"
 #include "atom_masks.h"
+#include "python_wrapper.h"
 #include "memory.h"
 #include "error.h"
-#include "force.h"
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
@@ -44,13 +45,13 @@ using namespace MathConst;
 #define MAXLEVEL 4
 #define MAXLINE 256
 #define CHUNK 1024
-#define VALUELENGTH 64
+#define VALUELENGTH 64               // also in python.cpp
 #define MAXFUNCARG 6
 
 #define MYROUND(a) (( a-floor(a) ) >= .5) ? ceil(a) : floor(a)
 
 enum{INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,GETENV,
-     SCALARFILE,ATOMFILE,FORMAT,EQUAL,ATOM};
+     SCALARFILE,ATOMFILE,FORMAT,EQUAL,ATOM,PYTHON};
 enum{ARG,OP};
 
 // customize by adding a function
@@ -106,6 +107,10 @@ Variable::Variable(LAMMPS *lmp) : Pointers(lmp)
   precedence[MULTIPLY] = precedence[DIVIDE] = precedence[MODULO] = 6;
   precedence[CARAT] = 7;
   precedence[UNARY] = precedence[NOT] = 8;
+
+  // Python wrapper, real or dummy
+
+  python = new Python(lmp);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -131,6 +136,8 @@ Variable::~Variable()
 
   delete randomequal;
   delete randomatom;
+
+  delete python;
 }
 
 /* ----------------------------------------------------------------------
@@ -413,9 +420,36 @@ void Variable::set(int narg, char **arg)
       copy(1,&arg[2],data[nvar]);
     }
 
+  // PYTHON
+  // replace pre-existing var if also style PYTHON (allows it to be reset)
+  // num = 2, which = 1st value
+  // data = 2 values, 1st is Python func to invoke, 2nd is filled by invoke
+
+  } else if (strcmp(arg[1],"python") == 0) {
+    if (narg != 3) error->all(FLERR,"Illegal variable command");
+    if (!python->python_exists)
+      error->all(FLERR,"LAMMPS is not built with Python embedded");
+    int ivar = find(arg[0]);
+    if (ivar >= 0) {
+      if (style[ivar] != PYTHON)
+        error->all(FLERR,"Cannot redefine variable as a different style");
+      delete [] data[ivar][0];
+      copy(1,&arg[2],data[ivar]);
+      replaceflag = 1;
+    } else {
+      if (nvar == maxvar) grow();
+      style[nvar] = PYTHON;
+      num[nvar] = 2;
+      which[nvar] = 1;
+      pad[nvar] = 0;
+      data[nvar] = new char*[num[nvar]];
+      copy(1,&arg[2],data[nvar]);
+      data[nvar][1] = new char[VALUELENGTH];
+    }
+
   } else error->all(FLERR,"Illegal variable command");
 
-  // set name of variable, if not replacing (STRING/EQUAL/ATOM)
+  // set name of variable, if not replacing (EQUAL/ATOM/STRING/PYTHON)
   // name must be all alphanumeric chars or underscores
 
   if (replaceflag) return;
@@ -446,6 +480,23 @@ void Variable::set(char *name, int narg, char **arg)
   delete [] newarg;
 }
 
+/* ----------------------------------------------------------------------
+   set existing STRING variable to str
+   return 0 if successful
+   return -1 if variable doesn't exist or isn't a STRING variable
+   called via library interface, so external programs can set variables
+------------------------------------------------------------------------- */
+
+int Variable::set_string(char *name, char *str)
+{
+  int ivar = find(name);
+  if (ivar < 0) return -1;
+  if (style[ivar] != STRING) return -1;
+  delete [] data[ivar][0];
+  copy(1,&str,data[ivar]);
+  return 0;
+}
+
 /* ----------------------------------------------------------------------
    increment variable(s)
    return 0 if OK if successfully incremented
@@ -470,11 +521,12 @@ int Variable::next(int narg, char **arg)
       error->all(FLERR,"All variables in next command must be same style");
   }
 
-  // invalid styles STRING or EQUAL or WORLD or ATOM or GETENV or FORMAT
+  // invalid styles: STRING, EQUAL, WORLD, ATOM, GETENV, FORMAT, PYTHON
 
   int istyle = style[find(arg[0])];
-  if (istyle == STRING || istyle == EQUAL || istyle == WORLD
-      || istyle == GETENV || istyle == ATOM || istyle == FORMAT)
+  if (istyle == STRING || istyle == EQUAL || istyle == WORLD || 
+      istyle == GETENV || istyle == ATOM || istyle == FORMAT || 
+      istyle == PYTHON)
     error->all(FLERR,"Invalid variable style with next command");
 
   // if istyle = UNIVERSE or ULOOP, insure all such variables are incremented
@@ -587,16 +639,96 @@ int Variable::next(int narg, char **arg)
   return flag;
 }
 
+/* ----------------------------------------------------------------------
+   search for name in list of variables names
+   return index or -1 if not found
+------------------------------------------------------------------------- */
+
+int Variable::find(char *name)
+{
+  for (int i = 0; i < nvar; i++)
+    if (strcmp(name,names[i]) == 0) return i;
+  return -1;
+}
+
+/* ----------------------------------------------------------------------
+   initialize one atom's storage values in all VarReaders via fix STORE
+   called when atom is created
+------------------------------------------------------------------------- */
+
+void Variable::set_arrays(int i)
+{
+  for (int i = 0; i < nvar; i++)
+    if (reader[i] && style[i] == ATOMFILE)
+      reader[i]->fixstore->vstore[i] = 0.0;
+}
+
+/* ----------------------------------------------------------------------
+   called by python command in input script
+   simply pass input script line args to Python class
+------------------------------------------------------------------------- */
+
+void Variable::python_command(int narg, char **arg)
+{
+  if (!python->python_exists)
+    error->all(FLERR,"LAMMPS is not built with Python embedded");
+  python->command(narg,arg);
+}
+
+/* ----------------------------------------------------------------------
+   return 1 if variable is EQUAL or PYTHON numeric style, 0 if not
+   this is checked before call to compute_equal() to return a double
+------------------------------------------------------------------------- */
+
+int Variable::equalstyle(int ivar)
+{
+  if (style[ivar] == EQUAL) return 1;
+  if (style[ivar] == PYTHON) {
+    int ifunc = python->variable_match(data[ivar][0],names[ivar],1);
+    if (ifunc < 0) return 0;
+    else return 1;
+  }
+  return 0;
+}
+
+/* ----------------------------------------------------------------------
+   return 1 if variable is ATOM or ATOMFILE style, 0 if not
+   this is checked before call to compute_atom() to return a vector of doubles
+------------------------------------------------------------------------- */
+
+int Variable::atomstyle(int ivar)
+{
+  if (style[ivar] == ATOM || style[ivar] == ATOMFILE) return 1;
+  return 0;
+}
+
+/* ----------------------------------------------------------------------
+   check if variable with name is PYTHON and matches funcname
+   called by Python class before it invokes a Python function
+   return data storage so Python function can return a value for this variable
+   return NULL if not a match
+------------------------------------------------------------------------- */
+
+char *Variable::pythonstyle(char *name, char *funcname)
+{
+  int ivar = find(name);
+  if (ivar == -1) return NULL;
+  if (style[ivar] != PYTHON) return NULL;
+  if (strcmp(data[ivar][0],funcname) != 0) return NULL;
+  return data[ivar][1];
+}
+
 /* ----------------------------------------------------------------------
    return ptr to the data text associated with a variable
-   if INDEX or WORLD or UNIVERSE or STRING or SCALARFILE var, 
+   if INDEX or WORLD or UNIVERSE or STRING or SCALARFILE, 
      return ptr to stored string
-   if LOOP or ULOOP var, write int to data[0] and return ptr to string
-   if EQUAL var, evaluate variable and put result in str
-   if FORMAT var, evaluate its variable and put formatted result in str
-   if GETENV var, query environment and put result in str
-   if ATOM or ATOMFILE var, return NULL
-   return NULL if no variable with name or which value is bad,
+   if LOOP or ULOOP, write int to data[0] and return ptr to string
+   if EQUAL, evaluate variable and put result in str
+   if FORMAT, evaluate its variable and put formatted result in str
+   if GETENV, query environment and put result in str
+   if PYTHON, evaluate Python function, it will put result in str
+   if ATOM or ATOMFILE, return NULL
+   return NULL if no variable with name, or which value is bad,
      caller must respond
 ------------------------------------------------------------------------- */
 
@@ -606,6 +738,10 @@ char *Variable::retrieve(char *name)
   if (ivar == -1) return NULL;
   if (which[ivar] >= num[ivar]) return NULL;
 
+  if (eval_in_progress[ivar]) 
+    error->all(FLERR,"Variable has circular dependency");
+  eval_in_progress[ivar] = 1;
+
   char *str;
   if (style[ivar] == INDEX || style[ivar] == WORLD ||
       style[ivar] == UNIVERSE || style[ivar] == STRING || 
@@ -625,16 +761,14 @@ char *Variable::retrieve(char *name)
     strcpy(data[ivar][0],result);
     str = data[ivar][0];
   } else if (style[ivar] == EQUAL) {
-    eval_in_progress[ivar] = 1;
     double answer = evaluate(data[ivar][0],NULL);
-    eval_in_progress[ivar] = 0;
     sprintf(data[ivar][1],"%.15g",answer);
     str = data[ivar][1];
   } else if (style[ivar] == FORMAT) {
     int jvar = find(data[ivar][0]);
     if (jvar == -1) return NULL;
     if (!equalstyle(jvar)) return NULL;
-    double answer = evaluate(data[jvar][0],NULL);
+    double answer = compute_equal(jvar);
     sprintf(data[ivar][2],data[ivar][1],answer);
     str = data[ivar][2];
   } else if (style[ivar] == GETENV) {
@@ -647,13 +781,24 @@ char *Variable::retrieve(char *name)
     }
     strcpy(data[ivar][1],result);
     str = data[ivar][1];
+  } else if (style[ivar] == PYTHON) {
+    int ifunc = python->variable_match(data[ivar][0],names[ivar],0);
+    if (ifunc < 0) 
+      error->all(FLERR,"Python variable does not match Python function");
+    python->invoke_function(ifunc,data[ivar][1]);
+    str = data[ivar][1];
   } else if (style[ivar] == ATOM || style[ivar] == ATOMFILE) return NULL;
 
+  eval_in_progress[ivar] = 0;
+
   return str;
 }
 
 /* ----------------------------------------------------------------------
    return result of equal-style variable evaluation
+   can be EQUAL style or PYTHON numeric style
+   for PYTHON, don't need to check python->variable_match() error return,
+     since caller will have already checked via equalstyle()
 ------------------------------------------------------------------------- */
 
 double Variable::compute_equal(int ivar)
@@ -662,7 +807,14 @@ double Variable::compute_equal(int ivar)
     error->all(FLERR,"Variable has circular dependency");
   eval_in_progress[ivar] = 1;
 
-  double value = evaluate(data[ivar][0],NULL);
+  double value;
+  if (style[ivar] == EQUAL) value = evaluate(data[ivar][0],NULL);
+  else if (style[ivar] == PYTHON) {
+    int ifunc = python->find(data[ivar][0]);
+    if (ifunc < 0) error->all(FLERR,"Python variable has no function");
+    python->invoke_function(ifunc,data[ivar][1]);
+    value = atof(data[ivar][1]);
+  }
 
   eval_in_progress[ivar] = 0;
   return value;
@@ -671,6 +823,7 @@ double Variable::compute_equal(int ivar)
 /* ----------------------------------------------------------------------
    return result of immediate equal-style variable evaluation
    called from Input::substitute()
+   don't need to flag eval_in_progress since is an immediate variable
 ------------------------------------------------------------------------- */
 
 double Variable::compute_equal(char *str)
@@ -744,54 +897,11 @@ void Variable::compute_atom(int ivar, int igroup,
   eval_in_progress[ivar] = 0;
 }
 
-/* ----------------------------------------------------------------------
-   search for name in list of variables names
-   return index or -1 if not found
-------------------------------------------------------------------------- */
-
-int Variable::find(char *name)
-{
-  for (int i = 0; i < nvar; i++)
-    if (strcmp(name,names[i]) == 0) return i;
-  return -1;
-}
-
-/* ----------------------------------------------------------------------
-   initialize one atom's storage values in all VarReaders via fix STORE
-   called when atom is created
-------------------------------------------------------------------------- */
-
-void Variable::set_arrays(int i)
-{
-  for (int i = 0; i < nvar; i++)
-    if (reader[i] && style[i] == ATOMFILE)
-      reader[i]->fixstore->vstore[i] = 0.0;
-}
-
-/* ----------------------------------------------------------------------
-   return 1 if variable is EQUAL style, 0 if not
-------------------------------------------------------------------------- */
-
-int Variable::equalstyle(int ivar)
-{
-  if (style[ivar] == EQUAL) return 1;
-  return 0;
-}
-
-/* ----------------------------------------------------------------------
-   return 1 if variable is ATOM or ATOMFILE style, 0 if not
-------------------------------------------------------------------------- */
-
-int Variable::atomstyle(int ivar)
-{
-  if (style[ivar] == ATOM || style[ivar] == ATOMFILE) return 1;
-  return 0;
-}
-
 /* ----------------------------------------------------------------------
    save copy of EQUAL style ivar formula in copy
    allocate copy here, later equal_restore() call will free it
    insure data[ivar][0] is of VALUELENGTH since will be overridden
+   next 3 functions are used by create_atoms to temporarily override variables
 ------------------------------------------------------------------------- */
 
 void Variable::equal_save(int ivar, char *&copy)
@@ -2712,8 +2822,6 @@ tagint Variable::int_between_brackets(char *&ptr, int varallow)
     int ivar = find(id);
     if (ivar < 0)
       error->all(FLERR,"Invalid variable name in variable formula");
-    if (eval_in_progress[ivar])
-      error->all(FLERR,"Variable has circular dependency");
 
     char *var = retrieve(id);
     if (var == NULL)
diff --git a/src/variable.h b/src/variable.h
index d5f4e0f57d12efd6995cf0db7475400b568587ff..703071da65753aa3ca34f30b70e2dda762deb3e6 100644
--- a/src/variable.h
+++ b/src/variable.h
@@ -25,11 +25,17 @@ class Variable : protected Pointers {
   ~Variable();
   void set(int, char **);
   void set(char *, int, char **);
+  int set_string(char *, char *);
   int next(int, char **);
+
   int find(char *);
   void set_arrays(int);
+  void python_command(int, char **);
+
   int equalstyle(int);
   int atomstyle(int);
+  char *pythonstyle(char *, char *);
+
   char *retrieve(char *);
   double compute_equal(int);
   double compute_equal(char *);
@@ -64,6 +70,8 @@ class Variable : protected Pointers {
   int precedence[17];      // precedence level of math operators
                            // set length to include up to OR in enum
 
+  class Python *python;    // ptr to embedded Python interpreter
+
   struct Tree {            // parse tree for atom-style variables
     double value;          // single scalar  
     double *array;         // per-atom or per-type list of doubles
@@ -78,6 +86,7 @@ class Variable : protected Pointers {
     Tree **extra;          // ptrs further down tree for nextra args
   };
 
+  int compute_python(int);
   void remove(int);
   void grow();
   void copy(int, char **, char **);