force_constant_hps.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Code used to estimate the force constant to use for harmonic potentials given the number of ions in simulation and the desired maximum concentration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import sys\n",
    "import mpmath\n",
    "sys.modules['sympy.mpmath'] = mpmath\n",
    "from openmm.unit import AVOGADRO_CONSTANT_NA, BOLTZMANN_CONSTANT_kB\n",
    "from openmm.unit import kelvin, bar, litre, kilojoule_per_mole, mole, nanometer, angstrom, kilocalorie_per_mole,molar,molal\n",
    "from openmm.unit import Quantity, Unit\n",
    "import sympy as sp\n",
    "from sympy import symbols, Function, log, Sum\n",
    "from sympy.abc import x, y, z\n",
    "import scipy.integrate\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.special import erf\n",
    "from scipy.optimize import root"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constants and variables to modify according to simulation setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = BOLTZMANN_CONSTANT_kB * AVOGADRO_CONSTANT_NA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# System composition\n",
    "Cmax_target = 3.26 * mole* litre**-1 # maximum molal concentration desired in units of molar\n",
    "N_ions = 217 #number of ions corresponding to desired cmax\n",
    "\n",
    "# Simulation box\n",
    "L_z = 14.4 * nanometer # total length of box in z direction\n",
    "Lx, Ly, Lz = np.array([4.8, 4.8, 14.4]) * nanometer\n",
    "\n",
    "# Simulation temperature\n",
    "T = 300 * kelvin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical Approximation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_C_max_from_spring_const(k : Quantity) -> Quantity:\n",
    "    A = Lx * Ly\n",
    "    K = k / (2 * R * T)\n",
    "    Ksq = np.sqrt(K._value) * (1*K.unit).sqrt()\n",
    "\n",
    "    Cmax = (N_ions/A * Ksq) / (np.sqrt(np.pi) * erf(Lz / (2 * Ksq)))\n",
    "    Cmax = (Cmax / AVOGADRO_CONSTANT_NA).in_units_of(molar)\n",
    "\n",
    "    return Cmax\n",
    "\n",
    "def Cmax_for_roots(k : float) -> float:\n",
    "    return calc_C_max_from_spring_const(\n",
    "        k * kilojoule_per_mole/nanometer**2\n",
    "        ) - Cmax_target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated force constant:  [0.68095403] \n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/bamo6610/miniconda3/envs/calculus/lib/python3.10/site-packages/openmm/unit/quantity.py:414: RuntimeWarning: divide by zero encountered in reciprocal\n",
      "  return self * pow(other, -1.0)\n"
     ]
    }
   ],
   "source": [
    "#Estimate force constant\n",
    "soln = root(lambda k : Cmax_for_roots(k)._value, np.float64(0.0))\n",
    "print('Estimated force constant: ',soln.x,\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cmax from spring constant:  [3.26] M\n"
     ]
    }
   ],
   "source": [
    "# Check force constant returns desired maximum concentration\n",
    "Cmax = calc_C_max_from_spring_const(\n",
    "    k=soln.x * kilojoule_per_mole/nanometer**2\n",
    "    )\n",
    "print('Cmax from spring constant: ',Cmax)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "calculus",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}