Upload v0.5

CHEG231MD.py

@@ -0,0 +1,179 @@
+"""
+CHEG231MD.py
+Simple MD code in Python (it's slow!) for a Lenard-Jones fluid
+based on QuickBASIC code by Richard L. Rowley [1]
+and FORTRAN code by Smit and Frenkel [2]
+
+Eric M. Furst
+November 2023
+
+[1] Richard L. Rowley (1994). Statistical Mechanics for Thermophysical Property Calculations. Prentice Hall, New York.
+[2] Daan Frenkel and Berend Smit (2002). Understanding Molecular Simulation, 2nd ed. Academic Press, New York.
+
+Revisions
+October 2024 - sped up with numba jit of accel routine and LJ calc.
+
+Requires:
+    python (tested with ver. 3.12.4 build h99e199e_1 channel defaults)
+    numba  (tested with ver. 0.60.0 build py312hd77ebd4_0 channel defaults)
+    numpy  (tested with ver. 1.26.4 build py312h7f4fdc5_0 channel defaults)
+
+Needs:
+    nn - total number of particles per dimension (thus, N = nn^3)
+    rho - dimensionless number density of particles (units of rho/sigma^3)
+    vmax - maximum velocity (units of sqrt[epsilon/m])
+
+Next to do:
+    write a 2D animated version
+    a challenge for students: dimensionless variables
+    can a fast r_cutoff be implemented based on 
+    from scipy.spatial.distance import pdist, squareform?
+"""
+
+import numpy as np
+from numba import jit
+
+X, Y, Z = 0, 1, 2 # helpful for indexing
+
+class MDSimulation:
+    """
+    x - new position coordinate array (N x 3)
+    v - velocity array (N x 3)
+    """
+
+    def __init__(self, nn, rho, vmax):
+        """
+        Initialize the simulation
+        - particles start on a cubic lattice
+        - velocities are assigned a random fraction of the maximum
+        - velocities are corrected to eliminate net flow
+        """
+    
+        # initial variables
+        self.N = nn**3  # total number of particles
+        self.rho = rho
+        self.vol = self.N/rho # volume of the simulation 
+        self.L = self.vol**(1/3) # size of the simulation
+        self.dL = self.L/nn # initial cubic array lattice size
+
+        self.dt = 0.005 # dimensionless time step (0.005 is common)
+        self.dt2 = self.dt*self.dt
+        
+        # initialize position and acceleration arrays 
+        self.x = np.zeros((self.N,3)) # init N by 3 array of positions
+        self.xnew = np.zeros((self.N,3)) # init N by 3 array of positions
+        self.a = np.zeros((self.N,3)) # init N by 3 array of accelerations
+
+        # state of the simulation
+        self.pe = 0 # total potential energy
+        self.ke = 0 # total kinetic energy
+        self.vir = 0 # virial coef
+        self.T = 1 # temperature
+        self.P = 0 # pressure
+
+        # create a cubic array of particles
+        particle = 0
+        for i in range(0, nn):
+            for j in range(0,nn):
+                for k in range(0,nn):
+                    self.x[particle,X]=(i+0.5)*self.dL
+                    self.x[particle,Y]=(j+0.5)*self.dL
+                    self.x[particle,Z]=(k+0.5)*self.dL
+                    particle += 1
+        
+        # assign random velocities from uniform distribution (units?)
+        self.v = vmax*np.random.rand(self.N,3) # init N x 3 array of velocities
+
+        # normlize the velocities s/t net velocity is zero
+        vcm = np.sum(self.v,axis=0)/self.N # 1x3 array of <vx>, <vy>, <vz>    
+        for particle in range(0,self.N):
+            self.v[particle,:] -= vcm
+ 
+        # set the initial accelerations
+        self.accel()
+
+    def move(self):
+        """
+        Verlet velocity algorithm
+        """
+        # find new positions
+        self.xnew = self.x + self.v*self.dt + self.a*self.dt2/2.
+
+        # apply periodic boundary conditions
+        beyond_L = self.xnew[:,:] > self.L # N x 3 Boolean array
+        self.xnew[beyond_L] -= self.L
+        beyond_zero = self.xnew[:,:] < 0
+        self.xnew[beyond_zero] += self.L
+
+        self.x = self.xnew   # update positions
+        self.v = self.v + self.a*self.dt/2 #half update velocity
+
+        # call accelerate 
+        self.accel()
+
+        # finish velocity update
+        self.v = self.v + self.a*self.dt/2 
+
+        # update properties
+        self.ke = 0.5*np.sum(self.v*self.v)
+        self.T = 2*self.ke/self.N/3
+        self.P = self.T*self.rho + self.rho/self.N*self.vir/3
+
+
+    # Non-JIT-compiled class method
+    def accel(self):
+        """                                                                     
+        Acceleration calculated using F = ma                                    
+        """
+        self.a = np.zeros((self.N, 3))  # Zero out all accelerations
+        self.pe = 0
+        self.vir = 0
+
+        # Call the JIT-compiled function to calculate acceleration
+        self.pe, self.vir = accel_jit(self.N, self.x, self.a, self.L)
+
+# Functions are below
+# We're using numba jit to make the loop of the calculation faster
+
+# JIT-compiled function to calculate acceleration
+# Note: numba works by "pass in reference" with numpy arrays,
+# so self.a values are modified and used in self.move()
+# in addition to returning the values of potential energy and
+# the virial factor
+@jit(nopython=True, parallel=False)
+def accel_jit(N, x, a, L):
+    pe = 0
+    vir = 0
+
+    for i in range(0, N - 1):
+        for j in range(i + 1, N):
+            r2 = 0
+            dx = x[i, :] - x[j, :]  # 1x3 array
+            # Apply minimum image convention over periodic boundary conditions
+            for k in range(3):  # Iterate over X, Y, Z
+                if abs(dx[k]) > 0.5 * L:
+                    dx[k] -= np.sign(dx[k]) * L
+                r2 += dx[k] * dx[k]
+
+            # Calculate the force, potential, and virial contribution
+            forceri, potential, virij = LJ(r2)
+            a[i, :] += forceri * dx
+            a[j, :] -= forceri * dx
+            pe += potential
+            vir += virij
+
+    return pe, vir
+
+# JIT-compiled function to calculate interaction
+@jit(nopython=True)
+def LJ(r2):
+    """
+    Compute force and potential from LJ interaction 
+    r2 is r^2, forceri returns force/r
+    """
+    r2i = 1/r2
+    r6i = r2i**3
+    potential = 4*(r6i*r6i-r6i)
+    virij = 48*(r6i*r6i-0.5*r6i)
+    forceri = virij*r2i
+    return forceri, potential, virij

MDsim.ipynb

@@ -0,0 +1,150 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "6d6ace33-bd89-49ac-a148-06c99d7d071c",
+   "metadata": {},
+   "source": [
+    "# Getting Started with CHEG231MD\n",
+    "\n",
+    "Eric M. Furst\\\n",
+    "October 2024\n",
+    "\n",
+    "Test run of the simple MD simulation.\n",
+    "\n",
+    "The temperature will depend on density. At any density, change temp by \n",
+    "scaling max velocity by sqrt(Tnew/Told) \n",
+    "\n",
+    "$$ \\langle v^2 \\rangle = v_\\mathrm{max}^2/2 $$\n",
+    "\n",
+    "Without a thermostat, both the temperature and pressure equilibrate.\n",
+    "This makes calculating an isotherm difficult, but it shows the\n",
+    "combined contributions of potential and kinetic energy\n",
+    "\n",
+    "rho and vmax values that give results close to T+ = 2:\\\n",
+    "0.7 5.8\\\n",
+    "0.6 5.47\\\n",
+    "0.4 5.04"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c7f0016e-5ba0-45b9-9d0e-df7bde2c901f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import CHEG231MD as md\n",
+    "\n",
+    "# create a simulation with # particles, density, and max speed\n",
+    "# N^(1/3), rho+, max vel+\n",
+    "\n",
+    "#sim = md.MDSimulation(8,0.4,5.04)\n",
+    "#sim = md.MDSimulation(8,0.6,5.47)\n",
+    "sim = md.MDSimulation(8,0.7,5.8)\n",
+    "\n",
+    "# initialize variables and lists\n",
+    "pres=[]\n",
+    "temp=[]\n",
+    "time=[]\n",
+    "Pavg = 0\n",
+    "Tavg = 0"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a14e5dc1-7601-4972-911e-0f42eb32a6c3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Run time steps of the simulation\n",
+    "for timestep in range(1,1001):\n",
+    "\n",
+    "    # advance simulation one timestep\n",
+    "    sim.move()\n",
+    "    \n",
+    "    # only average after an equilibration time\n",
+    "    if timestep > 100:\n",
+    "        # production stage\n",
+    "        Pavg += (sim.P-Pavg)/(timestep+1)\n",
+    "        Tavg += (sim.T-Tavg)/(timestep+1)\n",
+    "    else:\n",
+    "        # equilibration stage\n",
+    "        Pavg = sim.P\n",
+    "        Tavg = sim.T\n",
+    "\n",
+    "   # print timestep, ke, pe, total e, Tavg, Pavg        \n",
+    "    if timestep%50==0: # only print every 50 timesteps\n",
+    "        # timestep, <ke>, <pe>, <e>, T, P, Tavg, Pavg\n",
+    "        print(\"%4d  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\" % \n",
+    "              (timestep,sim.ke/sim.N,sim.pe/sim.N,\n",
+    "                (sim.ke+sim.pe)/sim.N,sim.T,sim.P,Tavg,Pavg))\n",
+    "\n",
+    "    pres.append(sim.P)\n",
+    "    temp.append(sim.T)\n",
+    "    time.append(timestep)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d3993efe-4dc7-4b09-9630-90201f989487",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# A different way to calculate the average pressure and temperature\n",
+    "# with standard deviations\n",
+    "\n",
+    "import numpy as np\n",
+    "\n",
+    "print(\"Pressure:    {:.2f}+/-{:.2f}\"\n",
+    "      .format(np.mean(pres[200:]),np.std(pres[200:])))\n",
+    "print(\"Temperature: {:.2f}+/-{:.2f}\"\n",
+    "      .format(np.mean(temp[200:]),np.std(temp[200:])))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "77d0de84-eb9b-4ea5-b5d3-4780335f3d7c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot pressure and temperature\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(time[0:],pres[0:],label=\"P*\")\n",
+    "ax.plot(time[0:],temp[0:],label=\"T*\")\n",
+    "ax.set_xlabel(\"time step\")\n",
+    "ax.legend(frameon=False)\n",
+    "ax.set_xscale('linear')\n",
+    "plt.title(\"Dimensionless instantaneous pressure and temperature\")\n",
+    "\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "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.12.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}