{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Energy versus Temperature

\n", "

A Molecular Dynamics simulation of 16 particles in a box. The particles interact via the Lennard Jones potential and thus the force:\n", "$$\n", "\\mathbf{f} (r) = \\frac{48}{r^2}\n", "\\left[\\left( \\frac{1}{r}\\right)^{12} - \\frac{1}{2}\n", "\\left(\\frac{1}{r} \\right)^6 \\right] \\mathbf{r}.\n", "$$\n", "The positions and velocities of the particles obey the Verlet-velocity algorithm:\n", "\\begin{eqnarray}\n", "\\mathbf{r}_i(t+h) & \\simeq & \\mathbf{r}_i(t) + h \\mathbf{v}_i(t) + \\frac{h^2}\n", "{2} \\mathbf{F}_i(t) + \\mbox{O}(h^3), \\label{VV1}\\\\[6pt]\n", "\\mathbf{v}_i(t+h) & \\simeq & \\mathbf{v}_i(t) + h \\,\n", "\\overline{\\mathbf{a}(t)}+ \\mbox{O}(h^2) \\\\[6pt]\n", "& \\simeq &\\mathbf{v}_i(t) + h \\left[ \\frac{\\mathbf{F}_i(t+h)+\\mathbf{F}_i(t)}{2}\n", "\\right] + \\mbox{O}(h^2). \\label{VV2}\n", " \\end{eqnarray}\n", " The specific capacity at constant volume is:\n", "$$\n", "C_V= \\frac{\\partial E}{\\partial T}\\Bigr )_V\n", "$$\n", "The program computes the energy versus temperature and $C_V$.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/html": [ "

" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") { window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")};}else{ element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/glow.min\");}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/glowcomm\");}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/jquery-ui.custom.min\");}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/glow.min\"], function(){console.log(\"GLOW LOADED\");});}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/glowcomm\"], function(){console.log(\"GLOWCOMM LOADED\");});}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/jquery-ui.custom.min\"], function(){console.log(\"JQUERY LOADED\");});}else{element.textContent = ' ';}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "finished\n" ] } ], "source": [ "from vpython import *\n", "import numpy as np\n", "import random\n", "L=1 # side square\n", "time=5\n", "#scene = display(width=400,height=400,range=(1.3) )\n", "scenetemp = graph(xmin=0,xmax=80.0,\n", " width=400, height=400, xtitle='Temperature', ytitle='Total Energy',\n", " ymax=50,ymin=0, title=\"Total Energy vs Temperature\" )\n", "plotT=gcurve(color=color.blue)\n", "sceneEnergy=graph(x=400,xmin=0,xmax=80,ymax=1,ymin=0,\n", " width=400, height=400, xtitle='Temperature', ytitle='Cv',\n", " title=\"Heat capacity at constant Vol, vs T\")\n", "plotCv=gcurve(color=color.cyan)\n", "plotPE=gcurve(color=color.magenta)\n", "plotTE=gcurve(color=color.cyan)\n", "\n", "Natom=16 # number of atoms\n", "#Nr=0 # number particles right side\n", "dt=1e-6 # time step\n", "positions=[] # position of atoms\n", "dN=[] # will contain atoms in right half at each tieme interval\n", "fr=[0]*(Natom) # atoms (spheres)\n", "fr2=[0]*(Natom) # second force\n", "ar=0.03 # radius of atom # a reference velocity\n", "h=0.02\n", "factor=10**(-11.0) # for lennard jones\n", "deltaN=1 # for histogram\n", "Tem=[0]*(75)\n", "KE=[0]*(75)\n", "PtE=[0]*(75)\n", "\n", "def initcond(pref): # initial conditions \n", " Nr=0\n", " vel=[]\n", " global pos\n", " for i in range (0,Natom): # initial positions and velocities\n", " col=vec(1.3*random.random(),1.3*random.random(),1.3*random.random())\n", " x=2.*(L-ar)*random.random()-L+ar #positons of atoms\n", " y=2.*(L-ar)*random.random()-L+ar # in the border forbidden\n", " theta=2*pi*random.random() # select angle 0<=theta<= 2pi\n", " vx=pref*cos(theta) # x component velocity\n", " vy=pref*sin(theta)\n", " positions.append((x,y,0)) # add positions to list\n", " vel.append((vx,vy,0)) # add momentum to list\n", " posi=np.array(positions) # array with positions\n", " ddp=posi[i]\n", " if ddp[0]>=0 and ddp[0]<=L: # count initial atoms at right half\n", " Nr+=1 \n", " v=np.array(vel) # array of velocities\n", " return posi,v\n", "\n", "def sign(a, b): # sign function\n", " if (b >= 0.0):\n", " return abs(a)\n", " else:\n", " return - abs(a)\n", "def forces(fr): \n", " fr=[0]*(Natom)\n", " #PE=[0]*(Natom)\n", " PE=0\n", " for i in range( 0, Natom-1 ):\n", " for j in range(i + 1, Natom):\n", " dr=posi[i]-posi[j] # relative position between particles\n", " if (abs(dr[0]) > L): # smallest distance from part.or image\n", " dr[0] = dr[0] - sign(2*L, dr[0]) # interact with closer image\n", " if (abs(dr[1]) > L): # same for y\n", " dr[1] = dr[1] - sign(2*L, dr[1])\n", " #r2=mag2(dr)\n", " r2=dr[0]**2+dr[1]**2+dr[2]**2\n", " if r2 < 0.02: # to avoid 0 denominator\n", " r2 = 0.03 # minimum distance between 2 atoms\n", " invr2 = 1./r2 # compute this factor\n", " fij =invr2*factor* 48.*(invr2**3 - 0.5) *invr2**3 #\n", " fr[i]=fij*dr+ fr[i]\n", " fr[j]=-fij*dr +fr[j] # lennard jones force\n", " PE = PE + factor* 4.*(invr2**3)*((invr2**3) - 1.)\n", " return fr,PE \n", "nn=0 # counter\n", "\n", "for velo in arange(0.5,20,1.0):\n", " #count=0\n", " pref=velo\n", " posi,v=initcond(pref)\n", " for t in range (0,time): # time steps to make averages\n", " Nr=0 # begin at zero in each time\n", " v2=0 \n", " PE=0\n", " for i in range(0,Natom):\n", " fr,potE=forces(fr)\n", " rate(500)\n", " dpos=posi[i]\n", " if dpos[0] <= -L:\n", " posi[i] = [dpos[0]+2*L,dpos[1],0] # x periodic boundary conditions\n", " if dpos[0] >= L:\n", " posi[i] = [dpos[0]-2*L,dpos[1],0]\n", " if dpos[1] <= -L:\n", " posi[i] = [dpos[0],dpos[1]+2*L,0] # y periodic boundary conditions\n", " if dpos[1] >= L:\n", " posi[i] = [dpos[0],dpos[1]-2*L,0]\n", " dpos=posi[i]\n", " if dpos[0]>0 and dpos[0]1: # forward difference derivative to find: Cv=dE/dT\n", " plotCv.plot(pos=(Tem[nn],(PtE[nn]+KE[nn]-PtE[nn-1]-KE[nn-1])/(Tem[nn]-Tem[nn-1])))\n", " nn=nn+1\n", "print(\"finished\")\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.13" } }, "nbformat": 4, "nbformat_minor": 2 }