Updated code NAMD ion-flux demo to new version. Added run.sh files to run the NAMD simulation, but NAMD needs to be built to execute the command

Author: hcho38

namd-demos/exec-files/run.sh

@@ -0,0 +1,4 @@
+#!/bin/bash
+
+~/namd-3.0/Linux-x86_64-g++/namd3 +p8 150mV-2.conf
+

namd-demos/ion-flux/ion-flux.ipynb

@@ -2,12 +2,11 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Import required libraries\n",
-    "from imdclient.IMD import IMDReader\n",
+    "#Import required libraries\n",
     "import MDAnalysis as mda\n",
     "\n",
     "import numpy as np\n",
@@ -59,70 +58,69 @@
     }
    ],
    "source": [
-    "%matplotlib widget\n",
+    "#%matplotlib widget\n",
     "NAMD_TOPOL = \"../exec-files/T3_MNN.psf\"\n",
     "\n",
-    "# Initialize MDAnalysis Universe\n",
-    "u = mda.Universe(NAMD_TOPOL, \"imd://localhost:1027\")\n",
-    "\n",
-    "# Enable interactive plotting\n",
     "plt.ion()\n",
     "\n",
-    "# Initialize the plot\n",
     "fig, ax = plt.subplots(figsize=(10,6))\n",
-    "line, = ax.plot([], [], lw=2)  # Line for dynamic updates\n",
-    "ax.set_xlim(0, 1000)  # Adjust as necessary\n",
-    "ax.set_ylim(-150, 150)  # Adjust based on expected distance range\n",
+    "line, = ax.plot([], [], lw=2)\n",
+    "ax.set_xlim(0, 1000)\n",
+    "ax.set_ylim(-150, 150)\n",
     "ax.set_xlabel(\"Simulation Time (steps)\")\n",
     "ax.set_ylabel(\"Charge Flow (nA)\")\n",
     "ax.set_title(\"Ion Current vs. Time\")\n",
     "\n",
-    "# Start the processing thread\n",
     "x_data = np.zeros(144000) \n",
     "y_data = np.zeros(144000)\n",
     "\n",
     "i = 0\n",
     "\n",
-    "timestep = 2.0 # timestep passed to NAMD in conf file\n",
-    "framediff = 1 # IMDfreq\n",
-    "skip = 1 # if any frames need to be skipped\n",
+    "timestep = 2.0 \n",
+    "framediff = 1\n",
+    "skip = 1\n",
     "\n",
     "update_frequency=1\n",
-    "chargeFactor = 1.6022e-19/1.0e-9/1.0e-9   #It's in nanoamps\n",
+    "chargeFactor = 1.6022e-19/1.0e-9/1.0e-9\n",
     "dt = 1.0*timestep*framediff*skip/1.0e6\n",
-    "#ionText = \"name POT CLA CAL SOD LIT CES IOD BRO H11 H21 H31 H41\"\n",
-    "ionText = \"name POT\"\n",
-    "ions = u.select_atoms(ionText)\n",
-    "charge = ions.charges\n",
-    "dims = ions.dimensions\n",
-    "zdist = dims[2]\n",
+    "ionText = \"name POT CLA\"\n",
     "\n",
-    "for t in u.trajectory:\n",
-    "    newPos = ions.positions\n",
-    "    idx = i % 1000\n",
+    "try:\n",
+    "    u = mda.Universe(NAMD_TOPOL, \"imd://localhost:8888\", buffersize = 1000*1024*1024)\n",
+    "    ions = u.select_atoms(ionText)\n",
+    "    charge = ions.charges\n",
+    "    dims = ions.dimensions\n",
+    "    zdist = dims[2]\n",
+    "except Exception as e:\n",
+    "    print(f\"Error during connection: {e}\")\n",
+    "else:\n",
+    "    try:\n",
+    "        for t in u.trajectory:\n",
+    "            newPos = ions.positions\n",
+    "            idx = i % 1000\n",
     "\n",
-    "    if i == 0:\n",
-    "        ionPos = newPos\n",
-    "    else:\n",
-    "        deltas = newPos[:,2]-ionPos[:,2]\n",
-    "        deltas[deltas>zdist/2.] = deltas[deltas>zdist/2.] - zdist\n",
-    "        deltas[deltas<-zdist/2.] = deltas[deltas<-zdist/2.] + zdist\n",
-    "        cd = deltas*charge\n",
-    "        tCur = np.sum(cd)/zdist/dt*chargeFactor\n",
-    "        ionPos = newPos\n",
-    "        time = round(dt * i - dt/2, 4)\n",
-    "        x_data[i] = i + 1\n",
-    "        y_data[i] = tCur\n",
-    "    i = i + 1\n",
+    "            if i == 0:\n",
+    "                ionPos = newPos\n",
+    "            else:\n",
+    "                deltas = newPos[:,2]-ionPos[:,2]\n",
+    "                deltas[deltas>zdist/2.] = deltas[deltas>zdist/2.] - zdist\n",
+    "                deltas[deltas<-zdist/2.] = deltas[deltas<-zdist/2.] + zdist\n",
+    "                cd = deltas*charge\n",
+    "                tCur = np.sum(cd)/zdist/dt*chargeFactor\n",
+    "                ionPos = newPos\n",
+    "                time = round(dt * i - dt/2, 4)\n",
+    "                x_data[i] = i + 1\n",
+    "                y_data[i] = tCur\n",
+    "            i = i + 1\n",
     "\n",
-    "    valid_indices = range(0, idx + 1)\n",
+    "            valid_indices = range(0, idx + 1)\n",
     "\n",
-    "    # Update the plot periodically\n",
-    "    if i % update_frequency == 0:\n",
-    "        line.set_data(x_data[valid_indices], y_data[valid_indices])  # Efficiently update the line with all data\n",
-    "        # ax.set_xlim(0, max(200,max(x_data)+100))  # Adjust X-axis dynamically\n",
-    "        clear_output(wait=True)\n",
-    "        display(fig)"
+    "            if i % update_frequency == 0:\n",
+    "                line.set_data(x_data[valid_indices], y_data[valid_indices])\n",
+    "                clear_output(wait=True)\n",
+    "                display(fig)\n",
+    "    finally:\n",
+    "        u.trajectory.close()"
    ]
   }
  ],