Simulations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "98f83939",
   "metadata": {},
   "source": [
    "# Free-Energy Decomposition of Salt Effects on the Solubilities of Small Molecules and the Role of Excluded-Volume Effects\n",
    "Stefan Hervø-Hansen<sup>a,&#42;</sup>, Lin Daoyang<sup>a</sup>, Kento Kasahara<sup>a</sup>, and Nobuyuki Matubayasi<sup>a,&#42;</sup>.<br><br>\n",
    "<sup>a</sup> Division of Chemical Engineering, Graduate School of Engineering Science, Osaka University, Toyonaka, Osaka 560-8531, Japan.<br>\n",
    "<sup>*</sup> To whom correspondence may be addressed: stefan@cheng.es.osaka-u.ac.jp and nobuyuki@cheng.es.osaka-u.ac.jp.\n",
    "\n",
    "## PART 2: Simulations \n",
    "\n",
    "## Import of Python Modules & Auxiliary Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e5c7f80",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Notebook dependent libs\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os, re, time\n",
    "from shutil import copyfile\n",
    "from distutils.spawn import find_executable\n",
    "import parmed as pmd\n",
    "from openmm import app\n",
    "import openmm as mm\n",
    "import mdtraj as md\n",
    "\n",
    "# Simulation specific libs\n",
    "# from openmm import app\n",
    "# import openmm as mm\n",
    "# import mdtraj as md\n",
    "\n",
    "# Check for external programs\n",
    "if None in [find_executable('packmol')]:\n",
    "    print('WARNING: External program missing!')\n",
    "    \n",
    "# Physical constants & useful conversions\n",
    "Na = 6.02214076e23 # Avogadro constant [mol]\n",
    "liter_to_cubeangstrom = 1e27\n",
    "\n",
    "\n",
    "homedir = !pwd\n",
    "homedir = homedir[0]\n",
    "workdir = '/home/users/eau/SI-Cavity-Formation'\n",
    "print(homedir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2a0c580",
   "metadata": {},
   "source": [
    "## Molecular Dynamics Simulations\n",
    "### Simulation settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c31d1cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "states = { # State of simulations, (outFreq is steps per frame)\n",
    "          'vac': {'Nsteps':   1000000, 'OutFreq':  10}, #   1 nanoseconds,  100000 frames\n",
    "          'sol': {'Nsteps':  10000000, 'OutFreq': 100}, #  20 nanoseconds,  100000 frames\n",
    "          'ref': {'Nsteps':  25000000, 'OutFreq': 100}, #  50 nanoseconds,  250000 frames\n",
    "          'rdf': {'Nsteps':  10000000, 'OutFreq': 100}, #  20 nanoseconds,  100000 frames\n",
    "\n",
    "}\n",
    "\n",
    "salts = {\n",
    "    'NaF':  {'cation': 'Na', 'anion': 'F' },\n",
    "    'NaCl': {'cation': 'Na', 'anion': 'Cl'},\n",
    "    'NaI':  {'cation': 'Na', 'anion': 'I' },\n",
    "    'KF':   {'cation': 'K',  'anion': 'F' },\n",
    "    'KCl':  {'cation': 'K',  'anion': 'Cl'},\n",
    "    'KI':   {'cation': 'K',  'anion': 'I' },\n",
    "    'CsF':  {'cation': 'Cs', 'anion': 'F' },\n",
    "    'CsCl': {'cation': 'Cs', 'anion': 'Cl'},\n",
    "    'CsI':  {'cation': 'Cs', 'anion': 'I' },\n",
    "}\n",
    "\n",
    "solutes = {\n",
    "    'methanol':   {'SMILES': 'CO',      'Category': 'n-alcohols', 'Resname': 'MOH'},\n",
    "    'ethanol' :   {'SMILES': 'CCO',     'Category': 'n-alcohols', 'Resname': 'EOH'},\n",
    "    'n-propanol': {'SMILES': 'CCCO',    'Category': 'n-alcohols', 'Resname': 'PHO'},\n",
    "    'n-butanol' : {'SMILES': 'CCCCO',   'Category': 'n-alcohols', 'Resname': 'BOH'},\n",
    "    'n-pentanol': {'SMILES': 'CCCCCO',  'Category': 'n-alcohols', 'Resname': 'POH'},\n",
    "    'n-hexanol':  {'SMILES': 'CCCCCCO', 'Category': 'n-alcohols', 'Resname': 'HHO'},\n",
    "    \n",
    "    'methane':    {'SMILES': 'C',       'Category': 'n-alkanes',  'Resname': 'MTH'},\n",
    "    'ethane' :    {'SMILES': 'CC',      'Category': 'n-alkanes',  'Resname': 'ETH'},\n",
    "    'propane':    {'SMILES': 'CCC',     'Category': 'n-alkanes',  'Resname': 'PRP'},\n",
    "    'n-butane' :  {'SMILES': 'CCCC',    'Category': 'n-alkanes',  'Resname': 'BUT'},\n",
    "    'n-pentane':  {'SMILES': 'CCCCC',   'Category': 'n-alkanes',  'Resname': 'PEN'},\n",
    "    'n-hexane':   {'SMILES': 'CCCCCC',  'Category': 'n-alkanes',  'Resname': 'HEX'},\n",
    "}\n",
    "\n",
    "concentrations = {\n",
    "     0.00: {'N_solute':1, 'N_water':10000, 'N_cation':0,   'N_anion':0},\n",
    "     0.50: {'N_solute':1, 'N_water':10000, 'N_cation':90,  'N_anion':90},\n",
    "     1.00: {'N_solute':1, 'N_water':10000, 'N_cation':180, 'N_anion':180}\n",
    "}\n",
    "\n",
    "GENERATE_SOLUTES    = False  # Generate the structual ensemble of a solute in vacuum\n",
    "GENERATE_REFERENCES = False  # Generate the reference conditions\n",
    "GENERATE_SOLUTIONS  = False  # Generate the solution conditions\n",
    "GENERATE_RDF        = False  # Generate simulations for RDF calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60386750",
   "metadata": {},
   "source": [
    "## Step 1. Generation of Vacuum Solute Configurations\n",
    "### Construction of initial configurations and topology file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dae8a118",
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd -q $homedir\n",
    "\n",
    "packmol_script=\"\"\"\n",
    "tolerance 2.0\n",
    "filetype pdb\n",
    "output Initial.pdb\n",
    "\n",
    "add_box_sides 1.0\n",
    "\n",
    "structure {homedir}/PDB_files/{category}/{solute}.pdb\n",
    "        number 1\n",
    "        fixed 10. 10. 10. 0. 0. 0.\n",
    "        centerofmass\n",
    "end structure\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "topology_file = \"\"\"\n",
    "; Include forcefield files\n",
    "#include \"{homedir}/Force_fields/{solute}.itp\"\n",
    "\n",
    "[ defaults ]\n",
    "; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ\n",
    "1               2               no              1            0.83333333\n",
    "\n",
    "[ system ]\n",
    "{solute} in vacuum\n",
    "\n",
    "[ molecules ]\n",
    ";molecule name   nr.\n",
    "{solute_resname}             1\n",
    "\"\"\"\n",
    "\n",
    "if GENERATE_SOLUTES:\n",
    "    for solute, solute_details in solutes.items():\n",
    "        %cd -q $homedir/Simulations/$solute/Solute\n",
    "          \n",
    "        # Packmol input\n",
    "        with open('packmol.in', 'w') as f:\n",
    "            f.write(packmol_script.format(homedir=homedir, category=solute_details['Category'], solute=solute))\n",
    "        f.close()\n",
    "        !packmol < packmol.in\n",
    "        \n",
    "        # Topology input\n",
    "        with open('Topology.top', 'w') as f:\n",
    "            f.write(topology_file.format(homedir=workdir, solute=solute, solute_resname=solute_details['Resname']))\n",
    "        f.close()\n",
    "        \n",
    "        print('Wrote initial configurations and topology files to'+os.getcwd())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63ba023a",
   "metadata": {},
   "source": [
    "### Molecular dynamics setup using OpenMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "611dd8b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd -q $homedir\n",
    "\n",
    "openmm_script=\"\"\"# Imports\n",
    "import sys\n",
    "import os\n",
    "import openmm as mm\n",
    "from openmm import app\n",
    "from openmm import unit as u\n",
    "from mdtraj.reporters import XTCReporter\n",
    "\n",
    "print('Loading initial configuration and toplogy')\n",
    "pdb = app.PDBFile('Initial.pdb')\n",
    "top = app.GromacsTopFile('Topology.top', periodicBoxVectors=pdb.topology.getPeriodicBoxVectors())\n",
    "\n",
    "# Creating system\n",
    "print('Creating OpenMM System')\n",
    "system = top.createSystem(nonbondedMethod=app.NoCutoff, constraints=app.HBonds,\n",
    "                          rigidWater=True, removeCMMotion=True)\n",
    "        \n",
    "# Temperature-coupling by leap frog (BAOAB) Langevin integrator (NVT)\n",
    "integrator = mm.LangevinMiddleIntegrator(298.15*u.kelvin, 1.0/u.picoseconds, 1.0*u.femtoseconds)\n",
    "\n",
    "platform = mm.Platform.getPlatformByName('CUDA')\n",
    "properties = {{'CudaPrecision': 'mixed', 'CudaDeviceIndex': '0'}}\n",
    "\n",
    "# Create the Simulation object\n",
    "sim = app.Simulation(top.topology, system, integrator, platform, properties)\n",
    "\n",
    "# Set the particle positions\n",
    "sim.context.setPositions(pdb.positions)\n",
    "\n",
    "# Minimize the energy\n",
    "print('Minimizing energy')\n",
    "sim.minimizeEnergy()\n",
    "    \n",
    "# Draw initial MB velocities\n",
    "sim.context.setVelocitiesToTemperature(298.15*u.kelvin)\n",
    "\n",
    "# Equlibrate simulation\n",
    "print('Equilibrating...')\n",
    "sim.step(500000)  # 500000*1 fs = 0.5 ns\n",
    "\n",
    "# Set up the reporters\n",
    "sim.reporters.append(app.StateDataReporter('output.dat', {outFreq}, totalSteps={Nsteps}+500000,\n",
    "    time=True, potentialEnergy=True, kineticEnergy=True, temperature=True, volume=True, density=True,\n",
    "    remainingTime=True, speed=True, separator='\\t'))\n",
    "\n",
    "# Set up trajectory reporter\n",
    "sim.reporters.append(XTCReporter('trajectory.xtc', reportInterval={outFreq}, append=False))\n",
    "\n",
    "# Run dynamics\n",
    "print('Running dynamics! (NPT)')\n",
    "sim.step({Nsteps})\n",
    "\"\"\"           \n",
    "\n",
    "if GENERATE_SOLUTES:\n",
    "    N_simulations=0\n",
    "    for solute, solute_details in solutes.items():\n",
    "        %cd -q $homedir/Simulations/$solute/Solute\n",
    "        \n",
    "        with open('openMM.py', 'w') as f:\n",
    "            f.write(openmm_script.format(Nsteps=states['vac']['Nsteps'], outFreq=states['vac']['OutFreq']))\n",
    "        f.close()\n",
    "        print('Wrote run_openMM.py files to '+os.getcwd())\n",
    "        N_simulations+=1\n",
    "\n",
    "    print('Simulations about to be submitted: {}'.format(N_simulations))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58e1b181",
   "metadata": {},
   "source": [
    "### Submit script"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6586f241",
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd -q $homedir\n",
    "\n",
    "submit_script=\"\"\"#!/bin/bash\n",
    "#PBS -l select=1:ncpus=1:mpiprocs=1:ompthreads=1:jobtype=gpu:ngpus=1\n",
    "#PBS -l walltime=4:00:00\n",
    "#PBS -N Solute\n",
    "#PBS -e run.err\n",
    "#PBS -o run.out\n",
    "\n",
    "source ~/.bashrc\n",
    "source ~/.bash_profile\n",
    "\n",
    "cd {path}\n",
    "\n",
    "python openMM.py\n",
    "\"\"\"\n",
    "\n",
    "master_submit_script=\"\"\"#!/bin/bash\n",
    "read -p \"Will submit {N_simulations} simulations. Do you want to proceed? (yes/no) \" yn\n",
    "\n",
    "case $yn in \n",
    "    yes ) echo submitting...;;\n",
    "    no ) echo exiting...;\n",
    "         exit;;\n",
    "    * ) echo invalid response;\n",
    "        exit 1;;\n",
    "esac\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "if GENERATE_SOLUTES:     \n",
    "    with open('Simulations/master_solute.sh', 'w') as ff:\n",
    "        ff.write(master_submit_script.format(N_simulations=N_simulations))\n",
    "        \n",
    "        for solute in solutes:\n",
    "            %cd -q $homedir/Simulations/$solute/Solute\n",
    "            \n",
    "            simpath = os.getcwd()\n",
    "            simpath = simpath.split(sep='SI-Cavity-Formation')[1]\n",
    "            simpath = workdir+simpath\n",
    "            with open('submit.pbs', 'w') as f:\n",
    "                f.write(submit_script.format(solute=solute, path=simpath))\n",
    "            f.close()\n",
    "            s = 'cd {}\\njsub submit.pbs\\n\\n'.format(simpath)\n",
    "            ff.write(s)\n",
    "    ff.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9aa8dda",
   "metadata": {},
   "source": [
    "## Step 2. Generation of Solvent Configurations - Reference System (λ = 0)\n",
    "### Construction of initial configurations and topology file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a40c17c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd -q $homedir\n",
    "\n",
    "packmol_script=\"\"\"\n",
    "tolerance 2.0\n",
    "filetype pdb\n",
    "output Initial.pdb\n",
    "\n",
    "add_box_sides 0.1\n",
    "\n",
    "{saltfix}structure {homedir}/PDB_files/Ions/{cation}.pdb\n",
    "{saltfix}        number {N_cation}\n",
    "{saltfix}        inside cube 0. 0. 0. 67.\n",
    "{saltfix}end structure\n",
    "        \n",
    "{saltfix}structure {homedir}/PDB_files/Ions/{anion}.pdb\n",
    "{saltfix}        number {N_anion}\n",
    "{saltfix}        inside cube 0. 0. 0. 67.\n",
    "{saltfix}end structure\n",
    "\n",
    "structure {homedir}/PDB_files/spce.pdb\n",
    "        number {N_water}\n",
    "        inside cube 0. 0. 0. 67.\n",
    "end structure\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "topology_file = \"\"\"\n",
    "; Include forcefield files\n",
    "#include \"{homedir}/Force_fields/ions.itp\"\n",
    "#include \"{homedir}/Force_fields/spce.itp\"\n",
    "\n",
    "[ defaults ]\n",
    "; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ\n",
    "1               2               no              1            0.83333333\n",
    "\n",
    "[ system ]\n",
    "Salt water\n",
    "\n",
    "[ molecules ]\n",
    ";molecule name   nr.\n",
    "{saltfix}{cation_resname}          {N_cation}\n",
    "{saltfix}{anion_resname}          {N_anion}\n",
    "SOL          {N_water}\n",
    "\"\"\"\n",
    "\n",
    "if GENERATE_REFERENCES:\n",
    "    for i, (salt, salt_details) in enumerate(salts.items()):\n",
    "        for conc, conc_details in concentrations.items():\n",
    "            if conc == 0 and i == 0:\n",
    "                %cd -q $homedir/Simulations/References/No_salt\n",
    "            elif conc == 0 and i > 0:\n",
    "                continue\n",
    "            else:\n",
    "                conc = '{0:.2f}'.format(conc)\n",
    "                %cd -q $homedir/Simulations/References/$salt/$conc\n",
    "                \n",
    "            # Construction of initial configuration\n",
    "            with open('packmol.in', 'w') as f:\n",
    "                # Fix for no salt\n",
    "                if float(conc) == 0:\n",
    "                    saltfix='#'\n",
    "                else:\n",
    "                    saltfix=''\n",
    "                f.write(packmol_script.format(homedir=homedir, saltfix=saltfix, N_water=conc_details['N_water'],\n",
    "                                              cation=salt_details['cation'], anion=salt_details['anion'],\n",
    "                                              N_cation=conc_details['N_cation'], N_anion=conc_details['N_anion']))\n",
    "            f.close()\n",
    "            !packmol < packmol.in\n",
    "            \n",
    "            # Construction of topology file\n",
    "            with open('Topology.top', 'w') as f:\n",
    "                if float(conc) == 0:\n",
    "                    saltfix=';'\n",
    "                else:\n",
    "                    saltfix=''\n",
    "                f.write(topology_file.format(homedir=workdir, saltfix=saltfix,\n",
    "                                             cation_resname=salt_details['cation'].upper(),\n",
    "                                             anion_resname=salt_details['anion'].upper(),\n",
    "                                             N_water=conc_details['N_water'],\n",
    "                                             N_cation=conc_details['N_cation'],\n",
    "                                             N_anion=conc_details['N_anion']\n",
    "                                             ))\n",
    "            f.close()\n",
    "            print('Wrote initial configurations and topology files to'+os.getcwd())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed3c8840",
   "metadata": {},
   "source": [
    "### Molecular dynamics setup using OpenMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1450185f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd -q $homedir\n",
    "\n",
    "openmm_script=\"\"\"# Imports\n",
    "import sys\n",
    "import os\n",
    "import openmm as mm\n",
    "from openmm import app\n",
    "from openmm import unit as u\n",
    "from mdtraj.reporters import XTCReporter\n",
    "\n",
    "print('Loading initial configuration and toplogy')\n",
    "pdb = app.PDBFile('Initial.pdb')\n",
    "top = app.GromacsTopFile('Topology.top', periodicBoxVectors=pdb.topology.getPeriodicBoxVectors())\n",
    "\n",
    "# Creating system\n",
    "print('Creating OpenMM System')\n",
    "system = top.createSystem(nonbondedMethod=app.PME, ewaldErrorTolerance=0.0005, switchDistance=1*u.nanometer,\n",
    "                          nonbondedCutoff=1.2*u.nanometers, constraints=app.HBonds, rigidWater=True, \n",
    "                          removeCMMotion=True)\n",
    "        \n",
    "# Temperature-coupling by leap frog (BAOAB) Langevin integrator (NVT)\n",
    "integrator = mm.LangevinMiddleIntegrator(298.15*u.kelvin, 1.0/u.picoseconds, 2.0*u.femtoseconds)\n",
    "\n",
    "# Pressure-coupling by a Monte Carlo Barostat (NPT)\n",
    "system.addForce(mm.MonteCarloBarostat(1*u.bar, 298.15*u.kelvin, 25))\n",
    "\n",
    "platform = mm.Platform.getPlatformByName('CUDA')\n",
    "properties = {{'CudaPrecision': 'mixed', 'CudaDeviceIndex': '0'}}\n",
    "\n",
    "# Create the Simulation object\n",
    "sim = app.Simulation(top.topology, system, integrator, platform, properties)\n",
    "\n",
    "# Create a dictionary to store all forces\n",
    "forces = {{sim.system.getForce(index).__class__.__name__: sim.system.getForce(index) for index in range(sim.system.getNumForces())}}\n",
    "\n",
    "# Set the particle positions\n",
    "sim.context.setPositions(pdb.positions)\n",
    "\n",
    "# Minimize the energy\n",
    "print('Minimizing energy')\n",
    "sim.minimizeEnergy()\n",
    "    \n",
    "# Draw initial MB velocities\n",
    "sim.context.setVelocitiesToTemperature(298.15*u.kelvin)\n",
    "\n",
    "# Equlibrate simulation\n",
    "print('Equilibrating...')\n",
    "sim.step(500000)  # 500000*2 fs = 1.0 ns\n",
    "\n",
    "# Set up the reporters\n",
    "sim.reporters.append(app.StateDataReporter('output.dat', {outFreq}, totalSteps={Nsteps}+500000,\n",
    "    time=True, potentialEnergy=True, kineticEnergy=True, temperature=True, volume=True, density=True,\n",
    "    remainingTime=True, speed=True, separator='\\t'))\n",
    "\n",
    "# Set up trajectory reporter\n",
    "sim.reporters.append(XTCReporter('trajectory.xtc', reportInterval={outFreq}, append=False))\n",
    "\n",
    "# Run dynamics\n",
    "print('Running dynamics! (NPT)')\n",
    "sim.step({Nsteps})\n",
    "\n",
    "# Print PME information\n",
    "print('''\n",
    "PARTICLE MESH EWALD PARAMETERS\n",
    "Separation parameter: {{}}\n",
    "Number of grid points along the X axis: {{}}\n",
    "Number of grid points along the Y axis: {{}}\n",
    "Number of grid points along the Z axis: {{}}\n",
    "'''.format(*forces['NonbondedForce'].getPMEParametersInContext(sim.context)))\n",
    "\"\"\"           \n",
    "\n",
    "if GENERATE_REFERENCES:\n",
    "    N_simulations = 0\n",
    "    for i, salt in enumerate(salts):\n",
    "        for conc in concentrations:\n",
    "            if conc == 0 and i == 0:\n",
    "                %cd -q $homedir/Simulations/References/No_salt\n",
    "            elif conc == 0 and i > 0:\n",
    "                continue\n",
    "            else:\n",
    "                conc = '{0:.2f}'.format(conc)\n",
    "                %cd -q $homedir/Simulations/References/$salt/$conc\n",
    "        \n",
    "            with open('openMM.py', 'w') as f:\n",
    "                f.write(openmm_script.format(Nsteps=states['ref']['Nsteps'], outFreq=states['ref']['OutFreq']))\n",
    "            f.close()\n",
    "            print('Wrote run_openMM.py files to '+os.getcwd())\n",
    "            N_simulations+=1\n",
    "\n",
    "print('Simulations about to be submitted: {}'.format(N_simulations))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9562c69",
   "metadata": {},
   "source": [
    "### Submit script"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d8e0465",
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd -q $homedir\n",
    "\n",
    "submit_script=\"\"\"#!/bin/bash\n",
    "#PBS -l select=1:ncpus=1:mpiprocs=1:ompthreads=1:jobtype=gpu:ngpus=1\n",
    "#PBS -l walltime=4:00:00\n",
    "#PBS -N Reference\n",
    "#PBS -e run.err\n",
    "#PBS -o run.out\n",
    "\n",
    "source ~/.bashrc\n",
    "source ~/.bash_profile\n",
    "\n",
    "cd {path}\n",
    "\n",
    "python openMM.py\n",
    "\"\"\"\n",
    "\n",
    "master_submit_script=\"\"\"#!/bin/bash\n",
    "read -p \"Will submit {N_simulations} simulations. Do you want to proceed? (yes/no) \" yn\n",
    "\n",
    "case $yn in \n",
    "    yes ) echo submitting...;;\n",
    "    no ) echo exiting...;\n",
    "         exit;;\n",
    "    * ) echo invalid response;\n",
    "        exit 1;;\n",
    "esac\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "if GENERATE_REFERENCES:\n",
    "    with open('Simulations/master_reference.sh', 'w') as ff:\n",
    "        ff.write(master_submit_script.format(N_simulations=N_simulations))\n",
    "        for i, salt in enumerate(salts):\n",
    "            for conc in concentrations:\n",
    "                if conc == 0 and i == 0:\n",
    "                    %cd -q $homedir/Simulations/References/No_salt\n",
    "                elif conc == 0 and i > 0:\n",
    "                    continue\n",
    "                else:\n",
    "                    conc = '{0:.2f}'.format(conc)\n",
    "                    %cd -q $homedir/Simulations/References/$salt/$conc\n",
    "                \n",
    "                simpath = os.getcwd()\n",
    "                simpath = simpath.split(sep='SI-Cavity-Formation')[1]\n",
    "                simpath = workdir+simpath\n",
    "                with open('submit.pbs', 'w') as f:\n",
    "                    f.write(submit_script.format(path=simpath))\n",
    "                f.close()\n",
    "                s = 'cd {}\\njsub submit.pbs\\n\\n'.format(simpath)\n",
    "                ff.write(s)\n",
    "        f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfd4c6af",
   "metadata": {},
   "source": [
    "## Step 3. Generation of Solvent Configurations in Presence of Solute -<br> Solution System (λ = 1)\n",
    "### Construction of initial configurations and topology file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3fc001c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd -q $homedir\n",
    "\n",
    "packmol_script=\"\"\"\n",
    "tolerance 2.0\n",
    "filetype pdb\n",
    "output Initial.pdb\n",
    "\n",
    "add_box_sides 0.1\n",
    "\n",
    "structure {homedir}/PDB_files/{category}/{solute}.pdb\n",
    "        number 1\n",
    "        fixed 33.5 33.5 33.5 0. 0. 0.\n",
    "        centerofmass\n",
    "end structure\n",
    "\n",
    "{saltfix}structure {homedir}/PDB_files/Ions/{cation}.pdb\n",
    "{saltfix}        number {N_cation}\n",
    "{saltfix}        inside cube 0. 0. 0. 67.\n",
    "{saltfix}end structure\n",
    "        \n",
    "{saltfix}structure {homedir}/PDB_files/Ions/{anion}.pdb\n",
    "{saltfix}        number {N_anion}\n",
    "{saltfix}        inside cube 0. 0. 0. 67.\n",
    "{saltfix}end structure\n",
    "\n",
    "structure {homedir}/PDB_files/spce.pdb\n",
    "        number {N_water}\n",
    "        inside cube 0. 0. 0. 67\n",
    "end structure\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "topology_file = \"\"\"\n",
    "; Include forcefield files\n",
    "#include \"{homedir}/Force_fields/ions.itp\"\n",
    "#include \"{homedir}/Force_fields/spce.itp\"\n",
    "#include \"{homedir}/Force_fields/{solute}.itp\"\n",
    "\n",
    "[ defaults ]\n",
    "; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ\n",
    "1               2               no              1            0.83333333\n",
    "\n",
    "[ system ]\n",
    "{solute} in salt water\n",
    "\n",
    "[ molecules ]\n",
    ";molecule name   nr.\n",
    "{solute_resname}          1\n",
    "{saltfix}{cation_resname}          {N_cation}\n",
    "{saltfix}{anion_resname}          {N_anion}\n",
    "SOL          {N_water}\n",
    "\"\"\"\n",
    "\n",
    "if GENERATE_SOLUTIONS:\n",
    "    water_done = False\n",
    "    for solute, solute_details in solutes.items():\n",
    "        for salt, salt_details in salts.items():\n",
    "            for conc, conc_details in concentrations.items():\n",
    "                if conc == 0 and water_done == False:\n",
    "                    %cd -q $homedir/Simulations/$solute/No_salt\n",
    "                    water_done = True\n",
    "                elif conc == 0 and water_done == True:\n",
    "                    continue\n",
    "                else:\n",
    "                    conc = '{0:.2f}'.format(conc)\n",
    "                    %cd -q $homedir/Simulations/$solute/$salt/$conc\n",
    "                \n",
    "                # Construction of initial configuration\n",
    "                with open('packmol.in', 'w') as f:\n",
    "                    # Fix for no salt\n",
    "                    if float(conc) == 0:\n",
    "                        saltfix='#'\n",
    "                    else:\n",
    "                        saltfix=''\n",
    "                    f.write(packmol_script.format(homedir=homedir, saltfix=saltfix,\n",
    "                                                  category=solute_details['Category'], solute=solute, N_water=conc_details['N_water'],\n",
    "                                                  cation=salt_details['cation'], N_cation=conc_details['N_cation'], \n",
    "                                                  anion=salt_details['anion'], N_anion=conc_details['N_anion']))\n",
    "                f.close()\n",
    "                !packmol < packmol.in\n",
    "                \n",
    "                # Construction of topology file\n",
    "                with open('Topology.top', 'w') as f:\n",
    "                    if float(conc) == 0:\n",
    "                        saltfix=';'\n",
    "                    else:\n",
    "                        saltfix=''\n",
    "                    f.write(topology_file.format(homedir=homedir, saltfix=saltfix,\n",
    "                                                 solute=solute, solute_resname=solute_details['Resname'],\n",
    "                                                 cation_resname=salt_details['cation'].upper(),\n",
    "                                                 anion_resname=salt_details['anion'].upper(),\n",
    "                                                 N_water=conc_details['N_water'],\n",
    "                                                 N_cation=conc_details['N_cation'],\n",
    "                                                 N_anion=conc_details['N_anion']\n",
    "                                                 ))\n",
    "                f.close()\n",
    "                \n",
    "                mol = pmd.load_file('Topology.top', xyz='Initial.pdb')\n",
    "                mol.save('Topology_stand_alone.top', overwrite=True)\n",
    "                \n",
    "                print('Wrote initial configurations and topology files to '+os.getcwd())\n",
    "        water_done = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36629028",
   "metadata": {},
   "source": [
    "### Molecular dynamics setup using OpenMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dfc85c8",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%cd -q $homedir\n",
    "\n",
    "openmm_script=\"\"\"# Imports\n",
    "import sys\n",
    "import os\n",
    "import openmm as mm\n",
    "from openmm import app\n",
    "from openmm import unit as u\n",
    "from mdtraj.reporters import XTCReporter\n",
    "\n",
    "print('Loading initial configuration and toplogy')\n",
    "pdb = app.PDBFile('Initial.pdb')\n",
    "top = app.GromacsTopFile('Topology_stand_alone.top', periodicBoxVectors=pdb.topology.getPeriodicBoxVectors())\n",
    "\n",
    "# Creating system\n",
    "print('Creating OpenMM System')\n",
    "system = top.createSystem(nonbondedMethod=app.PME, ewaldErrorTolerance=0.0005, switchDistance=1*u.nanometer,\n",
    "                          nonbondedCutoff=1.2*u.nanometers, constraints=app.HBonds, rigidWater=True, \n",
    "                          removeCMMotion=True)\n",
    "        \n",
    "# Temperature-coupling by leap frog (BAOAB) Langevin integrator (NVT)\n",
    "integrator = mm.LangevinMiddleIntegrator(298.15*u.kelvin, 1.0/u.picoseconds, 2.0*u.femtoseconds)\n",
    "\n",
    "# Pressure-coupling by a Monte Carlo Barostat (NPT)\n",
    "system.addForce(mm.MonteCarloBarostat(1*u.bar, 298.15*u.kelvin, 25))\n",
    "\n",
    "platform = mm.Platform.getPlatformByName('CUDA')\n",
    "properties = {{'CudaPrecision': 'mixed', 'CudaDeviceIndex': '0'}}\n",
    "\n",
    "# Create the Simulation object\n",
    "sim = app.Simulation(top.topology, system, integrator, platform, properties)\n",
    "\n",
    "# Create a dictionary to store all forces\n",
    "forces = {{sim.system.getForce(index).__class__.__name__: sim.system.getForce(index) for index in range(sim.system.getNumForces())}}\n",
    "\n",
    "# Set the particle positions\n",
    "sim.context.setPositions(pdb.positions)\n",
    "\n",
    "# Minimize the energy\n",
    "print('Minimizing energy')\n",
    "sim.minimizeEnergy()\n",
    "    \n",
    "# Draw initial MB velocities\n",
    "sim.context.setVelocitiesToTemperature(298.15*u.kelvin)\n",
    "\n",
    "# Equlibrate simulation\n",
    "print('Equilibrating...')\n",
    "sim.step(500000)  # 500000*2 fs = 1.0 ns\n",
    "\n",
    "# Set up the reporters\n",
    "sim.reporters.append(app.StateDataReporter('output.dat', {outFreq}, totalSteps={Nsteps}+500000,\n",
    "    time=True, potentialEnergy=True, kineticEnergy=True, temperature=True, volume=True, density=True,\n",
    "    remainingTime=True, speed=True, separator='\\t'))\n",
    "\n",
    "# Set up trajectory reporter\n",
    "sim.reporters.append(XTCReporter('trajectory.xtc', reportInterval={outFreq}, append=False))\n",
    "\n",
    "# Run dynamics\n",
    "print('Running dynamics! (NPT)')\n",
    "sim.step({Nsteps})\n",
    "\n",
    "# Print PME information\n",
    "print('''\n",
    "PARTICLE MESH EWALD PARAMETERS\n",
    "Separation parameter: {{}}\n",
    "Number of grid points along the X axis: {{}}\n",
    "Number of grid points along the Y axis: {{}}\n",
    "Number of grid points along the Z axis: {{}}\n",
    "'''.format(*forces['NonbondedForce'].getPMEParametersInContext(sim.context)))\n",
    "\"\"\"           \n",
    "\n",
    "N_simulations = 0\n",
    "if GENERATE_SOLUTIONS:\n",
    "    water_done = False\n",
    "    for solute, solute_details in solutes.items():\n",
    "        for salt, salt_details  in salts.items():\n",
    "            for conc, conc_details in concentrations.items():\n",
    "                if conc == 0 and water_done == False:\n",
    "                    %cd -q $homedir/Simulations/$solute/No_salt\n",
    "                    water_done = True\n",
    "                elif conc == 0 and water_done == True:\n",
    "                    continue\n",
    "                else:\n",
    "                    conc = '{0:.2f}'.format(conc)\n",
    "                    %cd -q $homedir/Simulations/$solute/$salt/$conc\n",
    "        \n",
    "                with open('openMM.py', 'w') as f:\n",
    "                    f.write(openmm_script.format(Nsteps=states['sol']['Nsteps'], outFreq=states['sol']['OutFreq']))\n",
    "                f.close()\n",
    "                print('Wrote run_openMM.py files to '+os.getcwd())\n",
    "                N_simulations+=1\n",
    "        water_done = False\n",
    "\n",
    "print('Simulations about to be submitted: {}'.format(N_simulations))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f398e2a",
   "metadata": {},
   "source": [
    "### Submit script"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e13f71f4",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%cd -q $homedir\n",
    "\n",
    "submit_script=\"\"\"#!/bin/bash\n",
    "#PBS -l select=1:ncpus=1:mpiprocs=1:ompthreads=1:jobtype=gpu:ngpus=1\n",
    "#PBS -l walltime=4:00:00\n",
    "#PBS -N Solution\n",
    "#PBS -e run.err\n",
    "#PBS -o run.out\n",
    "\n",
    "source ~/.bashrc\n",
    "source ~/.bash_profile\n",
    "\n",
    "cd {path}\n",
    "\n",
    "python openMM.py\n",
    "\"\"\"\n",
    "\n",
    "master_submit_script=\"\"\"#!/bin/bash\n",
    "read -p \"Will submit {N_simulations} simulations. Do you want to proceed? (yes/no) \" yn\n",
    "\n",
    "case $yn in \n",
    "    yes ) echo submitting...;;\n",
    "    no ) echo exiting...;\n",
    "         exit;;\n",
    "    * ) echo invalid response;\n",
    "        exit 1;;\n",
    "esac\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "if GENERATE_SOLUTIONS:\n",
    "    with open('Simulations/master_solution.sh', 'w') as ff:\n",
    "        ff.write(master_submit_script.format(N_simulations=N_simulations))\n",
    "        water_done = False\n",
    "        for solute, solute_details in solutes.items():\n",
    "            for salt, salt_details  in salts.items():\n",
    "                for conc, conc_details in concentrations.items():\n",
    "                    if conc == 0 and water_done == False:\n",
    "                        %cd -q $homedir/Simulations/$solute/No_salt\n",
    "                        water_done = True\n",
    "                    elif conc == 0 and water_done == True:\n",
    "                        continue\n",
    "                    else:\n",
    "                        conc = '{0:.2f}'.format(conc)\n",
    "                        %cd -q $homedir/Simulations/$solute/$salt/$conc\n",
    "                        \n",
    "                    simpath = os.getcwd()\n",
    "                    simpath = simpath.split(sep='SI-Cavity-Formation')[1]\n",
    "                    simpath = workdir+simpath\n",
    "            \n",
    "                    with open('submit.pbs', 'w') as f:\n",
    "                        f.write(submit_script.format(path=simpath))\n",
    "                    f.close()\n",
    "                    s = 'cd {}\\njsub submit.pbs\\n\\n'.format(simpath)\n",
    "                    ff.write(s)\n",
    "            water_done = False\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91d01ec1",
   "metadata": {},
   "source": [
    "## EXTRA: RDF Simulations\n",
    "### Conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d90d0fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "salts = {\n",
    "    'KF':   {'cation': 'K',  'anion': 'F' },\n",
    "    'KCl':  {'cation': 'K',  'anion': 'Cl'},\n",
    "    'KI':   {'cation': 'K',  'anion': 'I' },\n",
    "}\n",
    "\n",
    "solutes = {\n",
    "    'n-hexane':   {'SMILES': 'CCCCCC',  'Category': 'n-alkanes',  'Resname': 'HEX'},\n",
    "}\n",
    "\n",
    "concentrations = {\n",
    "     1.00: {'N_solute':1, 'N_water':10000, 'N_cation':180, 'N_anion':180}\n",
    "}\n",
    "\n",
    "Nreplica = 50"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "709d03a5",
   "metadata": {},
   "source": [
    "### Construction of initial configurations and topology file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e34e53a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd -q $homedir\n",
    "\n",
    "packmol_script=\"\"\"\n",
    "tolerance 2.0\n",
    "filetype pdb\n",
    "output Initial.pdb\n",
    "\n",
    "add_box_sides 0.1\n",
    "\n",
    "structure {homedir}/PDB_files/{category}/{solute}.pdb\n",
    "        number 1\n",
    "        fixed 33.5 33.5 33.5 0. 0. 0.\n",
    "        centerofmass\n",
    "end structure\n",
    "\n",
    "structure {homedir}/PDB_files/Ions/{cation}.pdb\n",
    "        number {N_cation}\n",
    "        inside cube 0. 0. 0. 67.\n",
    "end structure\n",
    "\n",
    "structure {homedir}/PDB_files/Ions/{anion}.pdb\n",
    "        number {N_anion}\n",
    "        inside cube 0. 0. 0. 67.\n",
    "end structure\n",
    "\n",
    "structure {homedir}/PDB_files/spce.pdb\n",
    "        number {N_water}\n",
    "        inside cube 0. 0. 0. 67\n",
    "end structure\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "topology_file = \"\"\"\n",
    "; Include forcefield files\n",
    "#include \"{homedir}/Force_fields/ions.itp\"\n",
    "#include \"{homedir}/Force_fields/spce.itp\"\n",
    "#include \"{homedir}/Force_fields/{solute}.itp\"\n",
    "\n",
    "[ defaults ]\n",
    "; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ\n",
    "1               2               no              1            0.83333333\n",
    "\n",
    "[ system ]\n",
    "{solute} in salt water\n",
    "\n",
    "[ molecules ]\n",
    ";molecule name   nr.\n",
    "{solute_resname}          1\n",
    "{cation_resname}          {N_cation}\n",
    "{anion_resname}          {N_anion}\n",
    "SOL          {N_water}\n",
    "\"\"\"\n",
    "\n",
    "if GENERATE_RDF:\n",
    "    for solute, solute_details in solutes.items():\n",
    "        for salt, salt_details in salts.items():\n",
    "            for conc, conc_details in concentrations.items():\n",
    "                conc = '{0:.2f}_rdf'.format(conc)\n",
    "                for replica in range(Nreplica):\n",
    "                    %cd -q $homedir/Simulations/$solute/$salt/$conc/$replica\n",
    "                \n",
    "                    # Construction of initial configuration\n",
    "                    with open('packmol.in', 'w') as f:\n",
    "                        f.write(packmol_script.format(homedir=homedir, \n",
    "                                                      category=solute_details['Category'], solute=solute, N_water=conc_details['N_water'],\n",
    "                                                      cation=salt_details['cation'], N_cation=conc_details['N_cation'], \n",
    "                                                      anion=salt_details['anion'], N_anion=conc_details['N_anion']))\n",
    "                    f.close()\n",
    "                    !packmol < packmol.in\n",
    "                    \n",
    "                    # Construction of topology file\n",
    "                    with open('Topology.top', 'w') as f:\n",
    "                        f.write(topology_file.format(homedir=homedir,\n",
    "                                                     solute=solute, solute_resname=solute_details['Resname'],\n",
    "                                                     cation_resname=salt_details['cation'].upper(),\n",
    "                                                     anion_resname=salt_details['anion'].upper(),\n",
    "                                                     N_water=conc_details['N_water'],\n",
    "                                                     N_cation=conc_details['N_cation'],\n",
    "                                                     N_anion=conc_details['N_anion']\n",
    "                                                     ))\n",
    "                    f.close()\n",
    "                    \n",
    "                    mol = pmd.load_file('Topology.top', xyz='Initial.pdb')\n",
    "                    mol.save('Topology_stand_alone.top', overwrite=True)\n",
    "                    \n",
    "                    print('Wrote initial configurations and topology files to '+os.getcwd())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "089b3dc6",
   "metadata": {},
   "source": [
    "### Molecular dynamics setup using OpenMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4f2a010",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%cd -q $homedir\n",
    "\n",
    "openmm_script=\"\"\"# Imports\n",
    "import sys\n",
    "import os\n",
    "import openmm as mm\n",
    "from openmm import app\n",
    "from openmm import unit as u\n",
    "from mdtraj.reporters import XTCReporter\n",
    "\n",
    "print('Loading initial configuration and toplogy')\n",
    "pdb = app.PDBFile('Initial.pdb')\n",
    "top = app.GromacsTopFile('Topology_stand_alone.top', periodicBoxVectors=pdb.topology.getPeriodicBoxVectors())\n",
    "\n",
    "# Creating system\n",
    "print('Creating OpenMM System')\n",
    "system = top.createSystem(nonbondedMethod=app.PME, ewaldErrorTolerance=0.0005, switchDistance=1*u.nanometer,\n",
    "                          nonbondedCutoff=1.2*u.nanometers, constraints=app.HBonds, rigidWater=True, \n",
    "                          removeCMMotion=True)\n",
    "        \n",
    "# Temperature-coupling by leap frog (BAOAB) Langevin integrator (NVT)\n",
    "integrator = mm.LangevinMiddleIntegrator(298.15*u.kelvin, 1.0/u.picoseconds, 2.0*u.femtoseconds)\n",
    "\n",
    "# Pressure-coupling by a Monte Carlo Barostat (NPT)\n",
    "system.addForce(mm.MonteCarloBarostat(1*u.bar, 298.15*u.kelvin, 25))\n",
    "\n",
    "platform = mm.Platform.getPlatformByName('CUDA')\n",
    "properties = {{'CudaPrecision': 'mixed', 'CudaDeviceIndex': '0'}}\n",
    "\n",
    "# Create the Simulation object\n",
    "sim = app.Simulation(top.topology, system, integrator, platform, properties)\n",
    "\n",
    "# Create a dictionary to store all forces\n",
    "forces = {{sim.system.getForce(index).__class__.__name__: sim.system.getForce(index) for index in range(sim.system.getNumForces())}}\n",
    "\n",
    "# Set the particle positions\n",
    "sim.context.setPositions(pdb.positions)\n",
    "\n",
    "# Minimize the energy\n",
    "print('Minimizing energy')\n",
    "sim.minimizeEnergy()\n",
    "    \n",
    "# Draw initial MB velocities\n",
    "sim.context.setVelocitiesToTemperature(298.15*u.kelvin)\n",
    "\n",
    "# Equlibrate simulation\n",
    "print('Equilibrating...')\n",
    "sim.step(500000)  # 500000*2 fs = 1.0 ns\n",
    "\n",
    "# Set up the reporters\n",
    "sim.reporters.append(app.StateDataReporter('output.dat', {outFreq}, totalSteps={Nsteps}+500000,\n",
    "    time=True, potentialEnergy=True, kineticEnergy=True, temperature=True, volume=True, density=True,\n",
    "    remainingTime=True, speed=True, separator='\\t'))\n",
    "\n",
    "# Set up trajectory reporter\n",
    "sim.reporters.append(XTCReporter('trajectory.xtc', reportInterval={outFreq}, append=False))\n",
    "\n",
    "# Run dynamics\n",
    "print('Running dynamics! (NPT)')\n",
    "sim.step({Nsteps})\n",
    "\n",
    "# Print PME information\n",
    "print('''\n",
    "PARTICLE MESH EWALD PARAMETERS\n",
    "Separation parameter: {{}}\n",
    "Number of grid points along the X axis: {{}}\n",
    "Number of grid points along the Y axis: {{}}\n",
    "Number of grid points along the Z axis: {{}}\n",
    "'''.format(*forces['NonbondedForce'].getPMEParametersInContext(sim.context)))\n",
    "\"\"\"           \n",
    "\n",
    "N_simulations = 0\n",
    "if GENERATE_RDF:\n",
    "    for solute, solute_details in solutes.items():\n",
    "        for salt, salt_details in salts.items():\n",
    "            for conc, conc_details in concentrations.items():\n",
    "                conc = '{0:.2f}_rdf'.format(conc)\n",
    "                for replica in range(Nreplica):\n",
    "                    %cd -q $homedir/Simulations/$solute/$salt/$conc/$replica\n",
    "        \n",
    "                    with open('openMM.py', 'w') as f:\n",
    "                        f.write(openmm_script.format(Nsteps=states['sol']['Nsteps'], outFreq=states['sol']['OutFreq']))\n",
    "                    f.close()\n",
    "                    print('Wrote run_openMM.py files to '+os.getcwd())\n",
    "                    N_simulations+=1\n",
    "\n",
    "print('Simulations about to be submitted: {}'.format(N_simulations))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88d89295",
   "metadata": {},
   "source": [
    "### Submit script"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b1ecdca",
   "metadata": {},
   "outputs": [],
   "source": [
    "submit_script=\"\"\"#!/bin/bash\n",
    "#PBS -q F1accs\n",
    "#PBS -N MD_RDF\n",
    "#PBS -l select=1:ncpus=16:mpiprocs=16:ompthreads=1\n",
    "#PBS -l walltime=2:00:00\n",
    "#PBS -e run.err\n",
    "#PBS -o run.out\n",
    "\n",
    "source ~/.bashrc\n",
    "source ~/.bash_profile\n",
    "\n",
    "module load nvhpc-nompi/22.2_cuda11.6\n",
    "\n",
    "cd {path}\n",
    "\n",
    "python openMM.py\n",
    "\"\"\"\n",
    "\n",
    "master_submit_script=\"\"\"#!/bin/bash\n",
    "read -p \"Will submit {N_simulations} MD simulations. Do you want to proceed? (yes/no) \" yn\n",
    "\n",
    "case $yn in \n",
    "    yes ) echo submitting...;;\n",
    "    no ) echo exiting...;\n",
    "         exit;;\n",
    "    * ) echo invalid response;\n",
    "        exit 1;;\n",
    "esac\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "%cd -q $homedir\n",
    "workdir = '/home/m0053/m005302/work/SI-Cavity-Formation'\n",
    "\n",
    "if GENERATE_RDF:\n",
    "    with open('Simulations/master_RDF.sh', 'w') as ff:\n",
    "        ff.write(master_submit_script.format(N_simulations=N_simulations))\n",
    "        for solute, solute_details in solutes.items():\n",
    "            for salt, salt_details in salts.items():\n",
    "                for conc, conc_details in concentrations.items():\n",
    "                    conc = '{0:.2f}_rdf'.format(conc)\n",
    "                    for replica in range(Nreplica):\n",
    "                        %cd -q $homedir/Simulations/$solute/$salt/$conc/$replica\n",
    "                \n",
    "                        simpath = os.getcwd()\n",
    "                        simpath = simpath.split(sep='SI-Cavity-Formation')[1]\n",
    "                        simpath = workdir+simpath\n",
    "                        with open('submit_md.pbs', 'w') as f:\n",
    "                            f.write(submit_script.format(path=simpath))\n",
    "                        f.close()\n",
    "                        s = 'cd {}\\nqsub submit_md.pbs\\n\\n'.format(simpath)\n",
    "                        ff.write(s)\n",
    "    ff.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89c30d03",
   "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.9"
  },
  "vscode": {
   "interpreter": {
    "hash": "00291345d083f1793ec18cf48e40fc20d9b764d462ab337f113eb5e44a23dbea"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}