simulation_changed_2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "39f53e51-6e85-4378-8d70-bc9ab0e9ac2b",
   "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",
      "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",
      "Model files being copied:\n",
      "src = /home/xsu2/Downloads/hydro/Model1\n",
      "dst = /home/xsu2/Downloads/hydro/output/simulation_hydro1/hydro_test001/\n",
      "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",
      "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",
      "Writing model data .pickle\n",
      "filename = model_setup.pickle\n",
      "Writing model data .read_data\n",
      "filename = Model.LAMMPS_read_data\n"
     ]
    }
   ],
   "source": [
    "#!/usr/bin/env python\n",
    "# coding: utf-8\n",
    "\n",
    "# ### SELM via PyLAMMPs for Simulations\n",
    "# Author: Paul Atzberger <br>\n",
    "# http://atzberger.org/\n",
    "# \n",
    "\n",
    "# In[1]:\n",
    "\n",
    "\n",
    "import os;\n",
    "script_base_name = \"simulation_hydro1\";\n",
    "script_dir = os.getcwd();\n",
    "\n",
    "\n",
    "# In[2]:\n",
    "\n",
    "\n",
    "# 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;\n",
    "\n",
    "\n",
    "# ### Setup SELM Simulation\n",
    "\n",
    "# In[3]:\n",
    "\n",
    "\n",
    "# @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",
    "\n",
    "# ### Setup LAMMPs\n",
    "\n",
    "# In[4]:\n",
    "\n",
    "\n",
    "L = IPyLammps();\n",
    "atz_util.print_version_info(L);    \n",
    "\n",
    "\n",
    "# ### Copy files to the output directory\n",
    "\n",
    "# In[5]:\n",
    "\n",
    "\n",
    "# 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));\n",
    "\n",
    "\n",
    "# In[6]:\n",
    "\n",
    "\n",
    "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));\n",
    "\n",
    "\n",
    "# # Common Physical Parameters (nano units)\n",
    "\n",
    "# In[7]:\n",
    "\n",
    "\n",
    "# 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']});\n",
    "\n",
    "\n",
    "# Set up Lammps\n",
    "\n",
    "# ### Setup the Model and Simulation Files (such as .read_data)\n",
    "\n",
    "# In[8]:\n",
    "\n",
    "\n",
    "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});\n",
    "\n",
    "\n",
    "# In[9]:\n",
    "\n",
    "\n",
    "# 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);\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "960076a8-a13f-437c-bb27-6ef29e909e64",
   "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": [
    "# In[10]:\n",
    "\n",
    "\n",
    "#!cat Model.LAMMPS_read_data\n",
    "\n",
    "\n",
    "# ### Perform the simulation\n",
    "\n",
    "# In[11]:\n",
    "\n",
    "\n",
    "# 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": 3,
   "id": "20d0c515-ab2e-4dfb-ba42-bbc167fef829",
   "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.006 seconds\n",
      "  read_data CPU = 0.068 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 = 01-03-2024 22:42:21\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": [
    "# In[12]:\n",
    "\n",
    "\n",
    "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);\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9fa32dc7-0219-4ffa-a9d1-5dfa8722e029",
   "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": [
    "# In[14]:\n",
    "\n",
    "\n",
    "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": 5,
   "id": "9663acd4-fc4c-468f-b638-7159bcaceedc",
   "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": [
    "# In[15]:\n",
    "\n",
    "\n",
    "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": 6,
   "id": "30358068-95f4-4c4e-a28a-1a5c26d30f71",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "delete_atoms region 7\n",
      "fix 7 flow indent 50 sphere 93.78038903435228 -89.72223612055441 152.58549857746436 41\n",
      "delete_atoms region 8\n",
      "fix 8 flow indent 50 sphere 24.201880519590702 11.625118128709033 -3.5541762631515326 41\n",
      "delete_atoms region 9\n",
      "fix 9 flow indent 50 sphere -81.91439217975162 -80.80059712799749 -119.88832420677906 41\n",
      "delete_atoms region 10\n",
      "fix 10 flow indent 50 sphere -24.011448152552674 -8.613257454153029 -45.31507042727988 41\n",
      "delete_atoms region 11\n",
      "fix 11 flow indent 50 sphere 117.16402816505939 75.51836082103563 64.53048153602847 41\n",
      "delete_atoms region 12\n",
      "fix 12 flow indent 50 sphere -159.82548923837714 -22.020670951971795 121.71142663363139 41\n",
      "delete_atoms region 13\n",
      "fix 13 flow indent 50 sphere -122.88580657606836 69.79822636490226 74.96398718278817 41\n",
      "delete_atoms region 14\n",
      "fix 14 flow indent 50 sphere -32.48186117255716 -44.077001982571005 -124.92032314356344 41\n",
      "delete_atoms region 15\n",
      "fix 15 flow indent 50 sphere 146.2725340914957 -20.118609700566708 -93.66516838404152 41\n",
      "delete_atoms region 16\n",
      "fix 16 flow indent 50 sphere 149.09915722987017 -5.997246266992363 -56.60976411421106 41\n"
     ]
    }
   ],
   "source": [
    "x=np.random.uniform(-160,160,10)\n",
    "y=np.random.uniform(-160,160,10)\n",
    "z=np.random.uniform(-160,160,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],40)\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(50)+\" sphere \"+ str(points[j-7][0])+\" \"+str(points[j-7][1])+\" \"+str(points[j-7][2])+\" \"+str(41)\n",
    "    L.command(line1)\n",
    "    print_log(line1)\n",
    "    L.command(line2)\n",
    "    print_log(line2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "118dac29-92ca-48f6-81b6-e29a4e378a90",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(93.78038903435228, -89.72223612055441, 152.58549857746436),\n",
       " (24.201880519590702, 11.625118128709033, -3.5541762631515326),\n",
       " (-81.91439217975162, -80.80059712799749, -119.88832420677906),\n",
       " (-24.011448152552674, -8.613257454153029, -45.31507042727988),\n",
       " (117.16402816505939, 75.51836082103563, 64.53048153602847),\n",
       " (-159.82548923837714, -22.020670951971795, 121.71142663363139),\n",
       " (-122.88580657606836, 69.79822636490226, 74.96398718278817),\n",
       " (-32.48186117255716, -44.077001982571005, -124.92032314356344),\n",
       " (146.2725340914957, -20.118609700566708, -93.66516838404152),\n",
       " (149.09915722987017, -5.997246266992363, -56.60976411421106)]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a6636b53-e9a8-47ff-b60a-fff74d6442a0",
   "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",
    "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": 9,
   "id": "e1d3423f-5fff-421a-b7b8-94787ed6efcb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sending commands to LAMMPs\n",
      "\n",
      "group forcing type 3\n",
      "Sending commands to LAMMPs\n",
      "\n",
      "fix 5 forcing setforce 20.0 0.0 0.0\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);\n",
    "\n",
    "\n",
    "# In[19]:\n",
    "\n",
    "\n",
    "s=\"\"\"\n",
    "fix 5 forcing setforce 20.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);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f8379b76-7586-43b2-bfaa-d08cc8cfce7b",
   "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",
      "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",
      "2664 atoms in group flow\n",
      "WARNING: Ignoring 'compress yes' for molecular system (../delete_atoms.cpp:125)\n",
      "Deleted 11 atoms, new total = 2653\n",
      "WARNING: Ignoring 'compress yes' for molecular system (../delete_atoms.cpp:125)\n",
      "Deleted 11 atoms, new total = 2642\n",
      "WARNING: Ignoring 'compress yes' for molecular system (../delete_atoms.cpp:125)\n",
      "Deleted 11 atoms, new total = 2631\n",
      "WARNING: Ignoring 'compress yes' for molecular system (../delete_atoms.cpp:125)\n",
      "Deleted 10 atoms, new total = 2621\n",
      "WARNING: Ignoring 'compress yes' for molecular system (../delete_atoms.cpp:125)\n",
      "Deleted 9 atoms, new total = 2612\n",
      "WARNING: Ignoring 'compress yes' for molecular system (../delete_atoms.cpp:125)\n",
      "Deleted 13 atoms, new total = 2599\n",
      "WARNING: Ignoring 'compress yes' for molecular system (../delete_atoms.cpp:125)\n",
      "Deleted 8 atoms, new total = 2591\n",
      "WARNING: Ignoring 'compress yes' for molecular system (../delete_atoms.cpp:125)\n",
      "Deleted 10 atoms, new total = 2581\n",
      "WARNING: Ignoring 'compress yes' for molecular system (../delete_atoms.cpp:125)\n",
      "Deleted 9 atoms, new total = 2572\n",
      "WARNING: Ignoring 'compress yes' for molecular system (../delete_atoms.cpp:125)\n",
      "Deleted 6 atoms, new total = 2566\n",
      "1280 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 Wed Jan  3 22:46:43 2024\n",
      "\n",
      "\n",
      "LAMMPS version: 29 Oct 2020 / 20201029\n",
      "\n",
      "OS information: Linux 6.2.0-39-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.63 Mbyte\n",
      "Maximum resident set size: 115.3 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     =         2566,  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]:  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[  1]:  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[  2]:  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[  3]:  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[  4]:  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[  5]:  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",
      "Region[  6]:  7,                style = sphere,           side = in\n",
      "   Boundary:  lo 53.780389 -129.72224 112.58550  hi 133.78039 -49.722236 192.58550\n",
      "Region[  7]:  8,                style = sphere,           side = in\n",
      "   Boundary:  lo -15.798119 -28.374882 -43.554176  hi 64.201881 51.625118 36.445824\n",
      "Region[  8]:  9,                style = sphere,           side = in\n",
      "   Boundary:  lo -121.91439 -120.80060 -159.88832  hi -41.914392 -40.800597 -79.888324\n",
      "Region[  9]:  10,               style = sphere,           side = in\n",
      "   Boundary:  lo -64.011448 -48.613257 -85.315070  hi 15.988552 31.386743 -5.3150704\n",
      "Region[ 10]:  11,               style = sphere,           side = in\n",
      "   Boundary:  lo 77.164028 35.518361 24.530482  hi 157.16403 115.51836 104.53048\n",
      "Region[ 11]:  12,               style = sphere,           side = in\n",
      "   Boundary:  lo -199.82549 -62.020671 81.711427  hi -119.82549 17.979329 161.71143\n",
      "Region[ 12]:  13,               style = sphere,           side = in\n",
      "   Boundary:  lo -162.88581 29.798226 34.963987  hi -82.885807 109.79823 114.96399\n",
      "Region[ 13]:  14,               style = sphere,           side = in\n",
      "   Boundary:  lo -72.481861 -84.077002 -164.92032  hi 7.5181388 -4.0770020 -84.920323\n",
      "Region[ 14]:  15,               style = sphere,           side = in\n",
      "   Boundary:  lo 106.27253 -60.118610 -133.66517  hi 186.27253 19.881390 -53.665168\n",
      "Region[ 15]:  16,               style = sphere,           side = in\n",
      "   Boundary:  lo 109.09916 -45.997246 -96.609764  hi 189.09916 34.002754 -16.609764\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]:      7,                style = indent,           group = flow\n",
      "Fix[  2]:      8,                style = indent,           group = flow\n",
      "Fix[  3]:      9,                style = indent,           group = flow\n",
      "Fix[  4]:      10,               style = indent,           group = flow\n",
      "Fix[  5]:      11,               style = indent,           group = flow\n",
      "Fix[  6]:      12,               style = indent,           group = flow\n",
      "Fix[  7]:      13,               style = indent,           group = flow\n",
      "Fix[  8]:      14,               style = indent,           group = flow\n",
      "Fix[  9]:      15,               style = indent,           group = flow\n",
      "Fix[ 10]:      16,               style = indent,           group = flow\n",
      "Fix[ 11]:      4,                style = spring,           group = all\n",
      "Fix[ 12]:      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:03\n",
      " Wall time:    0:05:24\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        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/"
     ]
    }
   ],
   "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);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "5bf6b2c7-2138-45e2-b934-07529f895465",
   "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",
      "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/\n",
      "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/occasional/extra = 0 0 0\n",
      "Setting up Verlet run ...\n",
      "  Unit style    : nano\n",
      "  Current step  : 0\n",
      "  Time step     : 6.0\n",
      "Per MPI rank memory allocation (min/avg/max) = 8.104 | 8.104 | 8.104 Mbytes\n",
      "Step Temp E_pair E_mol TotEng Press Volume \n",
      "       0            0            0            0            0            0     65939264 \n",
      "     100   0.14164409            0            0     7.524208 7.6072106e-08     65939264 \n",
      "Loop time of 60.8502 on 1 procs for 100 steps with 2566 atoms\n",
      "\n",
      "Performance: 851928.819 ns/day, 0.000 hours/ns, 1.643 timesteps/s\n",
      "95.2% CPU use with 1 MPI tasks x no OpenMP threads\n",
      "\n",
      "MPI task timing breakdown:\n",
      "Section |  min time  |  avg time  |  max time  |%varavg| %total\n",
      "---------------------------------------------------------------\n",
      "Pair    | 0          | 0          | 0          |   0.0 |  0.00\n",
      "Bond    | 3.2425e-05 | 3.2425e-05 | 3.2425e-05 |   0.0 |  0.00\n",
      "Neigh   | 0.00029302 | 0.00029302 | 0.00029302 |   0.0 |  0.00\n",
      "Comm    | 0.13872    | 0.13872    | 0.13872    |   0.0 |  0.23\n",
      "Output  | 2.1302     | 2.1302     | 2.1302     |   0.0 |  3.50\n",
      "Modify  | 58.573     | 58.573     | 58.573     |   0.0 | 96.26\n",
      "Other   |            | 0.007928   |            |       |  0.01\n",
      "\n",
      "Nlocal:        2566.00 ave        2566 max        2566 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "Nghost:        17962.0 ave       17962 max       17962 min\n",
      "Histogram: 1 0 0 0 0 0 0 0 0 0\n",
      "Neighs:        0.00000# == Write restart data\n",
      "write_restart ${baseFilename}.LAMMPS_restart_data\n"
     ]
    }
   ],
   "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);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "108258ea-40f4-4639-87e5-c606ae7f8358",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n",
      "<SELM_Info>\n",
      "<LAMMPS_Version value=\"29 Oct 2020\"/>\n",
      "<SELM_SVN_Version value=\"3389:3407\"/>\n",
      "<SELM_Compile_Date_Time value=\"Tue Jun 28 16:58:40 UTC 2022\"/>\n",
      "<Simulation_Start_Date_Time value=\"01-03-2024 22:42:22\"/>\n",
      "</SELM_Info>\n"
     ]
    }
   ],
   "source": [
    "!cat Model.SELM_Info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0441c61-bd66-4951-9650-c3dfca502b8a",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}