simulation_hydro1_changed.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SELM via PyLAMMPs for Simulations\n",
    "Author: Paul Atzberger <br>\n",
    "http://atzberger.org/\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os;\n",
    "script_base_name = \"simulation_hydro1\";\n",
    "script_dir = os.getcwd();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the lammps module\n",
    "try:  \n",
    "  from selm_lammps.lammps import IPyLammps # use this for the pip install of pre-built package\n",
    "  lammps_import_comment = \"from selm_lammps.lammps import IPyLammps\";  \n",
    "  from selm_lammps import util as atz_util;\n",
    "except Exception as e:  \n",
    "  from lammps import IPyLammps # use this for direct install of package\n",
    "  lammps_import_comment = \"from lammps import IPyLammps\";\n",
    "  from atz_lammps import util as atz_util;\n",
    "except Exception as e: # if fails to import, report the exception   \n",
    "  print(e);\n",
    "  lammps_import_comment = \"import failed\";\n",
    "\n",
    "import numpy as np;\n",
    "import matplotlib;\n",
    "import matplotlib.pyplot as plt;\n",
    "import sys,shutil,pickle,pdb;\n",
    "\n",
    "import logging;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup SELM Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Setting up log file in filename = /home/xsu2/Downloads/hydro/output/simulation_hydro1/hydro_test001/main.log\n",
      "from selm_lammps.lammps import IPyLammps\n",
      "For running LAMMPS changing the current working directory to:\n",
      "/home/xsu2/Downloads/hydro/output/simulation_hydro1/hydro_test001\n"
     ]
    }
   ],
   "source": [
    "# @base_dir\n",
    "base_dir_output   = '%s/output/%s'%(script_dir,script_base_name);\n",
    "atz_util.create_dir(base_dir_output);\n",
    "\n",
    "dir_run_name = 'hydro';\n",
    "base_dir = '%s/%s_test001'%(base_dir_output,dir_run_name);\n",
    "\n",
    "# remove all data from dir\n",
    "atz_util.rm_dir(base_dir);\n",
    "\n",
    "# setup the directories\n",
    "base_dir_fig    = '%s/fig'%base_dir;\n",
    "atz_util.create_dir(base_dir_fig);\n",
    "\n",
    "base_dir_vtk    = '%s/vtk'%base_dir;\n",
    "atz_util.create_dir(base_dir_vtk);\n",
    "\n",
    "# setup logging\n",
    "atzLog = atz_util.AtzLogging(print,base_dir);\n",
    "print_log = atzLog.print_log;\n",
    "\n",
    "# print the import comment\n",
    "print_log(lammps_import_comment);\n",
    "\n",
    "# change directory for running LAMMPS in output\n",
    "print_log(\"For running LAMMPS changing the current working directory to:\\n%s\"%base_dir);\n",
    "os.chdir(base_dir); # base the current working directory\n",
    "#os.chdir(script_dir); # base the current working directory\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup LAMMPs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LAMMPS output is captured by PyLammps wrapper\n",
      "LAMMPS version: 29 Oct 2020 / 20201029\n",
      "MPI version: LAMMPS MPI STUBS for LAMMPS version 29 Oct 2020\n"
     ]
    }
   ],
   "source": [
    "L = IPyLammps();\n",
    "atz_util.print_version_info(L);    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Copy files to the output directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model files being copied:\n",
      "src = /home/xsu2/Downloads/hydro/Model1\n",
      "dst = /home/xsu2/Downloads/hydro/output/simulation_hydro1/hydro_test001/\n"
     ]
    }
   ],
   "source": [
    "# copy the model files to the destination\n",
    "src = script_dir + '/' + \"Model1\";\n",
    "dst = base_dir + '/';\n",
    "atz_util.copytree2(src,dst,symlinks=False,ignore=None);\n",
    "\n",
    "print_log(\"Model files being copied:\\n\" + \"src = \" + str(src) + \"\\n\" + \"dst = \" + str(dst));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Copying notebook to archive:\n",
      "src = /home/xsu2/Downloads/hydro/simulation_hydro1.ipynb\n",
      "dst = /home/xsu2/Downloads/hydro/output/simulation_hydro1/hydro_test001/archive__simulation_hydro1.ipynb\n"
     ]
    }
   ],
   "source": [
    "flag_copy_notebook_to_output = True;\n",
    "if flag_copy_notebook_to_output:\n",
    "  #cur_dir = os.getcwd();\n",
    "  #src = cur_dir + '/' + script_base_name + '.ipynb';\n",
    "  src = script_dir + '/' + script_base_name + '.ipynb';    \n",
    "  dst = base_dir + '/' + 'archive__' + script_base_name + '.ipynb';\n",
    "  shutil.copyfile(src, dst);\n",
    "  print_log(\"Copying notebook to archive:\\n\" + \"src = \" + str(src) + \"\\n\" + \"dst = \" + str(dst));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Common Physical Parameters (nano units)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reference units and parameters\n",
    "units = {'name':'nano','mu':1.0,'rho':0.001,\n",
    "         'KB':0.01380651,'T':298.15};\n",
    "units.update({'KBT':units['KB']*units['T']});"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up Lammps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup the Model and Simulation Files (such as .read_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "atom_name = structure_pts\n",
      "num_pts = 2\n",
      "atom_name = tracer_pts\n",
      "num_pts = 1331\n",
      "atom_name = forcing terms\n",
      "num_pts = 1331\n"
     ]
    }
   ],
   "source": [
    "num_dim = 3;\n",
    "box = np.zeros((num_dim,2));\n",
    "LL = 202.5; box[:,0] = -LL; box[:,1] = LL;\n",
    "\n",
    "# setup atoms\n",
    "I_id = 1; I_type = 1; atom_types = [];\n",
    "atom_list = []; atom_mass_list = []; atom_id_list = []; \n",
    "atom_mol_list = []; atom_name_list = [];\n",
    "\n",
    "# structure atoms\n",
    "atom_name = \"structure_pts\";\n",
    "atom_name_list.append(atom_name);\n",
    "atom_types.append(I_type); \n",
    "\n",
    "# number of atoms in observation\n",
    "num_pts = 2; \n",
    "x1 = np.array([-100,100]);\n",
    "x2 = np.array([0,0]);\n",
    "x3 = np.array([0,0]);\n",
    "x = np.stack((x1,x2,x3),axis=1); # shape = [num_pts,num_dim]\n",
    "num_pts = x.shape[0]; m0 = 1.123; \n",
    "atom_id = np.arange(I_id + 0,I_id + num_pts,dtype=int);\n",
    "mol_id = 1; atom_mol = np.ones(num_pts,dtype=int)*mol_id;\n",
    "atom_list.append(x); atom_mass_list.append(m0); \n",
    "atom_id_list.append(atom_id); atom_mol_list.append(atom_mol);\n",
    "I_type += 1; I_id += num_pts;\n",
    "print_log(\"atom_name = \" + str(atom_name));\n",
    "print_log(\"num_pts = \" + str(num_pts));\n",
    "\n",
    "#number of obstacle points 11*11=121\n",
    "#num_pts = 11; \n",
    "#x1_raw= np.linspace(-200,200,num_pts,endpoint=True)#x1 = np.array([-100,100]);\n",
    "#x1=np.tile(x1_raw,num_pts)\n",
    "#x2_raw= np.linspace(-200,200,num_pts,endpoint=True)\n",
    "#x2=np.tile(x1_raw,num_pts)#x2= np.linspace(-200,200,num_pts,endpoint=True)#x2 = np.array([0,0]);\n",
    "#x3=np.zeros(num_pts**2);#x3 = np.array([0,0]);\n",
    "#x = np.stack((x1,x2,x3),axis=1); # shape = [num_pts,num_dim]\n",
    "#num_pts = x.shape[0]; m0 = 1.123; \n",
    "#atom_id = np.arange(I_id + 0,I_id + num_pts,dtype=int);\n",
    "#mol_id = 2; atom_mol = np.ones(num_pts,dtype=int)*mol_id;\n",
    "#atom_list.append(x); atom_mass_list.append(m0); \n",
    "#atom_id_list.append(atom_id); atom_mol_list.append(atom_mol);\n",
    "#I_type += 1; I_id += num_pts;\n",
    "#print_log(\"atom_name = \" + str(atom_name));\n",
    "#print_log(\"num_pts = \" + str(num_pts));\n",
    "\n",
    "\n",
    "# tracer atoms\n",
    "flag_tracer = True;\n",
    "if flag_tracer:\n",
    "  atom_name = \"tracer_pts\";\n",
    "  atom_name_list.append(atom_name);\n",
    "  atom_types.append(I_type); \n",
    "  atom_types[I_type - 1] = I_type;  \n",
    "  num_pts_dir = 10; m0 = 1.123; \n",
    "  x1 = np.linspace(-LL,LL,num_pts_dir + 1,endpoint=False); dx = x1[1] - x1[0];\n",
    "  x1 = x1 + 0.5*dx;\n",
    "  xx = np.meshgrid(x1,x1,x1);\n",
    "  x = np.stack((xx[0].flatten(),xx[1].flatten(),xx[2].flatten()),axis=1); # shape = [num_pts,num_dim]\n",
    "  #ipdb.set_trace();\n",
    "  num_pts = x.shape[0];\n",
    "  atom_id = np.arange(I_id + 0,I_id + num_pts,dtype=int);\n",
    "  mol_id = 2; atom_mol = np.ones(x.shape[0],dtype=int)*mol_id;\n",
    "  atom_list.append(x); atom_mass_list.append(m0); \n",
    "  atom_id_list.append(atom_id); atom_mol_list.append(atom_mol);\n",
    "  I_type += 1; I_id += num_pts;\n",
    "  print_log(\"atom_name = \" + str(atom_name));\n",
    "  print_log(\"num_pts = \" + str(num_pts));\n",
    "\n",
    "# atoms serve as forcing terms\n",
    "flag_force = True;\n",
    "if flag_force:\n",
    "  atom_name = \"forcing terms\";\n",
    "  atom_name_list.append(atom_name);\n",
    "  atom_types.append(I_type); \n",
    "  atom_types[I_type - 1] = I_type;  \n",
    "  num_pts_dir = 10; m0 = 1.123; \n",
    "  x1 = np.linspace(-LL,LL,num_pts_dir + 1,endpoint=False); dx = x1[1] - x1[0];\n",
    "  x1 = x1 + 0.3*dx;\n",
    "  xx = np.meshgrid(x1,x1,x1);\n",
    "  x = np.stack((xx[0].flatten(),xx[1].flatten(),xx[2].flatten()),axis=1); # shape = [num_pts,num_dim]\n",
    "  #ipdb.set_trace();\n",
    "  num_pts = x.shape[0];\n",
    "  atom_id = np.arange(I_id + 0,I_id + num_pts,dtype=int);\n",
    "  mol_id = 3; atom_mol = np.ones(x.shape[0],dtype=int)*mol_id;\n",
    "  atom_list.append(x); atom_mass_list.append(m0); \n",
    "  atom_id_list.append(atom_id); atom_mol_list.append(atom_mol);\n",
    "  I_type += 1; I_id += num_pts;\n",
    "  print_log(\"atom_name = \" + str(atom_name));\n",
    "  print_log(\"num_pts = \" + str(num_pts));\n",
    "\n",
    "# summary data    \n",
    "# get total number of atoms\n",
    "atom_types = np.array(atom_types,dtype=int);\n",
    "num_atoms = I_id - 1; # total number of atoms\n",
    "\n",
    "# setup bonds\n",
    "I_id = 1; I_type = 1; bond_types = []; bond_name_list = [];\n",
    "bond_list = []; bond_coeff_list = []; bond_id_list = [];\n",
    "\n",
    "flag_bond_1 = False;\n",
    "if flag_bond_1:\n",
    "  bond_types.append(I_type);\n",
    "  #bond_name_list.append(\"fene_1\");\n",
    "  bond_name_list.append(\"harmonic_1\");\n",
    "  \n",
    "  KBT = units['KBT']; ell = 5.0; K = 0.5*KBT/(ell*ell); r0 = 20;\n",
    "  b = \"harmonic %.7f %7f\"%(K,r0);\n",
    "  print_log(\"bond:\");\n",
    "  print_log(\"KBT = \" + str(KBT)); print_log(\"r0 = \" + str(r0));  \n",
    "  print_log(\"bond coeff = \" + str(b));\n",
    "  bond_coeff_list.append(b);\n",
    "\n",
    "  I0 = atz_util.atz_find_name(atom_name_list,\"structure_pts\"); I_atom_type = atom_types[I0];\n",
    "  atom_id = atom_id_list[I_atom_type - 1]; nn = atom_id.shape[0];\n",
    "  bonds = np.zeros((nn,2),dtype=int); \n",
    "  bond_id = np.zeros(bonds.shape[0],dtype=int);\n",
    "  for i in range(0,nn):\n",
    "    i1 = atom_id[i % nn]; i2 = atom_id[(i + 1)%nn];  # base 1 indexing, closed loop\n",
    "    bonds[i,0] = i1; bonds[i,1] = i2;\n",
    "    bond_id[i] = I_id;\n",
    "    I_id += 1;  \n",
    "  bond_list.append(bonds); bond_id_list.append(bond_id);\n",
    "  I_type += 1;\n",
    "\n",
    "# summary data    \n",
    "num_bonds = I_id - 1;\n",
    "bond_types = np.array(bond_types,dtype=int);\n",
    "\n",
    "# setup angles\n",
    "I_id = 1; I_type =1 ; angle_types = []; angle_name_list = [];\n",
    "angle_list = []; angle_coeff_list = []; angle_id_list = [];\n",
    "\n",
    "flag_angles_1 = False;\n",
    "if flag_angles_1:\n",
    "  angle_name_list.append(\"atom_type_1\");\n",
    "  angle_types.append(I_type);\n",
    "  #KBT = 2478959.87; K = 10*KBT; theta_0 = 180.0; # degrees\n",
    "  KBT = units['KBT']; K = 5*KBT; theta_0 = 180.0; # degrees\n",
    "  b = \"harmonic %.7f %.7f\"%(K,theta_0);\n",
    "  angle_coeff_list.append(b);\n",
    "\n",
    "# build angle bonds for type 1 atoms with type 1 atoms, closed loop\n",
    "if flag_angles_1:\n",
    "  I0 = atz_util.atz_find_name(atom_name_list,\"polymer_pts\"); I_atom_type = atom_types[I0];\n",
    "  atom_id = atom_id_list[I_atom_type - 1]; nn = atom_id.shape[0];\n",
    "  angles = np.zeros((nn,3),dtype=int);\n",
    "  angle_id = np.zeros(angles.shape[0],dtype=int);\n",
    "  for i in range(0,nn):\n",
    "    i1 = atom_id[i]; i2 = atom_id[(i + 1)%nn]; i3 = atom_id[(i + 2)%nn]; # base 1 indexing\n",
    "    angles[i,0] = i1; angles[i,1] = i2; angles[i,2] = i3;\n",
    "    angle_id[i] = I_id; I_id += 1;\n",
    "  angle_list.append(angles); angle_id_list.append(angle_id);\n",
    "  I_type += 1;\n",
    "\n",
    "# summary data    \n",
    "num_angles = I_id - 1;\n",
    "angle_types = np.array(angle_types,dtype=int);\n",
    "\n",
    "# store the model information\n",
    "model_info = {};\n",
    "model_info.update({'num_dim':num_dim,'box':box,'atom_types':atom_types,\n",
    "          'atom_list':atom_list,'atom_mass_list':atom_mass_list,'atom_name_list':atom_name_list,\n",
    "          'atom_id_list':atom_id_list,'atom_mol_list':atom_mol_list,\n",
    "          'bond_types':bond_types,'bond_list':bond_list,'bond_id_list':bond_id_list,\n",
    "          'bond_coeff_list':bond_coeff_list,'bond_name_list':bond_name_list,\n",
    "          'angle_types':angle_types,'angle_list':angle_list,'angle_id_list':angle_id_list,\n",
    "          'angle_coeff_list':angle_coeff_list,'angle_name_list':angle_name_list});"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing model data .pickle\n",
      "filename = model_setup.pickle\n",
      "Writing model data .read_data\n",
      "filename = Model.LAMMPS_read_data\n"
     ]
    }
   ],
   "source": [
    "# write .pickle data with the model setup information\n",
    "filename = \"model_setup.pickle\";\n",
    "print_log(\"Writing model data .pickle\");\n",
    "print_log(\"filename = \" + filename);\n",
    "s = model_info;\n",
    "f = open(filename,'wb'); pickle.dump(s,f); f.close();\n",
    "\n",
    "# write the model .read_data file for lammps\n",
    "filename = \"Model.LAMMPS_read_data\";\n",
    "print_log(\"Writing model data .read_data\");\n",
    "print_log(\"filename = \" + filename);\n",
    "atz_util.write_read_data(filename=filename,print_log=print_log,**model_info);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!cat Model.LAMMPS_read_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform the simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sending commands to LAMMPs\n",
      "\n",
      "# =========================================================================\n",
      "# LAMMPS main parameter file and script                                    \n",
      "#                                                                          \n",
      "# Author: Paul J. Atzberger.               \n",
      "#\n",
      "# Based on script generated by SELM Model Builder.\n",
      "#                                                                          \n",
      "# =========================================================================\n",
      "\n",
      "# == Setup variables for the script \n",
      "\n",
      "variable dumpfreq         equal    1\n",
      "variable restart          equal    0\n",
      "variable neighborSkinDist equal    1.0 # distance for bins beyond force cut-off (1.0 = 1.0 Ang for units = real) \n",
      "variable baseFilename     universe Model\n",
      "\n",
      "# == Setup the log file\n",
      "#log         ${baseFilename}.LAMMPS_logFile\n",
      "\n",
      "#(New) setup the force term by adding a group of random atoms with specific force (Constant force for simple case)\n",
      "#create_atoms 1 random 100 127569 box\n",
      "#group obstacle type 1\n",
      "#fix 4 force_term setforce 3.0 0.0 0.0\n",
      "#(/New)\n",
      "\n",
      "# == Setup style of the run\n",
      "\n",
      "# type of units to use in the simulation (units used are in fact: amu, nm, ns, Kelvins)\n",
      "units       nano\n",
      "\n",
      "# indicates possible types allowed for interactions between the atoms\n",
      "atom_style  angle \n",
      "\n",
      "# indicates possible types allowed for bonds between the atoms \n",
      "# bond_style  hybrid harmonic\n",
      "\n",
      "# indicates possible types allowed for bond angles between the atoms \n",
      "angle_style none\n",
      "\n",
      "# indicates type of boundary conditions in each direction (p = periodic) \n",
      "boundary p p p \n",
      "read_data ${baseFilename}.LAMMPS_read_data # file of atomic coordinates and topology\n",
      "velocity all zero linear                   # initialize all atomic velocities initially to zero\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# We can send collection of commands using the triple quote notation\n",
    "s = \"\"\"\n",
    "# =========================================================================\n",
    "# LAMMPS main parameter file and script                                    \n",
    "#                                                                          \n",
    "# Author: Paul J. Atzberger.               \n",
    "#\n",
    "# Based on script generated by SELM Model Builder.\n",
    "#                                                                          \n",
    "# =========================================================================\n",
    "\n",
    "# == Setup variables for the script \n",
    "\n",
    "variable dumpfreq         equal    1\n",
    "variable restart          equal    0\n",
    "variable neighborSkinDist equal    1.0 # distance for bins beyond force cut-off (1.0 = 1.0 Ang for units = real) \n",
    "variable baseFilename     universe Model\n",
    "\n",
    "# == Setup the log file\n",
    "#log         ${baseFilename}.LAMMPS_logFile\n",
    "\n",
    "#(New) setup the force term by adding a group of random atoms with specific force (Constant force for simple case)\n",
    "#create_atoms 1 random 100 127569 box\n",
    "#group obstacle type 1\n",
    "#fix 4 force_term setforce 3.0 0.0 0.0\n",
    "#(/New)\n",
    "\n",
    "# == Setup style of the run\n",
    "\n",
    "# type of units to use in the simulation (units used are in fact: amu, nm, ns, Kelvins)\n",
    "units       nano\n",
    "\n",
    "# indicates possible types allowed for interactions between the atoms\n",
    "atom_style  angle \n",
    "\n",
    "# indicates possible types allowed for bonds between the atoms \n",
    "# bond_style  hybrid harmonic\n",
    "\n",
    "# indicates possible types allowed for bond angles between the atoms \n",
    "angle_style none\n",
    "\n",
    "# indicates type of boundary conditions in each direction (p = periodic) \n",
    "boundary p p p \n",
    "read_data ${baseFilename}.LAMMPS_read_data # file of atomic coordinates and topology\n",
    "velocity all zero linear                   # initialize all atomic velocities initially to zero\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "# feed commands to LAMMPs one line at a time\n",
    "print_log(\"Sending commands to LAMMPs\");\n",
    "for line in s.splitlines():\n",
    "  print_log(line);\n",
    "  L.command(line);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sending commands to LAMMPs\n",
      "\n",
      "# == Interactions \n",
      "pair_style none\n",
      "atom_modify sort 1000 ${neighborSkinDist}          # setup sort data explicitly since no interactions to set this data. \n",
      "\n",
      "# == Setup neighbor list distance\n",
      "comm_style tiled\n",
      "comm_modify mode single cutoff 202.0 vel yes\n",
      "\n",
      "neighbor ${neighborSkinDist} bin                    # first number gives a distance beyond the force cut-off ${neighborSkinDist}\n",
      "neigh_modify every 1\n",
      "atom_modify sort 0 ${neighborSkinDist}           # setup sort data explicitly since no interactions to set this data. \n",
      "\n",
      "\n",
      "\n",
      "# == Setup the SELM integrator\n",
      "fix 1 all selm Main.SELM_params\n",
      "LAMMPS (29 Oct 2020)\n",
      "Reading data file ...\n",
      "  triclinic box = (-202.00000 -202.00000 -202.00000) to (202.00000 202.00000 202.00000) with tilt (0.0000000 0.0000000 0.0000000)\n",
      "  1 by 1 by 1 MPI processor grid\n",
      "  reading atoms ...\n",
      "  2664 atoms\n",
      "Finding 1-2 1-3 1-4 neighbors ...\n",
      "  special bond factors lj:    0.0      0.0      0.0     \n",
      "  special bond factors coul:  0.0      0.0      0.0     \n",
      "     0 = max # of 1-2 neighbors\n",
      "     0 = max # of 1-3 neighbors\n",
      "     0 = max # of 1-4 neighbors\n",
      "     1 = max # of special neighbors\n",
      "  special bonds CPU = 0.005 seconds\n",
      "  read_data CPU = 0.014 seconds\n",
      "\n",
      "USER-SELM Package (SVN Version = 3389:3407, Compile Date-Time = Tue Jun 28 16:58:40 UTC 2022) \n",
      "For examples and additional information see http://mango-selm.org/\n",
      "Simulation Start Date-Time = 12-09-2023 22:14:34\n",
      "\n",
      "# note langevin just computes forces, nve integrates the motions\n",
      "#fix 1 all langevin 298.15 298.15 0.00001 48279\n",
      "#fix 2 all nve\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "s = \"\"\"\n",
    "# == Interactions \n",
    "pair_style none\n",
    "atom_modify sort 1000 ${neighborSkinDist}          # setup sort data explicitly since no interactions to set this data. \n",
    "\n",
    "# == Setup neighbor list distance\n",
    "comm_style tiled\n",
    "comm_modify mode single cutoff 202.0 vel yes\n",
    "\n",
    "neighbor ${neighborSkinDist} bin                    # first number gives a distance beyond the force cut-off ${neighborSkinDist}\n",
    "neigh_modify every 1\n",
    "atom_modify sort 0 ${neighborSkinDist}           # setup sort data explicitly since no interactions to set this data. \n",
    "\n",
    "\n",
    "\n",
    "# == Setup the SELM integrator\n",
    "fix 1 all selm Main.SELM_params\n",
    "\n",
    "# note langevin just computes forces, nve integrates the motions\n",
    "#fix 1 all langevin 298.15 298.15 0.00001 48279\n",
    "#fix 2 all nve\n",
    "\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "# feed commands to LAMMPs one line at a time\n",
    "print_log(\"Sending commands to LAMMPs\");\n",
    "for line in s.splitlines():\n",
    "  print_log(line);\n",
    "  L.command(line);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sending commands to LAMMPs\n",
      "\n",
      "# 2 obstacles\n",
      "\n",
      "region\t\tvoid1 sphere 10 4 0 70\n",
      "delete_atoms\tregion void1\n",
      "region\t\tvoid2 sphere 120 7 0 70\n",
      "delete_atoms\tregion void2\n"
     ]
    }
   ],
   "source": [
    "s=\"\"\"\n",
    "# 2 obstacles\n",
    "\n",
    "region\t\tvoid1 sphere 10 4 0 70\n",
    "delete_atoms\tregion void1\n",
    "region\t\tvoid2 sphere 120 7 0 70\n",
    "delete_atoms\tregion void2\n",
    "\"\"\"\n",
    "\n",
    "print_log(\"Sending commands to LAMMPs\");\n",
    "for line in s.splitlines():\n",
    "  print_log(line);\n",
    "  L.command(line);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sending commands to LAMMPs\n",
      "\n",
      "# define groups\n",
      "\n",
      "region 1 block INF INF INF -200 INF INF\n",
      "group lower_y region 1\n",
      "region 2 block INF -200 INF INF  INF INF\n",
      "group lower_x region 2\n",
      "region 3 block INF  INF INF  INF INF -200\n",
      "group lower_z region 3\n",
      "region 4 block INF INF 200 INF  INF INF\n",
      "group upper_y region 4\n",
      "region 5 block 200 INF INF INF  INF INF\n",
      "group upper_x region 5\n",
      "region 6 block INF  INF INF  INF 200 INF\n",
      "group upper_z region 6\n",
      "\n"
     ]
    }
   ],
   "source": [
    "s=\"\"\"\n",
    "# define groups\n",
    "\n",
    "region 1 block INF INF INF -200 INF INF\n",
    "group lower_y region 1\n",
    "region 2 block INF -200 INF INF  INF INF\n",
    "group lower_x region 2\n",
    "region 3 block INF  INF INF  INF INF -200\n",
    "group lower_z region 3\n",
    "region 4 block INF INF 200 INF  INF INF\n",
    "group upper_y region 4\n",
    "region 5 block 200 INF INF INF  INF INF\n",
    "group upper_x region 5\n",
    "region 6 block INF  INF INF  INF 200 INF\n",
    "group upper_z region 6\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "print_log(\"Sending commands to LAMMPs\");\n",
    "for line in s.splitlines():\n",
    "  print_log(line);\n",
    "  L.command(line);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x=np.random.uniform(3,97,10)\n",
    "y=np.random.uniform(3,97,10)\n",
    "z=np.random.uniform(3,97,10)\n",
    "\n",
    "points=list(zip(x,y,z))\n",
    "\n",
    "j=7\n",
    "for i in points:\n",
    "    L.region(j,\"sphere\",i[0],i[1],i[2],3)\n",
    "    j=j+1\n",
    "\n",
    "original=\"delete_atoms region \"\n",
    "for j in range(7,17):\n",
    "    line1=original+str(j)\n",
    "    line2=\"fix \"+ str(j) +\" flow indent \"+str(100)+\" sphere \"+ str(points[j-7][0])+\" \"+str(points[j-7][1])+\" \"+str(points[j-7][2])+\" \"+str(4)\n",
    "    L.command(line1)\n",
    "    L.command(line2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sending commands to LAMMPs\n",
      "\n",
      "# define groups\n",
      "\n",
      "group boundary union lower_y lower_x lower_z upper_y upper_x upper_z\n",
      "group flow subtract all boundary\n",
      "\n"
     ]
    }
   ],
   "source": [
    "s=\"\"\"\n",
    "# define groups\n",
    "\n",
    "group boundary union lower_y lower_x lower_z upper_y upper_x upper_z\n",
    "group flow subtract all boundary\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "print_log(\"Sending commands to LAMMPs\");\n",
    "for line in s.splitlines():\n",
    "  print_log(line);\n",
    "  L.command(line);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sending commands to LAMMPs\n",
      "\n",
      "fix 2 flow indent 100 sphere 10 4 0 4\n",
      "fix 3 flow indent 100 sphere 20 7 0 4\n"
     ]
    }
   ],
   "source": [
    "s=\"\"\"\n",
    "fix 2 flow indent 100 sphere 10 4 0 4\n",
    "fix 3 flow indent 100 sphere 20 7 0 4\n",
    "\"\"\"\n",
    "\n",
    "print_log(\"Sending commands to LAMMPs\");\n",
    "for line in s.splitlines():\n",
    "  print_log(line);\n",
    "  L.command(line);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sending commands to LAMMPs\n",
      "\n",
      "##(New) fix spring \n",
      "fix 4 all spring tether 5.0 0.0 0.0 0.0 5.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "s=\"\"\"\n",
    "##(New) fix spring \n",
    "fix 4 all spring tether 5.0 0.0 0.0 0.0 5.0\n",
    "\n",
    "\"\"\"\n",
    "#fix 5 flow spring tether 50.0 0.0 0.0 0.0 0.0\n",
    "\n",
    "print_log(\"Sending commands to LAMMPs\");\n",
    "for line in s.splitlines():\n",
    "  print_log(line);\n",
    "  L.command(line);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sending commands to LAMMPs\n",
      "\n",
      "group forcing type 3\n"
     ]
    }
   ],
   "source": [
    "s=\"\"\"\n",
    "group forcing type 3\n",
    "\"\"\"\n",
    "\n",
    "print_log(\"Sending commands to LAMMPs\");\n",
    "for line in s.splitlines():\n",
    "  print_log(line);\n",
    "  L.command(line);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sending commands to LAMMPs\n",
      "\n",
      "fix 5 forcing setforce 300.0 0.0 0.0\n"
     ]
    }
   ],
   "source": [
    "s=\"\"\"\n",
    "fix 5 forcing setforce 300.0 0.0 0.0\n",
    "\"\"\"\n",
    "\n",
    "print_log(\"Sending commands to LAMMPs\");\n",
    "for line in s.splitlines():\n",
    "  print_log(line);\n",
    "  L.command(line);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sending commands to LAMMPs\n",
      "\n",
      "info all out\n",
      "Making directory: ./ \n",
      "Making directory: .//sim_data \n",
      "WARNING: Ignoring 'compress yes' for molecular system (../delete_atoms.cpp:125)\n",
      "Deleted 56 atoms, new total = 2608\n",
      "WARNING: Ignoring 'compress yes' for molecular system (../delete_atoms.cpp:125)\n",
      "Deleted 53 atoms, new total = 2555\n",
      "0 atoms in group lower_y\n",
      "0 atoms in group lower_x\n",
      "0 atoms in group lower_z\n",
      "0 atoms in group upper_y\n",
      "0 atoms in group upper_x\n",
      "0 atoms in group upper_z\n",
      "0 atoms in group boundary\n",
      "2555 atoms in group flow\n",
      "1275 atoms in group forcing\n",
      "WARNING: Ignoring unknown or incorrect info command flag (../info.cpp:257)\n",
      "\n",
      "Info-Info-Info-Info-Info-Info-Info-Info-Info-Info-Info\n",
      "Printed on Sat Dec  9 22:14:36 2023\n",
      "\n",
      "\n",
      "LAMMPS version: 29 Oct 2020 / 20201029\n",
      "\n",
      "OS information: Linux 6.2.0-37-generic on x86_64\n",
      "\n",
      "sizeof(smallint): 32-bit\n",
      "sizeof(imageint): 32-bit\n",
      "sizeof(tagint):   32-bit\n",
      "sizeof(bigint):   64-bit\n",
      "\n",
      "Compiler: GNU C++ 6.3.0 20170516 with OpenMP not enabled\n",
      "C++ standard: C++14\n",
      "\n",
      "Active compile time flags:\n",
      "\n",
      "-DLAMMPS_GZIP\n",
      "-DLAMMPS_PNG\n",
      "-DLAMMPS_JPEG\n",
      "-DLAMMPS_FFMPEG\n",
      "-DLAMMPS_EXCEPTIONS\n",
      "-DLAMMPS_SMALLBIG\n",
      "\n",
      "Installed packages:\n",
      "\n",
      "MOLECULE USER-MISC USER-VTK \n",
      "\n",
      "Memory allocation information (MPI rank 0):\n",
      "\n",
      "Total dynamically allocated memory: 3.212 Mbyte\n",
      "Current reserved memory pool size: 39.69 Mbyte\n",
      "Maximum resident set size: 115.9 Mbyte\n",
      "\n",
      "Communication information:\n",
      "MPI library level: MPI v1.2\n",
      "MPI version: LAMMPS MPI STUBS for LAMMPS version 29 Oct 2020\n",
      "Comm style = tiled,  Comm layout = uniform\n",
      "Communicate velocities for ghost atoms = yes\n",
      "Communication mode = single\n",
      "Communication cutoff = 202.0\n",
      "Nprocs = 1,   Nthreads = 1\n",
      "Processor grid = 1 x 1 x 1\n",
      "\n",
      "System information:\n",
      "Units         = nano\n",
      "Atom style    = angle\n",
      "Atom map      = array\n",
      "Molecule type = standard\n",
      "Atoms     =         2555,  types =        3,  style = none\n",
      "Bonds     =            0,  types =        0,  style = none\n",
      "Angles    =            0,  types =        0,  style = none\n",
      "Dihedrals =            0,  types =        0,  style = none\n",
      "Impropers =            0,  types =        0,  style = none\n",
      "Special bond factors lj =    0.0      0.0      0.0     \n",
      "Special bond factors coul =  0.0      0.0      0.0     \n",
      "Kspace style = none\n",
      "\n",
      "Dimensions = 3\n",
      "Triclinic box = 404.00000 x 404.00000 x 404.00000\n",
      "Boundaries = p,p p,p p,p\n",
      "xlo, xhi = -202.00000, 202.00000\n",
      "ylo, yhi = -202.00000, 202.00000\n",
      "zlo, zhi = -202.00000, 202.00000\n",
      "Xy, xz, yz = 0.0000000, 0.0000000, 0.0000000\n",
      "\n",
      "Coeff status information:\n",
      "\n",
      "Group information:\n",
      "Group[ 0]:     all              (static)\n",
      "Group[ 1]:     lower_y          (static)\n",
      "Group[ 2]:     lower_x          (static)\n",
      "Group[ 3]:     lower_z          (static)\n",
      "Group[ 4]:     upper_y          (static)\n",
      "Group[ 5]:     upper_x          (static)\n",
      "Group[ 6]:     upper_z          (static)\n",
      "Group[ 7]:     boundary         (static)\n",
      "Group[ 8]:     flow             (static)\n",
      "Group[ 9]:     forcing          (static)\n",
      "\n",
      "Region information:\n",
      "Region[  0]:  void1,            style = sphere,           side = in\n",
      "   Boundary:  lo -60.000000 -66.000000 -70.000000  hi 80.000000 74.000000 70.000000\n",
      "Region[  1]:  void2,            style = sphere,           side = in\n",
      "   Boundary:  lo 50.000000 -63.000000 -70.000000  hi 190.00000 77.000000 70.000000\n",
      "Region[  2]:  1,                style = block,            side = in\n",
      "   Boundary:  lo -1.0000000e+20 -1.0000000e+20 -1.0000000e+20  hi 1.0000000e+20 -200.00000 1.0000000e+20\n",
      "Region[  3]:  2,                style = block,            side = in\n",
      "   Boundary:  lo -1.0000000e+20 -1.0000000e+20 -1.0000000e+20  hi -200.00000 1.0000000e+20 1.0000000e+20\n",
      "Region[  4]:  3,                style = block,            side = in\n",
      "   Boundary:  lo -1.0000000e+20 -1.0000000e+20 -1.0000000e+20  hi 1.0000000e+20 1.0000000e+20 -200.00000\n",
      "Region[  5]:  4,                style = block,            side = in\n",
      "   Boundary:  lo -1.0000000e+20 200.00000 -1.0000000e+20  hi 1.0000000e+20 1.0000000e+20 1.0000000e+20\n",
      "Region[  6]:  5,                style = block,            side = in\n",
      "   Boundary:  lo 200.00000 -1.0000000e+20 -1.0000000e+20  hi 1.0000000e+20 1.0000000e+20 1.0000000e+20\n",
      "Region[  7]:  6,                style = block,            side = in\n",
      "   Boundary:  lo -1.0000000e+20 -1.0000000e+20 200.00000  hi 1.0000000e+20 1.0000000e+20 1.0000000e+20\n",
      "\n",
      "Compute information:\n",
      "Compute[  0]:  thermo_temp,      style = temp,             group = all\n",
      "Compute[  1]:  thermo_press,     style = pressure,         group = all\n",
      "Compute[  2]:  thermo_pe,        style = pe,               group = all\n",
      "\n",
      "Dump information:\n",
      "\n",
      "Fix information:\n",
      "Fix[  0]:      1,                style = selm,             group = all\n",
      "Fix[  1]:      2,                style = indent,           group = flow\n",
      "Fix[  2]:      3,                style = indent,           group = flow\n",
      "Fix[  3]:      4,                style = spring,           group = all\n",
      "Fix[  4]:      5,                style = setforce,         group = forcing\n",
      "\n",
      "Variable information:\n",
      "Variable[  0]: dumpfreq,         style = equal,            def = 1 (undefined)\n",
      "Variable[  1]: restart,          style = equal,            def = 0 (undefined)\n",
      "Variable[  2]: neighborSkinDist,  style = equal,            def = 1.0 1\n",
      "Variable[  3]: baseFilename,     style = universe,         def = Model\n",
      "\n",
      "Total time information (MPI rank 0):\n",
      "  CPU time:    0:00:02\n",
      " Wall time:    0:00:02\n",
      "\n",
      "Styles information:\n",
      "\n",
      "Atom styles:\n",
      "\n",
      "angle           atomic          body            bond            charge          \n",
      "ellipsoid       full            hybrid          line            molecular       \n",
      "sphere          template        tri             \n",
      "\n",
      "\n",
      "\n",
      "Integrate styles:\n",
      "\n",
      "respa           verlet          \n",
      "\n",
      "\n",
      "\n",
      "Minimize styles:\n",
      "\n",
      "cg              fire            fire/old        hftn            quickmin        \n",
      "sd              \n",
      "\n",
      "\n",
      "\n",
      "Pair styles:\n",
      "\n",
      "agni            beck            born            born/coul/dsf   born/coul/wolf  \n",
      "buck            buck/coul/cut   buck/mdf        cosine/squared  coul/cut        \n",
      "coul/debye      coul/diel       coul/dsf        coul/shield     coul/slater/cut \n",
      "coul/slater/long                coul/streitz    coul/wolf       dpd             \n",
      "dpd/tstat       drip            e3b             edip            edip/multi      \n",
      "extep           gauss           gauss/cut       hbond/dreiding/lj               \n",
      "hbond/dreiding/morse            hybrid          hybrid/overlay  ilp/graphene/hbn                \n",
      "kolmogorov/crespi/full          kolmogorov/crespi/z             lebedeva/z      \n",
      "lennard/mdf     list            lj/charmm/coul/charmm           \n",
      "lj/charmm/coul/charmm/implicit  lj/charmmfsw/coul/charmmfsh     lj/cubic        \n",
      "lj/cut          lj/cut/coul/cut lj/cut/coul/debye               lj/cut/coul/dsf \n",
      "lj/cut/coul/wolf                lj/cut/tip4p/cut                lj/expand       \n",
      "lj/expand/coul/long             lj/gromacs      lj/gromacs/coul/gromacs         \n",
      "lj/mdf          lj/sf           lj/sf/dipole/sf lj/smooth       lj/smooth/linear                \n",
      "lj96/cut        local/density   meam/spline     meam/sw/spline  mie/cut         \n",
      "momb            morse           morse/smooth/linear             reax            \n",
      "soft            srp             table           tersoff/table   tip4p/cut       \n",
      "ufm             yukawa          zbl             zero            \n",
      "\n",
      "\n",
      "\n",
      "Bond styles:\n",
      "\n",
      "fene            fene/expand     gromos          harmonic        harmonic/shift  \n",
      "harmonic/shift/cut              hybrid          morse           nonlinear       \n",
      "quartic         special         table           zero            \n",
      "\n",
      "\n",
      "\n",
      "Angle styles:\n",
      "\n",
      "charmm          cosine          cosine/delta    cosine/periodic cosine/shift    \n",
      "cosine/shift/exp                cosine/squared  dipole          fourier         \n",
      "fourier/simple  harmonic        hybrid          quartic         table           \n",
      "zero            \n",
      "\n",
      "\n",
      "\n",
      "Dihedral styles:\n",
      "\n",
      "charmm          charmmfsw       cosine/shift/exp                fourier         \n",
      "harmonic        helix           hybrid          multi/harmonic  nharmonic       \n",
      "opls            quadratic       spherical       table           table/cut       \n",
      "zero            \n",
      "\n",
      "\n",
      "\n",
      "Improper styles:\n",
      "\n",
      "cossq           cvff            distance        fourier         harmonic        \n",
      "hybrid          ring            umbrella        zero            \n",
      "\n",
      "\n",
      "\n",
      "KSpace styles:\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "Fix styles:\n",
      "\n",
      "accelerate/cos  adapt           addforce    "
     ]
    }
   ],
   "source": [
    "s=\"\"\"\n",
    "info all out\n",
    "\"\"\"\n",
    "\n",
    "print_log(\"Sending commands to LAMMPs\");\n",
    "for line in s.splitlines():\n",
    "  print_log(line);\n",
    "  L.command(line);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sending commands to LAMMPs\n",
      "\n",
      "# == Setup output data write to disk\n",
      "dump        dmp_vtk all vtk ${dumpfreq} ./vtk/Particles_*.vtp id type vx fx\n",
      "dump_modify dmp_vtk pad 8 # ensures filenames file_000000.data\n",
      "\n",
      "# == simulation time-stepping\n",
      "timestep 6\n",
      "\n",
      "# == Run the simulation\n",
      "run      100 upto\n",
      "    addtorque       ave/atom        \n",
      "ave/chunk       ave/correlate   ave/correlate/long              ave/histo       \n",
      "ave/histo/weight                ave/spatial     ave/spatial/sphere              \n",
      "ave/time        aveforce        balance         box/relax       cmap            \n",
      "controller      deform          drag            dt/reset        \n",
      "electron/stopping               enforce2d       external        ffl             \n",
      "filter/corotate flow/gauss      gle             gravity         grem            \n",
      "halt            heat            imd             indent          ipi             \n",
      "langevin        lineforce       momentum        momentum/chunk  move            \n",
      "nph             nph/sphere      npt             npt/cauchy      npt/sphere      \n",
      "numdiff         nve             nve/limit       nve/noforce     nve/sphere      \n",
      "nvk             nvt             nvt/sllod       nvt/sphere      orient/eco      \n",
      "pafi            pimd            planeforce      press/berendsen print           \n",
      "propel/self     property/atom   recenter        restrain        rhok            \n",
      "selm            setforce        smd             spring          spring/chunk    \n",
      "spring/rg       spring/self     store/force     store/state     temp/berendsen  \n",
      "temp/csld       temp/csvr       temp/rescale    ti/spring       tmd             \n",
      "ttm/mod         vector          viscous         wall/ees        wall/harmonic   \n",
      "wall/lj1043     wall/lj126      wall/lj93       wall/morse      wall/reflect    \n",
      "wall/reflect/stochastic         wall/region     wall/region/ees \n",
      "\n",
      "\n",
      "\n",
      "Compute styles:\n",
      "\n",
      "ackland/atom    adf             aggregate/atom  angle           angle/local     \n",
      "angmom/chunk    basal/atom      bond            bond/local      centro/atom     \n",
      "centroid/stress/atom            chunk/atom      chunk/spread/atom               \n",
      "cluster/atom    cna/atom        cnp/atom        com             com/chunk       \n",
      "contact/atom    coord/atom      dihedral        dihedral/local  dipole/chunk    \n",
      "displace/atom   entropy/atom    erotate/sphere  erotate/sphere/atom             \n",
      "fragment/atom   global/atom     group/group     gyration        gyration/chunk  \n",
      "gyration/shape  gyration/shape/chunk            heat/flux       hexorder/atom   \n",
      "hma             improper        improper/local  inertia/chunk   ke              \n",
      "ke/atom         momentum        msd             msd/chunk       omega/chunk     \n",
      "orientorder/atom                pair            pair/local      pe              \n",
      "pe/atom         pressure        pressure/cylinder               property/atom   \n",
      "property/chunk  property/local  rdf             reduce          reduce/chunk    \n",
      "reduce/region   slice           stress/atom     stress/mop      \n",
      "stress/mop/profile              stress/mop_selm stress/mop_selm/profile         \n",
      "temp            temp/chunk      temp/com        temp/deform     temp/partial    \n",
      "temp/profile    temp/ramp       temp/region     temp/rotate     temp/sphere     \n",
      "torque/chunk    vacf            vcm/chunk       viscosity/cos   \n",
      "\n",
      "\n",
      "\n",
      "Region styles:\n",
      "\n",
      "block           cone            cylinder        intersect       plane           \n",
      "prism           sphere          union           \n",
      "\n",
      "\n",
      "\n",
      "Dump styles:\n",
      "\n",
      "atom            cfg             custom          dcd             image           \n",
      "local           movie           vtk             xyz             \n",
      "\n",
      "\n",
      "\n",
      "Command styles (add-on input script commands):\n",
      "\n",
      "balance         change_box      create_atoms    create_bonds    create_box      \n",
      "delete_atoms    delete_bonds    displace_atoms  info            minimize        \n",
      "read_data       read_dump       read_restart    replicate       rerun           \n",
      "reset_atom_ids  reset_ids       reset_mol_ids   run             set             \n",
      "temper/grem     temper/npt      velocity        write_coeff     write_data      \n",
      "write_dump      write_restart   \n",
      "\n",
      "\n",
      "\n",
      "Info-Info-Info-Info-Info-Info-Info-Info-Info-Info-Info\n",
      "\n",
      "Neighbor list info ...\n",
      "  update every 1 steps, delay 10 steps, check yes\n",
      "  max neighbors/atom: 2000, page size: 100000\n",
      "  master list distance cutoff = 0\n",
      "  ghost atom cutoff = 202\n",
      "  binsize = 404, bins = 1 1 1\n",
      "  0 neighbor lists, perpetual/occasio"
     ]
    },
    {
     "ename": "Exception",
     "evalue": "ERROR: Lost atoms: original 2555 current 2091 (../thermo.cpp:427)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mException\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[21], line 19\u001b[0m\n\u001b[1;32m     17\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m line \u001b[38;5;129;01min\u001b[39;00m s\u001b[38;5;241m.\u001b[39msplitlines():\n\u001b[1;32m     18\u001b[0m   print_log(line);\n\u001b[0;32m---> 19\u001b[0m   \u001b[43mL\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcommand\u001b[49m\u001b[43m(\u001b[49m\u001b[43mline\u001b[49m\u001b[43m)\u001b[49m;\n",
      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/selm_lammps/lammps.py:2523\u001b[0m, in \u001b[0;36mPyLammps.command\u001b[0;34m(self, cmd)\u001b[0m\n\u001b[1;32m   2513\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mcommand\u001b[39m(\u001b[38;5;28mself\u001b[39m, cmd):\n\u001b[1;32m   2514\u001b[0m \u001b[38;5;250m  \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m   2515\u001b[0m \u001b[38;5;124;03m  Execute LAMMPS command\u001b[39;00m\n\u001b[1;32m   2516\u001b[0m \n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m   2521\u001b[0m \u001b[38;5;124;03m  :type: cmd: string\u001b[39;00m\n\u001b[1;32m   2522\u001b[0m \u001b[38;5;124;03m  \"\"\"\u001b[39;00m\n\u001b[0;32m-> 2523\u001b[0m   \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlmp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcommand\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcmd\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   2524\u001b[0m   \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_cmd_history\u001b[38;5;241m.\u001b[39mappend(cmd)\n",
      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/selm_lammps/lammps.py:651\u001b[0m, in \u001b[0;36mlammps.command\u001b[0;34m(self, cmd)\u001b[0m\n\u001b[1;32m    648\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlib\u001b[38;5;241m.\u001b[39mlammps_command(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlmp,cmd)\n\u001b[1;32m    650\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhas_exceptions \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlib\u001b[38;5;241m.\u001b[39mlammps_has_error(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlmp):\n\u001b[0;32m--> 651\u001b[0m   \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_lammps_exception\n",
      "\u001b[0;31mException\u001b[0m: ERROR: Lost atoms: original 2555 current 2091 (../thermo.cpp:427)"
     ]
    }
   ],
   "source": [
    "s=\"\"\"\n",
    "# == Setup output data write to disk\n",
    "dump        dmp_vtk all vtk ${dumpfreq} ./vtk/Particles_*.vtp id type vx fx\n",
    "dump_modify dmp_vtk pad 8 # ensures filenames file_000000.data\n",
    "\n",
    "# == simulation time-stepping\n",
    "timestep 6\n",
    "\n",
    "# == Run the simulation\n",
    "run      100 upto\n",
    "\n",
    "# == Write restart data\n",
    "write_restart ${baseFilename}.LAMMPS_restart_data\n",
    "\"\"\"\n",
    "# feed commands to LAMMPs one line at a time\n",
    "print_log(\"Sending commands to LAMMPs\");\n",
    "for line in s.splitlines():\n",
    "  print_log(line);\n",
    "  L.command(line);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!cat Model.SELM_Info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done\n"
     ]
    }
   ],
   "source": [
    "print_log(\"Done\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  },
  "vscode": {
   "interpreter": {
    "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}