Two new openmm forces: allow and forbidden z region

docs/contents/user/tools/thirds/openmm/forces/forbidden_z_region.ipynb

@@ -0,0 +1,284 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "653bc9c8-3e1b-45ee-944d-85874a8b33b3",
+   "metadata": {},
+   "source": [
+    "# Forbidden z region"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "64de0002-d07f-4bb1-8c28-b7c0f3e8762f",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "3b0bd911055c473a91e784a05fc00011",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": []
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "import molsysmt as msm"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "6e09021e-88d6-4ade-9ae5-dfa9ec540214",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import openmm as mm\n",
+    "from openmm import unit\n",
+    "from openmm import app\n",
+    "from tqdm import tqdm\n",
+    "import numpy as np\n",
+    "from matplotlib import pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "96637603-6cd0-49f4-8d8c-bd692f1d51e5",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "system = mm.System()\n",
+    "system.addParticle(39.948 * unit.amu) # masa del átomo de argón"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "991d580d-2f65-444e-a907-c0b94fb9de9b",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "msm.thirds.openmm.forces.forbidden_z_region(system, force_constant = '500 kilojoules/(mol*angstroms**2)',\n",
+    "                                          pbc=True, adding_force=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "2dd1d2d8-d85e-4396-8fa6-d3920bee98d1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Formalismo NVT\n",
+    "temperature = 300*unit.kelvin\n",
+    "pressure    = None"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ae200579-9af9-40de-aba3-7dce8fa33e54",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "integration_timestep = 2.0*unit.femtoseconds\n",
+    "saving_timestep = 1.00*unit.picoseconds\n",
+    "simulation_time = 2500.*unit.picoseconds\n",
+    "\n",
+    "saving_steps = int(saving_timestep/integration_timestep)\n",
+    "num_saving_steps = int(simulation_time/saving_timestep)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c3ca01e5-f1d4-4881-8b41-7f140dd441b7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "friction   = 5.0/unit.picoseconds\n",
+    "integrator = mm.LangevinIntegrator(temperature, friction, integration_timestep)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "43e38fd7-6f23-44e4-951d-9d9f8149cc0e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "platform = mm.Platform.getPlatformByName('CPU')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a3287fd3-cb9a-434e-b4e5-2b98a1b4478d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "times = np.zeros(num_saving_steps, np.float32) * unit.picoseconds\n",
+    "positions  = np.zeros([num_saving_steps,3], np.float32) * unit.angstroms\n",
+    "velocities = np.zeros([num_saving_steps,3], np.float32) * unit.angstroms/unit.picosecond\n",
+    "potential_energies   = np.zeros([num_saving_steps], np.float32) * unit.kilocalories_per_mole\n",
+    "kinetic_energies     = np.zeros([num_saving_steps], np.float32) * unit.kilocalories_per_mole"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5c69b0b6-93d0-4e93-901d-df21b34cf525",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "initial_positions  = [[0.0, 0.0, 2.0]] * unit.nanometers"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a10b47a1-fa8f-480c-9261-8ecb259a2352",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "context = mm.Context(system, integrator, platform)\n",
+    "context.setPositions(initial_positions)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8f358f9e-faab-4b93-9dc6-16f58dad933a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "L = 4.0\n",
+    "v1 = [L,0,0] * unit.nanometers\n",
+    "v2 = [0,L,0] * unit.nanometers\n",
+    "v3 = [0,0,L] * unit.nanometers\n",
+    "L = L * unit.nanometers\n",
+    "context.setPeriodicBoxVectors(v1, v2, v3)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b4f8fa76-c874-4fdc-9f75-ec40f4579943",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "state = context.getState(getEnergy=True, getPositions=True, getVelocities=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "97849832-985e-41ec-bf4b-bb18ebf817eb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "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": "code",
+   "execution_count": null,
+   "id": "4559e045-b60d-46d0-a0dc-2f988edfcb6d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for ii in tqdm(range(1,num_saving_steps)):\n",
+    "    context.getIntegrator().step(saving_steps)\n",
+    "    state_xx = context.getState(getEnergy=True, getPositions=True, getVelocities=True)\n",
+    "    times[ii] = state_xx.getTime()\n",
+    "    positions[ii] = state_xx.getPositions()[0]\n",
+    "    velocities[ii] = state_xx.getVelocities()[0]\n",
+    "    kinetic_energies[ii]=state_xx.getKineticEnergy()\n",
+    "    potential_energies[ii]=state_xx.getPotentialEnergy()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6f0c68d0-3191-4266-b338-789328039407",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#plt.plot(positions[:,0], label='X')\n",
+    "#plt.plot(positions[:,1], label='Y')\n",
+    "plt.plot(positions[:,2], label='Z')\n",
+    "plt.legend()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c5af5804-00df-44e1-8eef-d987bdf7d994",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2bc40b17-07c8-4835-b662-1ac15001fea9",
+   "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.10.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

docs/contents/user/tools/thirds/openmm/forces/index.md

@@ -0,0 +1,11 @@
+# Custom External Forces
+
+```{eval-rst}
+.. toctree::
+   :maxdepth: 2
+
+   allowed_z_region.ipynb
+   forbidden_z_region.ipynb
+   harmonic_potential_to_plane.ipynb
+   
+```

docs/contents/user/tools/thirds/openmm/index.md

@@ -4,4 +4,6 @@
 .. toctree::
    :maxdepth: 2
 
+   forces/index.md
+   
 ```

molsysmt/_private/digestion/argument/width.py

@@ -0,0 +1,12 @@
+import numpy as np
+from molsysmt import pyunitwizard as puw
+from ...exceptions import ArgumentError
+
+def digest_width(width, caller=None):
+
+    if puw.is_quantity(width):
+        if puw.check(width, dimensionality={'[L]':1}):
+            return puw.standardize(width)
+
+    raise ArgumentError('width', value=width, caller=caller, message=None)
+

molsysmt/_private/digestion/argument/z0.py

@@ -0,0 +1,12 @@
+import numpy as np
+from molsysmt import pyunitwizard as puw
+from ...exceptions import ArgumentError
+
+def digest_z0(z0, caller=None):
+
+    if puw.is_quantity(z0):
+        if puw.check(z0, dimensionality={'[L]':1}):
+            return puw.standardize(z0)
+
+    raise ArgumentError('z0', value=z0, caller=caller, message=None)
+

molsysmt/_version.py

@@ -1 +1 @@
-__version__ = "0.8.3+15.g410eec0f.dirty"
+__version__ = "0.8.3+17.gcda3efaf.dirty"

molsysmt/form/molsysmt_MSMH5FileHandler/get.py

@@ -14,11 +14,7 @@
 @digest(form=form)
 def get_atom_id_from_atom(item, indices='all'):
 
-    output = None
-
     output = item.file['topology']['atoms']['id'][:].astype('int64')
-    if is_all(indices):
-
 
     raise NotImplementedError
 

molsysmt/form/openmm_System/attributes.py

@@ -2,6 +2,7 @@
 
 attributes = {ii:False for ii in _all_attributes}
 
+attributes['n_atoms'] = True
 attributes['atom_index'] = True
 attributes['atom_id'] = True
 attributes['atom_name'] = True

molsysmt/form/openmm_System/get.py

@@ -305,13 +305,7 @@ def get_entity_type_from_entity(item, indices='all'):
 @digest(form=form)
 def get_n_atoms_from_system(item):
 
-    from . import to_openmm_Topology
-    from ..openmm_Topology import get_n_atoms_from_system as aux_get
-
-    tmp_item = to_openmm_Topology(item)
-    output = aux_get(tmp_item, indices=indices)
-
-    return output
+    return item.getNumParticles()
 
 @digest(form=form)
 def get_n_groups_from_system(item):

molsysmt/thirds/openmm/forces/__init__.py

@@ -3,3 +3,5 @@
 from .flat_harmonic_potential_to_plane import flat_harmonic_potential_to_plane
 from .harmonic_bonds_potential import harmonic_bonds_potential
 from .constant_pulling_force import constant_pulling_force
+from .allowed_z_region import allowed_z_region 
+from .forbidden_z_region import forbidden_z_region