diff --git a/.ipynb_checkpoints/1-basics-checkpoint.ipynb b/.ipynb_checkpoints/1-basics-checkpoint.ipynb
deleted file mode 100644
index c95e9e64a8911ae6af32756cacbdc4f2fbb3d2f6..0000000000000000000000000000000000000000
--- a/.ipynb_checkpoints/1-basics-checkpoint.ipynb
+++ /dev/null
@@ -1,311 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Basic usage of the `arbdmodel` package\n",
- "\n",
- "We'll learn how to build simple models using the arbdmodel Python package that will run on the Atomic Resolution Brownian Dynamics (ARBD/arbd) simulation engine. ARBD requires a GPU and linux operating system.\n",
- "\n",
- "This tutorial assumes you have a working knowledge of Python, including its object-oriented features, and the numpy package.\n",
- "\n",
- "Also, a word of caution: at the time of writing ARBD will not prevent you from simulating with poorly selected runtime parameters (e.g. too large a timestep; or too large a decomposition period with too small a pairlist distance; or simply potentials that prescribe unreasonably large forces). Hence it is crucial to have a good understanding of the MD/BD algorithm if you are using it in production settings."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Step 1: Create particle types "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "import ipdb # in case we need to debug, after a crash, run `ipdb.pm()`\n",
- "import numpy as np\n",
- "from arbdmodel import ParticleType, PointParticle, ArbdModel\n",
- "\n",
- "system_size = 620\n",
- "num_particles = 6400\n",
- "\n",
- "## diffusion of liquid argon: https://www.sciencedirect.com/science/article/pii/0375960171901290 \n",
- "## units \"60 * 3.41 AA * sqrt(119.5 k K / 40 dalton)\" \"AA**2/ns\" \n",
- "argon = ParticleType(name=\"Ar\",\n",
- " mass=39.95, # dalton\n",
- " damping_coefficient=5000, # per ns \n",
- " custom_property=\"my_value\",\n",
- " epsilon=0.177, # kcal_mol \n",
- " radius=3.345/2) # Angstroms\n",
- "\n",
- "print(argon)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Step 2: Build a system"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "## Generate Nx3 array of random cartesian coordinates\n",
- "positions = system_size*(np.random.random([num_particles,3])-0.5)\n",
- "\n",
- "## Create a list of point particles located at those positions\n",
- "points = [PointParticle(type_=argon, name=\"Ar\", position=positions[i])\n",
- " for i in range(num_particles)]\n",
- "\n",
- "model = ArbdModel(points, \n",
- " dimensions=[system_size for i in range(3)], # Ã…ngstroms\n",
- " timestep = 20e-6, # ns; can be specified below with engine instead\n",
- " )\n",
- "\n",
- "print(model)\n",
- "print(model.children[:2])"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Step 3: Describe the interactions between the particles \n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "from arbdmodel.interactions import LennardJones\n",
- "\n",
- "lj = LennardJones()\n",
- "model.add_nonbonded_interaction(lj)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Step 4: Run the simulation"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "from arbdmodel import ArbdEngine\n",
- "\n",
- "engine = ArbdEngine(output_period = 1e2, num_steps = 1e5)\n",
- "engine.simulate(model, output_name=\"1-lj\", directory='sims/1-argon',\n",
- " num_steps=1e3, gpu=1\n",
- ") # simulation parameters here override those in constructor \n"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Step 5: Visualize the results\n",
- "\n",
- "We recommend using [VMD](https://www.ks.uiuc.edu/Development/Download/download.cgi?PackageName=VMD) for visualizing the simulations trajectories with:\n",
- "```bash\n",
- "vmd sims/1-argon/1-lj.psf sims/1-argon/output/1-lj.dcd\n",
- "```\n",
- "However, for convenience we have provided commands to use the nglview module, which provides browser-based visualization.\n",
- "\n",
- "note: you must enable certain jupyter extensions to use nglview as follows before running jupyter-notebook\n",
- "```bash\n",
- "jupyter-nbextension enable --py --sys-prefix widgetsnbextension\n",
- "jupyter-nbextension enable nglview --py --sys-prefix\n",
- "```"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "import MDAnalysis as mda\n",
- "import nglview as nv\n",
- "\n",
- "u = mda.Universe(\"sims/1-argon/1-lj.pdb\", \"sims/1-argon/output/1-lj.dcd\")\n",
- "w = nv.show_mdanalysis(u)\n",
- "w.clear_representations()\n",
- "w.add_ball_and_stick(\"all\",radius=10)\n",
- "w"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Step 6: Customize the interactions"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "from arbdmodel.interactions import AbstractPotential\n",
- "\n",
- "## Clear nonbonded interactions \n",
- "model.nonbonded_interactions = []\n",
- "\n",
- "class BuckinghamPotential(AbstractPotential):\n",
- " def potential(self, r, types):\n",
- " ## https://royalsocietypublishing.org/doi/10.1098/rspa.1938.0173 \n",
- " ## optionally could pull information from typeA, typeB \n",
- " typeA, typeB = types\n",
- " return 143932.65 * ( 1.69 * np.exp( -r/0.273 ) - 102e-4 / r**6 )\n",
- "\n",
- "pot = BuckinghamPotential()\n",
- "model.add_nonbonded_interaction( pot )\n",
- "\n",
- "engine.simulate(model, output_name=\"2-buck\", directory='sims/1-argon',\n",
- " num_steps = 1e3, gpu=1\n",
- ")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "import MDAnalysis as mda\n",
- "import nglview as nv\n",
- "\n",
- "u = mda.Universe(\"sims/1-argon/2-buck.pdb\", \"sims/1-argon/output/2-buck.dcd\")\n",
- "w = nv.show_mdanalysis(u)\n",
- "w.clear_representations()\n",
- "w.add_ball_and_stick(\"all\",radius=10)\n",
- "w"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### Step 7: Add bonds\n",
- "\n",
- "Argon is, of course, a noble gas and will not form bonds; this is a demonstration only, and not representative of a physical system"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "from arbdmodel.interactions import HarmonicBond\n",
- "bond = HarmonicBond(k=100, # kcal/mol AA**2 \n",
- " r0 = 3 # AA \n",
- ")\n",
- "\n",
- "## Go through every other index \n",
- "for i in range(0,len(model.children),2):\n",
- " j = i+1\n",
- " model.add_bond( model.children[i],\n",
- " model.children[j],\n",
- " bond )\n",
- "\n",
- "print(\"numbonds:\", len(model.bonds))\n",
- "\n",
- "engine.simulate(model, output_name=\"3-bonds\", directory='sims/1-argon',\n",
- " num_steps = 1e3, gpu=1\n",
- ")\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "import MDAnalysis as mda\n",
- "import nglview as nv\n",
- "\n",
- "u = mda.Universe(\"sims/1-argon/3-bonds.psf\",\"sims/1-argon/3-bonds.pdb\", \"sims/1-argon/output/3-bonds.dcd\")\n",
- "w = nv.show_mdanalysis(u)\n",
- "w.clear_representations()\n",
- "w.add_ball_and_stick(\"index < 10\",radius=10)\n",
- "w"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "## Why are there so many long bonds?\n",
- "print(bond.range_) ## Distances are in Angstroms\n",
- "r = np.linspace(*bond.range_, 20)\n",
- "print(r)\n",
- "print(bond.potential(r)) ## Energies are kcal/mol\n",
- "\n",
- "## Outside of this range, the bond applies zero force"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### ARBD is configured through plain text files\n",
- "Use your favorite editor to examine the files named `sims/1-argon/*.bd`.\n",
- "ARBD uses *key-value* pairs to configure most features. The configuration parser does not use an interpretter, so math and string operations may not be used inside .bd files. For a list of the available keywords, see the documentation that comes with ARBD. \n",
- "\n",
- "Also, please note that ARBD is under development; some keywords have been deprecated, and the syntax is subject to evolve. The user-oriented interface of the arbdmodel package is less likely to change, which is one of the reasons we strongly recommend using arbdmodel to set up your simulation systems. Nevertheless, if you want to use ARBD, it's a good idea to familiarize yourself with the configuration file syntax."
- ]
- }
- ],
- "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.8.18"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/.ipynb_checkpoints/2-polymer-objects-checkpoint.ipynb b/.ipynb_checkpoints/2-polymer-objects-checkpoint.ipynb
deleted file mode 100644
index 9ed99f505a0198c48b5da1be83b9d010adfe24a8..0000000000000000000000000000000000000000
--- a/.ipynb_checkpoints/2-polymer-objects-checkpoint.ipynb
+++ /dev/null
@@ -1,597 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "58505291",
- "metadata": {},
- "outputs": [],
- "source": [
- "## These commands are useful if you are modifying the arbdmodel package while running the notebook\n",
- "%load_ext autoreload\n",
- "%autoreload 2"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "9461eaba",
- "metadata": {},
- "outputs": [],
- "source": [
- "import ipdb # in case we want to debug something\n",
- "\"\"\"\n",
- "If you want to enter function `fn` with the debugger, run:\n",
- "`ipdb.runcall(fn,*args,**kwargs)` instead of the usual `fn(*args,**kwargs)`\n",
- "\n",
- "If a cell crashes, you can (immediately after) create a new notebook\n",
- "cell, and add to it: `ipdb.pm()` to start the post-mortem debugger. \n",
- " \"\"\"\n",
- "\n",
- "import numpy as np\n",
- "from arbdmodel.polymer import PolymerSection, PolymerGroup\n",
- "\n",
- "polymer = PolymerSection(\"PROT\", num_monomers=70, monomer_length=3.8,\n",
- " start_position = np.array((20,0,-3*70/2)) )\n",
- "print(polymer)\n",
- "\n",
- "## Let's inspect the polymer object from within Python by printing it's dictionary of attributes\n",
- "for k,v in polymer.__dict__.items():\n",
- " print(f'{k}: {v}')\n",
- " "
- ]
- },
- {
- "cell_type": "markdown",
- "id": "cbf962fe",
- "metadata": {},
- "source": [
- "As you may be able to see above, the PolymerSection class describes a 1D polymeric object abstractly, tracing a linear path through 3D space using splines. It can also use splines to trace the orientation of the polymer through the `quaternion_spline_params` attribute.\n",
- "\n",
- "The object is agnostic to the underlying details of a bead-based model, which is useful when switching model resolutions.\n",
- "\n",
- "In order to do something with the polymer, we'll need to import or create a model. The arbdmodel package alreadly includes a few simple polymer models that we can use. The simplest is a freely-jointed chain model that connects beads by harmonic bonds. These classes use the basic elements of ParticleType and AbstractPotential to construct a bead-based model.\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "91973184",
- "metadata": {},
- "outputs": [],
- "source": [
- "from arbdmodel import ArbdEngine\n",
- "from arbdmodel.fjc_polymer_model import FjcModel\n",
- "\n",
- "model = FjcModel([polymer], monomers_per_bead_group=10)\n",
- "\n",
- "engine = ArbdEngine(output_period = 1e2, num_steps = 1e5)\n",
- "engine.simulate(model, output_name=\"1-linear-fjc\", directory='sims/2-polymer',\n",
- " num_steps=1e4, gpu=1\n",
- ")\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "258ee2b1",
- "metadata": {},
- "outputs": [],
- "source": [
- "import MDAnalysis as mda\n",
- "import nglview as nv\n",
- "\n",
- "u = mda.Universe(\"sims/2-polymer/1-linear-fjc.psf\", \"sims/2-polymer/output/1-linear-fjc.dcd\")\n",
- "w = nv.show_mdanalysis(u)\n",
- "w.clear_representations()\n",
- "w.add_ball_and_stick(\"all\",radius=1)\n",
- "w"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "dc99b017",
- "metadata": {},
- "source": [
- "Let's change the resolution of our model, being aware that the FJC properties are resolution dependent."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "045d0422",
- "metadata": {},
- "outputs": [],
- "source": [
- "from arbdmodel.coords import read_arbd_coordinates\n",
- "\n",
- "c = read_arbd_coordinates('sims/2-polymer/output/1-linear-fjc.restart')\n",
- "model.update_splines(c)\n",
- "\n",
- "model.monomers_per_bead_group = 1.0 ;# really per bead in this case\n",
- "model.generate_beads()\n",
- "\n",
- "engine.simulate(model, output_name=\"2-linear-fjc-hres\", directory='sims/2-polymer',\n",
- " num_steps=1e4, gpu=1\n",
- ")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "2609f481",
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "u = mda.Universe(\"sims/2-polymer/2-linear-fjc-hres.psf\", \"sims/2-polymer/output/2-linear-fjc-hres.dcd\")\n",
- "w = nv.show_mdanalysis(u, multipleBonds=False)\n",
- "w.clear_representations()\n",
- "w.add_ball_and_stick(\"all\",radius=1, multipleBonds=False)\n",
- "w.center()\n",
- "w"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "3df91fb0",
- "metadata": {},
- "source": [
- "Let's reuse our polymer object with a better model of a disordered peptide, the hydrophobicity scale model"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "9c7671a6",
- "metadata": {},
- "outputs": [],
- "source": [
- "print(f'Number of points used to represent the polymer path before update: {len(polymer.position_spline_params[0][0])}')\n",
- "c = read_arbd_coordinates('sims/2-polymer/output/2-linear-fjc-hres.restart')\n",
- "model.update_splines(c)\n",
- "print(f'Number of points used to represent the polymer path after update: {len(polymer.position_spline_params[0][0])}')\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "527cd03f",
- "metadata": {},
- "outputs": [],
- "source": [
- "from arbdmodel.hps_polymer_model import HpsModel\n",
- "\n",
- "## This model is strictly 1 bead per amino acid/monomer\n",
- "seq = 'A'*polymer.num_monomers ;# poly(alanine)\n",
- "hps_model = HpsModel([polymer],sequences = [seq])\n",
- "engine.simulate(hps_model, output_name=\"3-linear-hps\", directory='sims/2-polymer',\n",
- " num_steps=1e4, gpu=1\n",
- ")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "db639806",
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "u = mda.Universe(\"sims/2-polymer/3-linear-hps.psf\", \"sims/2-polymer/output/3-linear-hps.dcd\")\n",
- "w = nv.show_mdanalysis(u)\n",
- "w.clear_representations()\n",
- "w.add_ball_and_stick(\"all\",radius=2)\n",
- "w.center()\n",
- "w"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "354401cf",
- "metadata": {
- "scrolled": false
- },
- "outputs": [],
- "source": [
- "from arbdmodel.kh_polymer_model import KhModel\n",
- "from arbdmodel.onck_polymer_model import OnckModel\n",
- "\n",
- "kh_model = KhModel([polymer],sequences = [seq])\n",
- "onck_model= OnckModel([polymer],sequences = [seq])\n",
- "\n",
- "## It's possible to merge models, but after this it will not be possible to update splines any longer\n",
- "combined_model = model\n",
- "for i,m in enumerate((hps_model, kh_model, onck_model)):\n",
- " m.translate( ((i+1)*50, 0, 0) )\n",
- " combined_model.extend(m)\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "cf6d1af3",
- "metadata": {},
- "outputs": [],
- "source": [
- "groups = combined_model.children[:1] + [_g.children[0] for _g in combined_model.children[1:]]\n",
- "for g in groups:\n",
- " print(g, len(g.children))\n",
- " try:\n",
- " ngrids = len(g.children[0].type_.grid_potentials)\n",
- " except:\n",
- " ngrids = 0\n",
- " print(f'{ngrids} grids attached to {g.children[0]}')\n",
- " \n",
- "## Note that these models all use different particle types and no interactions are specified between them\n",
- "## If you want the polymers to \"feel\" each other, then you would need to add some custom interaction that \n",
- "## works on all the types using `add_nonbonded_interaction(interaction, types=(typeA,typeB))`\n",
- "## See the source files for the various models you've imported for examples of how to do this.\n",
- "\n",
- "engine.simulate(combined_model, output_name=\"4-linear-combined\", directory='sims/2-polymer',\n",
- " num_steps=1e4, gpu=1\n",
- ")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "8c965da5",
- "metadata": {},
- "outputs": [],
- "source": [
- "u = mda.Universe(\"sims/2-polymer/4-linear-combined.psf\", \"sims/2-polymer/output/4-linear-combined.dcd\")\n",
- "w = nv.show_mdanalysis(u)\n",
- "w.clear_representations()\n",
- "w.add_ball_and_stick(\"all\",radius=2)\n",
- "w.center()\n",
- "w"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8ee65f80",
- "metadata": {},
- "source": [
- "### Some polymer models may have multiple beads per monomer\n",
- "\n",
- "<img style=\"float: right; pad=10px;\" src=https://pubs.acs.org/cms/10.1021/ct500193u/asset/images/large/ct-2014-00193u_0001.jpeg width=300>\n",
- "\n",
- "One such model is a poly(dT) ssDNA model we built by optimizing the forces to match the structure observed in all-atom MD simulations, then adjusting the non-bonded forces to swell a DNA chain until the model matched experimental measurements (based on SAXS data) of the DNA size.\n",
- "https://pubs.acs.org/doi/10.1021/ct500193u\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "6212ec38",
- "metadata": {},
- "outputs": [],
- "source": [
- "# from importlib import reimport\n",
- "from arbdmodel.ssdna_two_bead import DnaModel\n",
- "# reimport(arbdmodel.ssdna_two_bead)\n",
- "\n",
- "## This model is strictly 2 bead per nt/monomer\n",
- "dna_model = DnaModel([polymer],sequences = None)\n",
- "engine.simulate(dna_model, output_name=\"5-linear-ssdna\", directory='sims/2-polymer',\n",
- " num_steps=1e4, gpu=1\n",
- ")\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "5b98ccfe",
- "metadata": {},
- "outputs": [],
- "source": [
- "u = mda.Universe(\"sims/2-polymer/5-linear-ssdna.psf\", \"sims/2-polymer/output/5-linear-ssdna.dcd\")\n",
- "w = nv.show_mdanalysis(u)\n",
- "w.clear_representations()\n",
- "w.add_ball_and_stick(\"all\",radius=1.0)\n",
- "w.center()\n",
- "w"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "1978d4a1",
- "metadata": {},
- "source": [
- "### Add a constant-force to the ends of the polymers to determine their elastic properties\n",
- "For the time being, all interactions in ARBD are tabulated, and forces are determined from potentials using linear interpolation. This is true for 1D and 3D potentials. It makes adding a constant force a _bit_ more tedious than it could be, but the grid-speficied potentials offer a versatile mechanism for importing data from experiments (e.g. cryo-EM) and continuum models."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "69f23403",
- "metadata": {},
- "outputs": [],
- "source": [
- "num_voxels = (4,4,4)\n",
- "\n",
- "## Size of each voxel along the three axes\n",
- "delta = [l/n for l,n in zip(combined_model.dimensions, num_voxels)]\n",
- "\n",
- "## Lower left corner of grid\n",
- "origin = combined_model.origin\n",
- "if origin is None: origin = -0.5 * np.array(combined_model.dimensions)\n",
- "\n",
- "## Center of voxels along each axis\n",
- "x,y,z = [(np.arange(n)+0.5)*d+o for n,d,o in zip(num_voxels,delta,origin)]\n",
- "print(f'x,y,z are 1D: x.shape = {x.shape}')\n",
- "print(f'x = {x}')\n",
- "\n",
- "## 3D arrays of cartesian coordinates\n",
- "X,Y,Z = np.meshgrid(x,y,z, indexing='ij')\n",
- "print(f'X,Y,Z are 3D (X.shape = {X.shape}) and can be used to compute potentials')\n",
- "print(f'X = \\n{X}\\n')\n",
- "# print(f'Y = \\n{Y}')\n",
- "print(f'X[:,:,0] = \\n{X[:,:,0]}')\n",
- "print(f'Y[:,:,0] = \\n{Y[:,:,0]}')\n",
- "\n",
- "## We can use these to compute other things, like the cylindrical or spherical radius:\n",
- "R = np.sqrt( X**2 + Y**2 + Z**2 )\n",
- "\n",
- "## units \"pN\" \"kcal_mol/AA\"\n",
- "U = 0.014393265 * Z;"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "84202784",
- "metadata": {},
- "outputs": [],
- "source": [
- "## The above technique is really flexible, which is why we guide you through it above, \n",
- "## but for simple cases, you might be able to use one of the functions packaged \n",
- "## with arbdmodel\n",
- "from arbdmodel.grid import slab_potential_z, constant_force, spherical_confinement\n",
- "U_pkg = constant_force( force=(0,0,0.014393265), dimensions=combined_model.dimensions, resolution = delta )\n",
- "assert( np.all(U == U_pkg) ) ## Check that we get the same thing!\n",
- "\n",
- "## Finally, we need to write the grid to disk to pass along to arbd\n",
- "from arbdmodel.grid import writeDx\n",
- "grid_filename = 'const_force_1pN_z.dx'\n",
- "writeDx(grid_filename, U, delta=delta, origin=origin)\n",
- "\n",
- "## It's a good idea to visualize the grid using VMD!"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "d05167cb",
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "## Let's add the DNA model for fun\n",
- "dna_model.translate( ((len(combined_model.children)+1)*50,0,0) )\n",
- "combined_model.extend(dna_model) "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "0b0a006e",
- "metadata": {},
- "outputs": [],
- "source": [
- "## Let's also update the segnames\n",
- "groups = combined_model.children[:1] + [_g.children[0] for _g in combined_model.children[1:]]\n",
- "for i,g in enumerate(groups):\n",
- " segname = '{:c}PRO'.format(i+65) # APRO,BPRO,CPRO,...\n",
- " g.segname = segname\n",
- "groups[-1].segname = 'ADNA'\n",
- "\n",
- "## And verify that it worked\n",
- "for g in groups[:-1]:\n",
- " print(g.children[0].segname)\n",
- "\n",
- "## The DNA model has a group for each nucleotide containing two beads, so it's a bit more complicated\n",
- "g = groups[-1] # Strand group associated with DNA model\n",
- "print(g.children[0].children[0].segname)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "a3c707b3",
- "metadata": {},
- "outputs": [],
- "source": [
- "## Now let's stretch our polymers by a 10pN force, first selectign the groups\n",
- "groups = combined_model.children[:1] + [_g.children[0] for _g in combined_model.children[1:]]\n",
- "for g in groups:\n",
- " print(g, len(g.children))\n",
- " try:\n",
- " ngrids = len(g.children[0].type_.grid_potentials)\n",
- " except:\n",
- " ngrids = 0\n",
- " print(f'{ngrids} grids attached to {g.children[0]}')\n",
- "\n",
- "## The 2-bead/nt DNA model is a bit different from the others\n",
- "try:\n",
- " ngrids = len(groups[-1].children[0].children[0].type_.grid_potentials)\n",
- "except:\n",
- " ngrids = 0\n",
- "print(f'{ngrids} grids attached to {groups[-1].children[0].children[0]}')\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "cf9d35b2",
- "metadata": {},
- "outputs": [],
- "source": [
- "## The thermostat can also cause some drift, so let's suppress that by creating a \"group site\" for each segment\n",
- "for g in groups:\n",
- " particles = [p for p in g]\n",
- " particles = particles[::10] ;# Skip every 10 to make the restraints a little more efficient\n",
- " site = combined_model.add_group_site(particles)\n",
- " ## 'site' can be used like any regular particle with grid potentials, bonds, angles, etc\n",
- " center = g.get_center()\n",
- " print(f'Restraining segment {g.segname} to center {center}')\n",
- " site.add_restraint( (10,center) ) # 10 kcal/mol/AA**2 spring constant"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "48055ae6",
- "metadata": {},
- "outputs": [],
- "source": [
- "for g in groups[:-1]:\n",
- " g.children[0].add_grid_potential(grid_filename, scale=-10)\n",
- " g.children[-1].add_grid_potential(grid_filename, scale=10)\n",
- "\n",
- "## The DNA model has a group for each nucleotide containing two beads, so it's a bit more complicated\n",
- "g = groups[-1] # Strand group associated with DNA model\n",
- "g.children[0].children[0].add_grid_potential(grid_filename, scale=-10)\n",
- "g.children[-1].children[0].add_grid_potential(grid_filename, scale=10)\n",
- " \n",
- "## If you made a mistake, you might need to fix the underlying data structures\n",
- "\"\"\"\n",
- "for g in groups[:-1]:\n",
- " for b in [g.children[i] for i in [0,-1]]:\n",
- " t0 = b.type_.parent\n",
- " assert(t0 is not None)\n",
- " b.type_ = t0 ## This will return the particle type to it's original value, clearing any applied grids\n",
- "\"\"\";"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f1400aee",
- "metadata": {},
- "outputs": [],
- "source": [
- "## Double-check that we added the grid potentials\n",
- "for g in groups[:]:\n",
- " try:\n",
- " ngrids = len(g.children[0].type_.grid_potentials)\n",
- " except:\n",
- " ngrids = 0\n",
- " print(f'{ngrids} grids attached to {g.children[0]}')\n",
- "\n",
- "## The 2-bead/nt DNA model is a bit different from the others\n",
- "try:\n",
- " ngrids = len(groups[-1].children[0].children[0].type_.grid_potentials)\n",
- "except:\n",
- " ngrids = 0\n",
- "print(f'{ngrids} grids attached to {groups[-1].children[0].children[0]}')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "c5492d02",
- "metadata": {},
- "outputs": [],
- "source": [
- "engine.simulate(combined_model, output_name=\"6-linear-combined-10pN\", directory='sims/2-polymer',\n",
- " num_steps=1e5, gpu=1, output_period=1e3\n",
- ")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "45fe39d9",
- "metadata": {},
- "outputs": [],
- "source": [
- "u = mda.Universe(\"sims/2-polymer/6-linear-combined-10pN.psf\",\n",
- " \"sims/2-polymer/6-linear-combined-10pN.pdb\",\n",
- " \"sims/2-polymer/output/6-linear-combined-10pN.dcd\")\n",
- "w = nv.show_mdanalysis(u)\n",
- "w.clear_representations()\n",
- "w.add_ball_and_stick(\"all\",radius=2, color_scheme='atomindex')\n",
- "w.center()\n",
- "w"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "45e13f9c",
- "metadata": {},
- "source": [
- "Notice that the models have obviously distinct relaxation timescales, primarily influenced by the damping coefficient prescribed by each model.\n",
- "\n",
- "Now let's plot the extension. Are the simulations long enough?"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "5dcd3c3e",
- "metadata": {},
- "outputs": [],
- "source": [
- "%matplotlib inline\n",
- "from matplotlib import pyplot as plt\n",
- "\n",
- "\n",
- "def get_extension(universe, segment_index, stride=1):\n",
- " seg = u.segments[segment_index].atoms\n",
- " data = np.empty( len(u.trajectory[::stride]) )\n",
- " for i,t in enumerate(u.trajectory[::stride]):\n",
- " r = seg.positions\n",
- " data[i] = r[-1,2] - r[0,2]\n",
- " return data\n",
- "\n",
- "for i,s in enumerate(u.segments):\n",
- " x = get_extension(u,i)\n",
- " plt.plot(x,label=s.segid)\n",
- " \n",
- "plt.ylabel('Extension (\\N{ANGSTROM SIGN})')\n",
- "plt.xlabel('Frame')\n",
- "plt.legend()\n",
- "plt.show()"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "658ec3df",
- "metadata": {},
- "source": [
- "That's all we mean to cover for now. \n",
- "\n",
- "We'll end by noting that the Polymer class provides some data structures for building branching polymers as well as linear ones covered in the tutorial. For more information about branching polymers, see the mrdna tutorial (for the time being)."
- ]
- }
- ],
- "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.8.18"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/.ipynb_checkpoints/3-run-and-customized-checkpoint.ipynb b/.ipynb_checkpoints/3-run-and-customized-checkpoint.ipynb
deleted file mode 100644
index 670aa50f2fa011b6afc1cfb30a314a3d18b98796..0000000000000000000000000000000000000000
--- a/.ipynb_checkpoints/3-run-and-customized-checkpoint.ipynb
+++ /dev/null
@@ -1,168 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "id": "024005b4",
- "metadata": {},
- "source": [
- "# Learn to run simulation and implement customized coarsed-grained polymer models\n",
- "\n",
- "As we've dived into the polymer objects in arbd and learned the basic usage of arbd, now it's time to put everything together and run coarsed-grained protein model using arbd! \n",
- "We will constructing a polymer model by coarse-graining the HpsModel.\n",
- "\n",
- "## Step 1: Model Construction\n",
- "\n",
- "We will start with the Hps Model we just saw in section 2 following and run the initial conformation with the same procedures. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "04d59bef",
- "metadata": {},
- "outputs": [],
- "source": [
- "import numpy as np\n",
- "from pathlib import Path\n",
- "from arbdmodel import ArbdEngine\n",
- "from arbdmodel.coords import readArbdCoords\n",
- "from arbdmodel.polymer import PolymerSection\n",
- "from arbdmodel.hps_polymer_model import _types\n",
- "from arbdmodel.hps_polymer_model import HpsModel as NupModel\n",
- "\n",
- "#from info import NupModel, create_arbd_polymer_objects, dimensions\n",
- "\n",
- "temperature = 298.15\n",
- "ion_concentration = 250 # in mM\n",
- "gpu = 0\n",
- "\n",
- "decomp_period = 50\n",
- "skin_depth = 8\n",
- "## 3.2\n",
- "\n",
- "model_name = 'HPS'\n",
- "step = 1\n",
- "\n",
- "## units \"sqrt(80 epsilon0 k K /(2 * (mM/particle) * e**2))\" AA\n",
- "debye_length = 10\n",
- "\n",
- "polymers, sequences = create_arbd_polymer_objects()\n",
- "\n",
- "model = NupModel( polymers, sequences,\n",
- " debye_length = debye_length,\n",
- " temperature=temperature,\n",
- " damping_coefficient = 50000, # units of 1/ns\n",
- " decomp_period = decomp_period,\n",
- " pairlist_distance = 50+skin_depth,\n",
- " dimensions = dimensions\n",
- " )\n",
- "engine = ArbdEngine( integrator = 'Langevin',\n",
- " num_steps=1e7,\n",
- " output_period=1e4,\n",
- " gpu = 0\n",
- " )\n",
- "directory = f'1-one_bead_per_res-{model_name}'\n",
- "outname = 'run'\n",
- "\n",
- "if not Path(f'{directory}/output/{outname}.dcd').exists():\n",
- " # \"\"\" \"Minimization\" \"\"\"\n",
- " # add_restraints(model,restraints)\n",
- " # model.simulate(output_name='{}-min-{}'.format(model_name,step),\n",
- " # num_steps=1e3, output_period=500,\n",
- " # gpu=gpu,\n",
- " # )\n",
- " # coords = readArbdCoords('output/{}-min-{}.restart'.format(model_name,step))\n",
- " # model.update_splines(coords)\n",
- "\n",
- " \"\"\" Production \"\"\"\n",
- " model.set_damping_coefficient( 100 )\n",
- " engine.simulate(model, output_name=outname,\n",
- " directory = directory\n",
- " )\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "5fb80332",
- "metadata": {},
- "source": [
- "Then we can run the simulation with "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "28a28e18",
- "metadata": {},
- "outputs": [],
- "source": [
- "\n",
- "\n",
- "_seq = 'GLFG' * 8 # 10 repeats per polymer\n",
- "## We're going to tile polymers along x and y in a square lattice\n",
- "Nx = 5\n",
- "Ny = 5\n",
- "N_polymers = Nx*Ny\n",
- "\n",
- "density = 55 # mg / mL\n",
- "\n",
- "\"\"\" Calculated parameters \"\"\"\n",
- "polymer_mass = sum( [_types[aa].mass for aa in _seq] ) # in daltons\n",
- "## units \"mg/ml\" \"dalton/AA**3\"\n",
- "_conversion = 0.0006022142\n",
- "\n",
- "dimensions = (Nx*Ny*polymer_mass/(density*_conversion))**(1/3) # in Angstroms, the unit of length in ARBD\n",
- "dimensions = [dimensions]*3 # along x,y,z\n",
- "\n",
- "xs = np.linspace( -dimensions[0]*0.5, dimensions[0]*0.5, Nx+1 ) # row edges (one extra value)\n",
- "ys = np.linspace( -dimensions[1]*0.5, dimensions[1]*0.5, Ny+1 ) # column edges\n",
- "xs = (xs[1:] + xs[:-1])*0.5 # row centers\n",
- "ys = (ys[1:] + ys[:-1])*0.5 # column centers\n",
- "\n",
- "z = np.arange(len(_seq))*3.8*0.5 # times coordinate for every amino acid in a polymer, compressed a bit\n",
- "\n",
- "\"\"\" Build peptide list consisting of sequence, coordinates\"\"\"\n",
- "peptides = []\n",
- "\n",
- "for x in xs:\n",
- " for y in ys:\n",
- " r = np.empty( (len(_seq),3) ) # allocate array for coordinates of each amino acid in polymer\n",
- " r[:,0] = x\n",
- " r[:,1] = y\n",
- " r[:,2] = z\n",
- "\n",
- " peptides.append( (_seq, r ) )\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "b5aef887",
- "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.8.18"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}