Check_tests2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e056a445-3aee-47f7-9f8c-fad94aca97bc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f8eaa4a2aaf1449f9e08d8e16291b5f9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pickle\n",
    "import gzip\n",
    "from matplotlib import pyplot as plt\n",
    "import networkx as nx\n",
    "from tqdm import tqdm\n",
    "import numpy as np\n",
    "import molsysmt as msm\n",
    "from openmm.app import *\n",
    "from openmm import *\n",
    "from openmm.unit import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "514b8168-3184-453d-8c23-d064f3e1a1e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "with gzip.open(f'tests/energies_0.pkl.gz', 'r') as fff:\n",
    "    energy_dbs=pickle.load(fff)\n",
    "        \n",
    "with gzip.open(f'tests/traj_inh_0.pkl.gz', 'r') as fff:\n",
    "    traj_inh=pickle.load(fff)\n",
    "\n",
    "with gzip.open(f'tests/coors_0.pkl.gz', 'r') as fff:\n",
    "    coors_dbs=pickle.load(fff)\n",
    "\n",
    "with gzip.open(f'tests/dihed_0.pkl.gz', 'r') as fff:\n",
    "    dihed_dbs=pickle.load(fff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9668adba-4359-4bd2-9dbc-def3180ad0e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_1 = 4486\n",
    "min_2 = 575"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4caf2eed-3fc1-4c37-b822-92c29bef7c27",
   "metadata": {},
   "outputs": [],
   "source": [
    "molsys1 = msm.convert(\"pdbs/V5.pdb\")\n",
    "new_coordinates = msm.pyunitwizard.quantity(value=coors_dbs[min_1], unit='nm')\n",
    "msm.set(molsys1, element='atom', coordinates=new_coordinates)\n",
    "\n",
    "molsys2 = msm.convert(\"pdbs/V5.pdb\")\n",
    "new_coordinates = msm.pyunitwizard.quantity(value=coors_dbs[min_2], unit='nm')\n",
    "msm.set(molsys2, element='atom', coordinates=new_coordinates)\n",
    "molsys2 = msm.structure.least_rmsd_fit(molecular_system=molsys2, reference_molecular_system=molsys1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c854a73b-27bd-4421-b3b4-9cf27bbaf4f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def interpolate_images(pos_min1, pos_min2, num_images):\n",
    "    images = [pos_min1 + i/(num_images-1)*(pos_min2 - pos_min1) for i in range(num_images)]\n",
    "    return images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64296c04-8998-4bca-94f6-240c8b5b9841",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_spring_forces(images, k_spring):\n",
    "    spring_forces = []\n",
    "    for i in range(1, len(images)-1):  # No modificar extremos\n",
    "        prev_image = images[i-1]\n",
    "        next_image = images[i+1]\n",
    "        current_image = images[i]\n",
    "        \n",
    "        tangent = next_image - prev_image\n",
    "        tangent /= np.linalg.norm(tangent)  # Normaliza el vector tangente\n",
    "\n",
    "        # Proyecta las fuerzas en componentes paralelas y perpendiculares\n",
    "        spring_force = k_spring * (np.linalg.norm(next_image - current_image) -\n",
    "                                   np.linalg.norm(current_image - prev_image)) * tangent\n",
    "        spring_forces.append(spring_force)\n",
    "    return spring_forces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0b78f28-515f-47bc-9615-d68da47aba73",
   "metadata": {},
   "outputs": [],
   "source": [
    "pdb_filename = \"pdbs/V5.pdb\"\n",
    "pdb = PDBFile(pdb_filename)\n",
    "forcefield = ForceField(\"amber14-all.xml\", \"amber14/tip3p.xml\")\n",
    "system = forcefield.createSystem(\n",
    "    pdb.topology,\n",
    "    nonbondedMethod=NoCutoff,  # Sin tratamiento de fronteras periódicas\n",
    "    constraints=HBonds         # Restringe enlaces de hidrógeno para un paso de integración más largo\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02985a8f-6689-418b-9862-d4ac35322578",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Crear simulaciones para cada imagen\n",
    "num_images = 10\n",
    "positions_min1 = msm.get(molsys1, coordinates=True)  # Conformación mínima 1\n",
    "positions_min2 = msm.get(molsys2, coordinates=True)  # Conformación mínima 2\n",
    "images = interpolate_images(positions_min1, positions_min2, num_images)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab417c91-0925-4aa3-87fd-29a559d4f035",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulations = []\n",
    "for i in range(num_images):\n",
    "    integrator = VerletIntegrator(1.0*femtoseconds)\n",
    "    simulation = Simulation(pdb.topology, system, integrator)\n",
    "    simulation.context.setPositions(msm.pyunitwizard.convert(images[i][0],to_form='openmm.unit'))\n",
    "    simulations.append(simulation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "532deb74-a7e2-49d4-836e-6ef2dc0977ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_positions_with_forces(positions, forces, step_size=0.001 * nanometer):\n",
    "    \"\"\"\n",
    "    Actualiza las posiciones en base a las fuerzas.\n",
    "    \n",
    "    Parámetros:\n",
    "        positions: numpy.array con las posiciones actuales (en nm).\n",
    "        forces: numpy.array con las fuerzas actuales (en kJ/mol/nm).\n",
    "        step_size: Escalamiento del desplazamiento (en nm).\n",
    "    \"\"\"\n",
    "    displacement = step_size * (msm.pyunitwizard.get_value(forces / np.linalg.norm(forces, axis=1)[:, None]))\n",
    "    return positions + displacement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "660b05b1-9813-4236-9ed8-4a4e523c4b3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_neb(simulations, k_spring=100*kilojoule_per_mole/nanometer**2, max_iterations=100):\n",
    "    for iteration in range(max_iterations):\n",
    "        # Obtener posiciones actuales de todas las imágenes\n",
    "        positions = [sim.context.getState(getPositions=True).getPositions(asNumpy=True) for sim in simulations]\n",
    "        \n",
    "        # Calcular fuerzas físicas y de resorte\n",
    "        spring_forces = calculate_spring_forces(positions, k_spring)\n",
    "        \n",
    "        for i, sim in enumerate(simulations[1:-1], start=1):  # Excluye los extremos\n",
    "            state = sim.context.getState(getForces=True)\n",
    "            physical_forces = state.getForces(asNumpy=True)\n",
    "            \n",
    "            # Combinar fuerzas físicas y de resorte\n",
    "            total_forces = physical_forces + spring_forces[i-1]\n",
    "            sim.context.setVelocitiesToTemperature(0.0)  # Congela el sistema\n",
    "            positions[i] = update_positions_with_forces(positions[i], total_forces)\n",
    "            sim.context.setPositions(positions[i])\n",
    "\n",
    "        # Criterio de convergencia (opcional)\n",
    "        max_force = max(np.linalg.norm(f) for f in spring_forces)\n",
    "        print(max_force)\n",
    "        if max_force < 1e-3:\n",
    "            break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2b64aeb-3c84-41ed-bdb9-8923c11d9858",
   "metadata": {},
   "outputs": [],
   "source": [
    "run_neb(simulations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8e0c867-c807-468b-a0dc-0da40d93c144",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63154c40-79cf-4897-9ad7-31c8bd03f7e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_dbs[min_1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c1404af-552c-4a4e-970b-7340fb48a471",
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_dbs[min_2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3afc6555-f062-486f-8710-60979e10ea1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "aux = msm.merge([molsys1,molsys2])\n",
    "msm.view(aux, standard=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1995c7b5-7c6e-45ce-8c70-629d2b9965d3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1b5a4cb-5d11-4190-b580-2178024b52cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "path = interpolate_images(molsys1, molsys2, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "863e1e0a-6af3-4f4c-8d9d-d5735b37b013",
   "metadata": {},
   "outputs": [],
   "source": [
    "msm.view(msm.merge(path), standard=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91b1198c-ac76-4fdc-9d17-8ae2b974401e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_spring_forces(path, k_spring):\n",
    "\n",
    "    images=[]\n",
    "    for sys_aux in path:\n",
    "        images.append(msm.get(sys_aux, coordinates=True))\n",
    "        \n",
    "    spring_forces = []\n",
    "    for i in range(1, len(images)-1):  # No modificar extremos\n",
    "        prev_image = images[i-1]\n",
    "        next_image = images[i+1]\n",
    "        current_image = images[i]\n",
    "        \n",
    "        tangent = next_image - prev_image\n",
    "        tangent /= np.linalg.norm(tangent)  # Normaliza el vector tangente\n",
    "\n",
    "        # Proyecta las fuerzas en componentes paralelas y perpendiculares\n",
    "        spring_force = k_spring * (np.linalg.norm(next_image - current_image) -\n",
    "                                   np.linalg.norm(current_image - prev_image)) * tangent\n",
    "        spring_forces.append(spring_force)\n",
    "        \n",
    "    return spring_forces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a16583f9-d562-48bf-9c0a-42f0169f1b52",
   "metadata": {},
   "outputs": [],
   "source": [
    "path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c403c8a-0f6b-42ca-8dab-40800e3fe786",
   "metadata": {},
   "outputs": [],
   "source": [
    "pesos_transiciones={}\n",
    "pesos_nodos={}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d217b75d-3422-43a9-9a45-32b964c132bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(len(traj_inh)-1):\n",
    "    nodo_1 = traj_inh[ii]\n",
    "    nodo_2 = traj_inh[ii+1]\n",
    "    if (nodo_1, nodo_2) in pesos_transiciones:\n",
    "        pesos_transiciones[nodo_1, nodo_2]+=1\n",
    "    else:\n",
    "        pesos_transiciones[nodo_1, nodo_2]=1\n",
    "    if nodo_1 in pesos_nodos:\n",
    "        pesos_nodos[nodo_1]+=1\n",
    "    else:\n",
    "        pesos_nodos[nodo_1]=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "296829d2-d836-48e1-a812-263b0bce6218",
   "metadata": {},
   "outputs": [],
   "source": [
    "G = nx.DiGraph()  # Para un grafo dirigido"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25cb0b85-56c0-455e-a89b-8f735b57b0f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Agregar las aristas desde el diccionario\n",
    "for nodes, peso in pesos_transiciones.items():\n",
    "    G.add_edge(nodes[0], nodes[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d54a19b1-30ae-4815-8087-15a2ec952fbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "G[4486]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96fd5d81-7266-4fd5-851d-01736d4c25b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.argmax(pesos_nodos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e7e7c7f-9df4-451c-8b20-139b5a4d90d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "peso_max=0\n",
    "nodo_max=0\n",
    "for ii,jj in pesos_nodos.items():\n",
    "    if jj>peso_max:\n",
    "        nodo_max=ii\n",
    "        peso_max=jj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28f457a2-f9ba-44b4-b675-c111e961c916",
   "metadata": {},
   "outputs": [],
   "source": [
    "nodo_max"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61c83f53-f406-4238-b031-17da68b0cadf",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_1 = np.argmin(energy_dbs[1])\n",
    "energy_dbs[1][min_1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d48932f-0052-4b65-8cd3-730a9f7cbcf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_2 = np.argmin(energy_dbs[2])\n",
    "energy_dbs[2][min_2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0f30f95-6b50-4790-a214-587467563228",
   "metadata": {},
   "outputs": [],
   "source": [
    "dihed_dbs[0][min_0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "990bd436-bb67-480a-8b2b-d5c437a0dfaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "dihed_dbs[1][min_1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae865a20-c577-4a62-8597-981ab08c55d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "dihed_dbs[2][min_2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81b6fbf7-eed7-49ba-a1d2-2bc1ecc859ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Nombre del archivo PDB del péptido\n",
    "molsys0 = msm.convert(\"pdbs/V5.pdb\")\n",
    "molsys1 = msm.convert(\"pdbs/V5.pdb\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73cbcab3-0dc5-4197-9151-5432a251c5bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "new_coordinates = msm.pyunitwizard.quantity(value=coors_dbs[min_0], unit='nm')\n",
    "msm.set(molsys0, element='atom', coordinates=new_coordinates)\n",
    "\n",
    "new_coordinates = msm.pyunitwizard.quantity(value=coors_dbs[575], unit='nm')\n",
    "msm.set(molsys1, element='atom', coordinates=new_coordinates)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aedb49be-8799-4151-a755-cc7fa5fb39f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "molsys01 = msm.structure.least_rmsd_fit(molecular_system=molsys0, reference_molecular_system=molsys1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c558ed4-f4aa-49d7-98fd-42270930dba1",
   "metadata": {},
   "outputs": [],
   "source": [
    "aux = msm.merge([molsys01,molsys1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb335392-5fb8-4efc-a136-1c079ef82d7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "msm.view(aux, standard=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b0cd6a9-d1d9-4810-957c-e3bcc2d74fac",
   "metadata": {},
   "outputs": [],
   "source": [
    "G = nx.DiGraph()  # Para un grafo dirigido"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b96b4044-b1cd-4c89-9f77-301d30e428a5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bd0f32e-9f06-4867-8b88-d915ae2e87d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def es_el_mismo(angs1, angs2):\n",
    "    if np.max(np.abs(angs1-angs2))<0.1:\n",
    "        return True\n",
    "    else:\n",
    "        return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccceb26b-e2c9-4228-9ad5-f53422b25668",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee9a995b-dbed-42a6-aeb5-f0e04caf2ee2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6748d8a1-0503-4b59-9780-445453c7b9af",
   "metadata": {},
   "outputs": [],
   "source": [
    "pesos_transiciones={}\n",
    "pesos_nodos={}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c774f5ad-ed4f-4263-a69f-e7bb1f7b99ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(len(traj_inh)-1):\n",
    "    nodo_1 = traj_inh[ii]\n",
    "    nodo_2 = traj_inh[ii+1]\n",
    "    if (nodo_1, nodo_2) in pesos_transiciones:\n",
    "        pesos_transiciones[nodo_1, nodo_2]+=1\n",
    "    else:\n",
    "        pesos_transiciones[nodo_1, nodo_2]=1\n",
    "    if nodo_1 in pesos_nodos:\n",
    "        pesos_nodos[nodo_1]+=1\n",
    "    else:\n",
    "        pesos_nodos[nodo_1]=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad61f60f-1f7b-41e4-b6c8-f790fbcbd06d",
   "metadata": {},
   "outputs": [],
   "source": [
    "G = nx.DiGraph()  # Para un grafo dirigido"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ee9ca50-98bf-4fd9-991b-c982c0474bfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Agregar las aristas desde el diccionario\n",
    "for nodes, peso in pesos_transiciones.items():\n",
    "    G.add_edge(nodes[0], nodes[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "125ea653-a074-415a-be40-5abf04988735",
   "metadata": {},
   "outputs": [],
   "source": [
    "nx.draw(G, with_labels=False, node_color=\"blue\", node_size=0.1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca95ce3e-f278-4214-a7bb-364aac352104",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}