Argon_melting.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "data = np.loadtxt('parameters.txt')\n",
    "\n",
    "def potentials(r):  \n",
    "    rij_intermediate = r.reshape((N, 1, 3)) - r.reshape((1, N, 3))\n",
    "    rij = np.sqrt((rij_intermediate**2).sum(axis = 2))\n",
    "    rij[np.arange(N), np.arange(N)] = np.inf\n",
    "\n",
    "    ri = np.sqrt((r**2).sum(axis=1))\n",
    "\n",
    "    Vs = np.zeros_like(ri)\n",
    "    conditionVs = ri > L\n",
    "\n",
    "    \n",
    "    Vs[conditionVs] = 0.5 * f * (L - ri[conditionVs])**2\n",
    "\n",
    "    auxiliary_to_6pow = (R/rij)**6\n",
    "    auxiliary_to_12pow = auxiliary_to_6pow**2\n",
    "    \n",
    "    Vp = e *(auxiliary_to_12pow - 2 * auxiliary_to_6pow)\n",
    "    Fp = 12*e*(auxiliary_to_12pow - auxiliary_to_6pow).reshape((N,N,1)) * rij_intermediate / rij.reshape((N,N,1))**2\n",
    "    Fp_summed = Fp.sum(axis=1)\n",
    "    Epot = Vs.sum()+Vp.sum() * 0.5\n",
    "\n",
    "    Fs = np.zeros_like(Fp_summed)\n",
    "    Fs[conditionVs] = f * (L - ri[conditionVs, np.newaxis]) * r[conditionVs] / ri[conditionVs, np.newaxis]\n",
    "\n",
    "    F = Fp_summed + Fs\n",
    "\n",
    "    mod_F = np.sqrt((Fs**2).sum(axis=1))\n",
    "    pres = 1/(4*np.pi*L**2)*mod_F.sum()\n",
    "    return F, Epot, pres\n",
    "\n",
    "n = data[0]\n",
    "a = data[6]\n",
    "m = int(data[1])\n",
    "T = data[7] \n",
    "e = data[2]\n",
    "R = data[3] \n",
    "f = int(data[4])\n",
    "L = data[5] \n",
    "tau = data[8]\n",
    "S_d = int(data[10])\n",
    "S_o = int(data[9])\n",
    "S_out = int(data[11])\n",
    "S_xyz = int(data[12])\n",
    "\n",
    "n = int(n)\n",
    "N = n*n*n\n",
    "k = 8.31e-3\n",
    "i = 0\n",
    "r = np.zeros((N, 3))\n",
    "p = np.zeros((N, 3))\n",
    "\n",
    "b0 = np.array([a, 0, 0])\n",
    "b1 = np.array([0.5*a, 0.5*a*np.sqrt(3), 0])\n",
    "b2 = np.array([0.5*a, a*np.sqrt(3)/6, a*np.sqrt(2/3)])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T_av 597.507758544 P_av 13.3580796424 mean_Energy 886.053039319\n"
     ]
    }
   ],
   "source": [
    "T = 1000 #set temperature in absolute units - if 0, crystal will only show gentle oscilations\n",
    "for i0 in range (n):\n",
    "    for i1 in range (n):\n",
    "        for i2 in range (n):\n",
    "            i = i0 + n*i1+n*n*i2\n",
    "            r[i,:] = (i0 - 0.5*(n-1))*b0+(i1 - 0.5*(n-1))*b1+(i2 - 0.5*(n-1))*b2\n",
    "\n",
    "energia = 1\n",
    "with open('output.xyz', 'w') as f_file: # File needs to be saved this way in order to open the output file in Jmoll\n",
    "    f_file.write(\"{}\\n\".format(N))      # This allows user to animate boiling of the crystal in this applet\n",
    "    f_file.write(\"Energy =     {} kcal/mol\\n\".format(energia))\n",
    "    \n",
    "    df = pd.DataFrame(data=r, columns=list('xyz'), index=n*n*n*['Ar'])\n",
    "    df.to_csv(f_file, header=False, sep=' ', float_format='%.5f')\n",
    "\n",
    "    energy = -0.5*k*T*np.log(np.random.random((N, 3)))\n",
    "    E_av_x = np.mean(energy,axis=0)\n",
    "    if (E_av_x != 0).all():\n",
    "        energy *=  (0.5*k*T)/E_av_x\n",
    "    p = np.random.choice([-1, 1], size=(N,3))*np.sqrt(2*m*energy)\n",
    "\n",
    "    p_cm = p.sum(axis=0)\n",
    "    p -= p_cm/N\n",
    "\n",
    "    F, Epot, pres = potentials(r)\n",
    "\n",
    "    mean_temperature = 0\n",
    "    mean_pressure = 0\n",
    "    mean_Energy = 0\n",
    "    s = 0\n",
    "    \n",
    "with open('output.xyz', 'w') as f_file:\n",
    "    for i in range(int(S_d+S_o)):\n",
    "        s+=1\n",
    "        p = p+0.5*F*tau\n",
    "        r = r+p*tau/m\n",
    "        F, Epot, pres = potentials(r)\n",
    "        p = p+0.5*F*tau\n",
    "        Ekin = (p**2).sum(axis=1)/(2*m)\n",
    "        Ekin_total = Ekin.sum()\n",
    "        E = Epot + Ekin_total\n",
    "        T = 2/(3*N*k)*Ekin_total\n",
    "        t = i*tau\n",
    "\n",
    "        if i % S_xyz==0: \n",
    "            \n",
    "            f_file.write(\"{}\\n\".format(N))\n",
    "            f_file.write(\"Energy =     {} kcal/mol\\n\".format(energia))\n",
    "\n",
    "            df = pd.DataFrame(data=r, columns=list('xyz'), index=n*n*n*['Ar'])\n",
    "            df.to_csv(f_file, header=False, sep=' ', float_format='%.5f')\n",
    "\n",
    "        if i >= S_o:\n",
    "            mean_temperature += T\n",
    "            mean_pressure += pres\n",
    "            mean_Energy += E\n",
    "        \n",
    "    \n",
    "        \n",
    "mean_temperature/=S_d\n",
    "mean_pressure/=S_d\n",
    "mean_Energy/=S_d\n",
    "print(\"T_av\", mean_temperature, \"P_av\", mean_pressure, \"mean_Energy\", mean_Energy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Presenting Energy(time step) plot\n",
    "from pylab import *\n",
    "\n",
    "s = 0\n",
    "E_tab = []\n",
    "t_tab = []\n",
    "\n",
    "for i in range(int(S_d+S_o)):\n",
    "    s+=1\n",
    "    p = p+0.5*F*tau\n",
    "    r = r+p*tau/m\n",
    "    F, Epot, pres = potentials(r)\n",
    "    p = p+0.5*F*tau\n",
    "    Ekin = (p**2).sum(axis=1)/(2*m)\n",
    "    Ekin_total = Ekin.sum()\n",
    "    E = Epot + Ekin_total\n",
    "    E_tab.append(E)\n",
    "    T = 2/(3*N*k)*Ekin_total\n",
    "    t = i*tau\n",
    "    t_tab.append(t)\n",
    "    \n",
    "plot(t_tab, E_tab)\n",
    "plt.title(data[7])\n",
    "\n",
    "locs,labels = xticks()\n",
    "xticks(locs, map(lambda x: \"%g\" % x, locs))\n",
    "xlabel('time step')\n",
    "\n",
    "\n",
    "locs,labels = yticks()\n",
    "yticks(locs, map(lambda x: \"%.1f\" % x, locs))\n",
    "ylabel('energy')\n",
    "\n",
    "show()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}