Dinámica_Molecular/Triple_Pozo_2D.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Triple Pozo en 2D\n",
    "\n",
    "El potencial está todavía por reajustar. Esto se hizo con apoyo del notebook Gaussian_kernel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sympy as sy\n",
    "import matplotlib.pyplot as plt\n",
    "import simtk.openmm as mm\n",
    "import simtk.unit as unit\n",
    "import simtk.openmm.app as app\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('ggplot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Definición del sistema"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = mm.System()\n",
    "mass   = 39.948 * unit.amu # masa del átomo Ar\n",
    "system.addParticle(mass)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Triple Pozo en 2D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pozo harmónico en Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# '0.5*K*x**2'\n",
    "\n",
    "def harmonic_well(x,K):\n",
    "    return 0.5*K*x**2\n",
    "\n",
    "K=1.0\n",
    "\n",
    "x_serie = np.arange(-5., 5., 0.05)\n",
    "plt.rcParams['figure.figsize'] = 6,4\n",
    "plt.plot(x_serie, harmonic_well(x_serie,K), 'r-')\n",
    "plt.ylim(-1,7)\n",
    "plt.xlim(-10,10)\n",
    "plt.xlabel(\"Y or Z\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculamos el periodo de oscilación natural para tener referencia temporal para la elección del paso de integración. Igual que en los notebooks anteriores. Pero esta vez, como ya sabemos que en el caso del pozo harmónico k_effectiva=k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "effective_k = K\n",
    "T = 2*np.pi*np.sqrt(mass/(effective_k * unit.kilocalories_per_mole/unit.angstrom**2))\n",
    "print('Periodo T de oscilación',T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Triple Pozo en 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def elliptical_gaussian_kernel_2D(x,y,xo,yo,A=-1.0,cxx=0.0,cxy=0.0,cyy=0.0):\n",
    "    inc_x = x-xo\n",
    "    inc_y = y-yo\n",
    "    B = cxx*inc_x*inc_x + 2*cxy*inc_x*inc_y + cyy*inc_y*inc_y\n",
    "    return A*np.exp(-B)\n",
    "\n",
    "def check_positive_2d(embudo):\n",
    "    pot_i=0\n",
    "    for pot in embudo:\n",
    "        xo, yo, A, cxx, cxy, cyy = pot\n",
    "        M = np.array([[cxx, cxy],[cxy, cyy]])\n",
    "        is_positive = np.all(np.linalg.eigvals(M+M.transpose()) > 0)\n",
    "        if not is_positive:\n",
    "            print('La gaussiana',pot_i, 'no es definida positiva')\n",
    "        pot_i+=1\n",
    "\n",
    "def effective_ks_2D(embudo):\n",
    "    \n",
    "    effective_ks = []\n",
    "    pot_i = 0\n",
    "    for pot in embudo:\n",
    "        xo, yo, A, cxx, cxy, cyy = pot\n",
    "        x = sy.symbols('x')\n",
    "        y = sy.symbols('y')\n",
    "        f = A*sy.exp(-(cxx*(x-xo)**2 + 2*cxy*(x-xo)*(y-yo) + cyy*(y-yo)**2))\n",
    "        gradient = sy.derive_by_array(f,(x,y))\n",
    "        hessian = sy.derive_by_array(gradient,(x,y))\n",
    "        M = hessian.subs(x,xo).subs(y,yo)\n",
    "        M = np.array(M.tolist(),dtype=float)\n",
    "        eigenvals, eigenvecs = np.linalg.eig(np.array(M))\n",
    "        #print('Gaussiana',pot_i)\n",
    "        for ii in range(eigenvals.shape[0]):\n",
    "            eff_k = eigenvals[ii]\n",
    "            if eff_k > 0.0:\n",
    "                effective_ks.append(eff_k)\n",
    "            #print('\\t effective k',eff_k,'en modo',eigenvecs[:,ii])\n",
    "        pot_i+=1\n",
    "        #print('')\n",
    "    return np.max(effective_ks)\n",
    "\n",
    "def plot_embudo_2D(embudo):\n",
    "    \n",
    "    xlim = [-2,1]\n",
    "    ylim = [-1,2]\n",
    "    xbins = 100\n",
    "    ybins = 100\n",
    "\n",
    "    x = np.linspace(xlim[0], xlim[1], xbins)\n",
    "    y = np.linspace(ylim[0], ylim[1], ybins)\n",
    "\n",
    "    X,Y = np.meshgrid(x, y)\n",
    "    Z = 0.0\n",
    "    for pot in embudo:\n",
    "        xo, yo, A, cxx, cxy, cyy = pot\n",
    "        Z += elliptical_gaussian_kernel_2D(X,Y,xo,yo,A,cxx,cxy,cyy)\n",
    "\n",
    "    plt.rcParams['figure.figsize'] = 12, 12\n",
    "    extent = [xlim[0], xlim[1], ylim[0], ylim[1]]\n",
    "    #im = plt.imshow(Z, extent=extent, origin='lower',\n",
    "    #                cmap=plt.cm.RdBu, interpolation='gaussian')\n",
    "\n",
    "    cset = plt.contour(Z, np.arange(-4,0,0.1), extent=extent, origin='lower',\n",
    "                       linewidths=1, cmap=plt.cm.Set1)\n",
    "\n",
    "    plt.clabel(cset,inline=True,fmt='%1.1f',fontsize=10)\n",
    "    #plt.colorbar(im, alpha=0.5)\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pot_i = [xo, yo, A, cxx, cxy, cyy]\n",
    "pot_0 =  [ 0.0,  0.0, -4, 6.0, 0.0, 6.0]\n",
    "pot_1 =  [ -1.0, 0.0, -2, 6.0, 0.0, 6.0]\n",
    "pot_2 =  [ 0.0,  0.8, -3, 6.0, 0.0, 6.0]\n",
    "\n",
    "embudo = [pot_0, pot_1, pot_2]\n",
    "\n",
    "check_positive_2d(embudo)\n",
    "plot_embudo_2D(embudo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculamos el periodo de oscilación natural en lo mínimos de los kernel gaussianos usados para construir el potencial embudo. Tomamos la máxima k efectiva para tener referencia temporal para la elección del paso de integración."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "effective_k = effective_ks_2D(embudo)\n",
    "T = 2*np.pi*np.sqrt(mass/(effective_k * unit.kilocalories_per_mole/unit.angstrom**2))\n",
    "print('Periodo T de oscilación',T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Añadiendo el potencial al sistema"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Harmónico en eje Z\n",
    "\n",
    "force = mm.CustomExternalForce('(K/2.0)*(z^2)')\n",
    "force.addGlobalParameter('K', K * unit.kilocalories_per_mole/unit.angstrom**2)\n",
    "force.addParticle(0, [])\n",
    "system.addForce(force)\n",
    "\n",
    "### Potencial tipo embudo en X e Y\n",
    "\n",
    "gaussian_kernel = 'A*exp(-(cxx*(x-xo)^2+2.0*cxy*(x-xo)*(y-yo)+cyy*(y-yo)^2))'\n",
    "\n",
    "for pot in embudo:\n",
    "    xo, yo, A, cxx, cxy, cyy = pot\n",
    "    force = mm.CustomExternalForce(gaussian_kernel)\n",
    "    force.addGlobalParameter('A', A * unit.kilocalories_per_mole)\n",
    "    force.addGlobalParameter('cxx', cxx * 1.0/unit.angstrom**2)\n",
    "    force.addGlobalParameter('cxy', cxy * 1.0/unit.angstrom**2)\n",
    "    force.addGlobalParameter('cyy', cyy * 1.0/unit.angstrom**2)\n",
    "    force.addGlobalParameter('xo', xo * unit.angstrom)\n",
    "    force.addGlobalParameter('yo', yo * unit.angstrom)\n",
    "    force.addParticle(0, [])\n",
    "    system.addForce(force)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estado termodinámico"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Formalismo NVT\n",
    "kB = unit.BOLTZMANN_CONSTANT_kB * unit.AVOGADRO_CONSTANT_NA\n",
    "temperature = 300*unit.kelvin\n",
    "pressure    = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parámetros de la simulación"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "step_size       = 0.020*unit.picoseconds\n",
    "num_steps       = 40000\n",
    "saving_period   = 100\n",
    "num_steps_saved = int(num_steps/saving_period)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Integrador"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "friction   = 10.0/unit.picosecond\n",
    "integrator = mm.LangevinIntegrator(temperature, friction, step_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plataforma de simulación"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "platform_name = 'CPU'  #platform:     A platform (CPU, OpenCL, CUDA, or reference); default is platform=OpenCL\"\n",
    "# for ii in range(mm.Platform.getNumPlatforms()):\n",
    "#     print(mm.Platform.getPlatform(ii).getName())\n",
    "platform = mm.Platform.getPlatformByName(platform_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reporteros y arrays de salida"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "times      = unit.Quantity(np.zeros([num_steps_saved], np.float32), unit.picoseconds)\n",
    "positions  = unit.Quantity(np.zeros([num_steps_saved,3], np.float32), unit.angstroms)\n",
    "velocities = unit.Quantity(np.zeros([num_steps_saved,3], np.float32), unit.angstroms/unit.picosecond)\n",
    "potential_energies   = unit.Quantity(np.zeros([num_steps_saved], np.float32), unit.kilocalories_per_mole)\n",
    "kinetic_energies     = unit.Quantity(np.zeros([num_steps_saved], np.float32), unit.kilocalories_per_mole)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Condiciones iniciales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_positions  = [[0.0, 0.0, 0.0]] * unit.angstroms\n",
    "#initial_velocities = None # Las velocidades serán adjudicadas aleatoriamente según la distribución Maxwell-Boltzmann del estado termodinámico\n",
    "\n",
    "context = mm.Context(system, integrator, platform)\n",
    "context.setPositions(initial_positions)\n",
    "context.setVelocitiesToTemperature(temperature)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "state = context.getState(getEnergy=True, getPositions=True, getVelocities=True)\n",
    "times[0] = state.getTime()\n",
    "positions[0] = state.getPositions()[0]\n",
    "velocities[0] = state.getVelocities()[0]\n",
    "kinetic_energies[0]=state.getKineticEnergy()\n",
    "potential_energies[0]=state.getPotentialEnergy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Corriendo la simulación"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in tqdm(range(num_steps_saved)):\n",
    "    context.getIntegrator().step(saving_period)\n",
    "    state = context.getState(getEnergy=True, getPositions=True, getVelocities=True)\n",
    "    times[ii] = state.getTime()\n",
    "    positions[ii] = state.getPositions()[0]\n",
    "    velocities[ii] = state.getVelocities()[0]\n",
    "    kinetic_energies[ii]=state.getKineticEnergy()\n",
    "    potential_energies[ii]=state.getPotentialEnergy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Análisis de resultados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = 18, 4\n",
    "for ii, ylabel in zip(range(3),['X','Y','Z']):\n",
    "    plt.plot(times,positions[:,ii])\n",
    "    plt.ylabel(ylabel)\n",
    "    plt.xlabel('time')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist, bin_edges = np.histogram(positions[:,0],bins=100, density=True)\n",
    "plt.plot((bin_edges[1:]+bin_edges[:-1])/2.0,hist)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot((bin_edges[1:]+bin_edges[:-1])/2.0,-np.log(hist))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_frames = positions.shape[0]\n",
    "averages   = unit.Quantity(np.zeros(num_frames, np.float32), unit.angstroms)\n",
    "std_devs   = unit.Quantity(np.zeros(num_frames, np.float32), unit.angstroms)\n",
    "    \n",
    "for ii in np.arange(num_frames):\n",
    "    averages[ii] = np.mean(positions[:(ii+1),0])\n",
    "    std_devs[ii]  = np.std(positions[:(ii+1),0])\n",
    "\n",
    "plt.plot(times,averages)\n",
    "plt.show()\n",
    "\n",
    "plt.plot(times,std_devs)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "writing = False\n",
    "\n",
    "if writing: \n",
    "    h5f = h5py.File('data_double_well_md.h5', 'w')\n",
    "    h5f.create_dataset('times', data=times)\n",
    "    h5f.create_dataset('positions', data=positions)\n",
    "    h5f.create_dataset('velocities', data=velocities)\n",
    "    h5f.create_dataset('kinetic_energies', data=kinetic_energies)\n",
    "    h5f.create_dataset('potential_energies', data=potential_energies)\n",
    "    h5f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reading =False\n",
    "\n",
    "if reading:\n",
    "    h5f = h5py.File('data_double_well_md.h5', 'r')\n",
    "    #for name in h5f:\n",
    "    #    print(name)\n",
    "    times              = h5f['times'][:]\n",
    "    positions          = h5f['positions'][:]\n",
    "    velocities         = h5f['velocities'][:]\n",
    "    kinetic_energies   = h5f['kinetic_energies'][:]\n",
    "    potential_energies = h5f['potential_energies'][:]\n",
    "    h5f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "T = sqrt(m/K)\n",
    "timestep smaller than ~ T/10\n",
    "standard deviation in each dimension sigma= (kT / K)^(1/2)\n",
    "Expectation and standard deviation of the potential energy of a 3D harmonic oscillator is (3/2)kT\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Todo:\n",
    "\n",
    "- Los gaussian kernel son demasiado abiertos, la k efectiva es muy baja en comparación con el potencial del eje z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.7.1"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}