Skip to content
Snippets Groups Projects
Commit 8e2c111a authored by cmaffeo2's avatar cmaffeo2
Browse files

Update 3 files 

- /.ipynb_checkpoints/1-basics-checkpoint.ipynb
- /.ipynb_checkpoints/2-polymer-objects-checkpoint.ipynb
- /.ipynb_checkpoints/3-run-and-customized-checkpoint.ipynb
parent dbc8129d
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
.ipynb_checkpoints/1-basics-checkpoint.ipynb deleted 100644 → 0
+ 0
− 311
View file @ dbc8129d
%% Cell type:markdown id: tags:
# Basic usage of the `arbdmodel` package
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.
This tutorial assumes you have a working knowledge of Python, including its object-oriented features, and the numpy package.
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 id: tags:
### Step 1: Create particle types
%% Cell type:code id: tags:
``` python
import ipdb # in case we need to debug, after a crash, run `ipdb.pm()`
import numpy as np
from arbdmodel import ParticleType, PointParticle, ArbdModel
system_size = 620
num_particles = 6400
## diffusion of liquid argon: https://www.sciencedirect.com/science/article/pii/0375960171901290
## units "60 * 3.41 AA * sqrt(119.5 k K / 40 dalton)" "AA**2/ns"
argon = ParticleType(name="Ar",
mass=39.95, # dalton
damping_coefficient=5000, # per ns
custom_property="my_value",
epsilon=0.177, # kcal_mol
radius=3.345/2) # Angstroms
print(argon)
```
%% Cell type:markdown id: tags:
### Step 2: Build a system
%% Cell type:code id: tags:
``` python
## Generate Nx3 array of random cartesian coordinates
positions = system_size*(np.random.random([num_particles,3])-0.5)
## Create a list of point particles located at those positions
points = [PointParticle(type_=argon, name="Ar", position=positions[i])
for i in range(num_particles)]
model = ArbdModel(points,
dimensions=[system_size for i in range(3)], # Ångstroms
timestep = 20e-6, # ns; can be specified below with engine instead
)
print(model)
print(model.children[:2])
```
%% Cell type:markdown id: tags:
### Step 3: Describe the interactions between the particles
%% Cell type:code id: tags:
``` python
from arbdmodel.interactions import LennardJones
lj = LennardJones()
model.add_nonbonded_interaction(lj)
```
%% Cell type:markdown id: tags:
### Step 4: Run the simulation
%% Cell type:code id: tags:
``` python
from arbdmodel import ArbdEngine
engine = ArbdEngine(output_period = 1e2, num_steps = 1e5)
engine.simulate(model, output_name="1-lj", directory='sims/1-argon',
num_steps=1e3, gpu=1
) # simulation parameters here override those in constructor
```
%% Cell type:markdown id: tags:
### Step 5: Visualize the results
We recommend using [VMD](https://www.ks.uiuc.edu/Development/Download/download.cgi?PackageName=VMD) for visualizing the simulations trajectories with:
```bash
vmd sims/1-argon/1-lj.psf sims/1-argon/output/1-lj.dcd
```
However, for convenience we have provided commands to use the nglview module, which provides browser-based visualization.
note: you must enable certain jupyter extensions to use nglview as follows before running jupyter-notebook
```bash
jupyter-nbextension enable --py --sys-prefix widgetsnbextension
jupyter-nbextension enable nglview --py --sys-prefix
```
%% Cell type:code id: tags:
``` python
import MDAnalysis as mda
import nglview as nv
u = mda.Universe("sims/1-argon/1-lj.pdb", "sims/1-argon/output/1-lj.dcd")
w = nv.show_mdanalysis(u)
w.clear_representations()
w.add_ball_and_stick("all",radius=10)
w
```
%% Cell type:markdown id: tags:
### Step 6: Customize the interactions
%% Cell type:code id: tags:
``` python
from arbdmodel.interactions import AbstractPotential
## Clear nonbonded interactions
model.nonbonded_interactions = []
class BuckinghamPotential(AbstractPotential):
def potential(self, r, types):
## https://royalsocietypublishing.org/doi/10.1098/rspa.1938.0173
## optionally could pull information from typeA, typeB
typeA, typeB = types
return 143932.65 * ( 1.69 * np.exp( -r/0.273 ) - 102e-4 / r**6 )
pot = BuckinghamPotential()
model.add_nonbonded_interaction( pot )
engine.simulate(model, output_name="2-buck", directory='sims/1-argon',
num_steps = 1e3, gpu=1
)
```
%% Cell type:code id: tags:
``` python
import MDAnalysis as mda
import nglview as nv
u = mda.Universe("sims/1-argon/2-buck.pdb", "sims/1-argon/output/2-buck.dcd")
w = nv.show_mdanalysis(u)
w.clear_representations()
w.add_ball_and_stick("all",radius=10)
w
```
%% Cell type:markdown id: tags:
### Step 7: Add bonds
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 id: tags:
``` python
from arbdmodel.interactions import HarmonicBond
bond = HarmonicBond(k=100, # kcal/mol AA**2
r0 = 3 # AA
)
## Go through every other index
for i in range(0,len(model.children),2):
j = i+1
model.add_bond( model.children[i],
model.children[j],
bond )
print("numbonds:", len(model.bonds))
engine.simulate(model, output_name="3-bonds", directory='sims/1-argon',
num_steps = 1e3, gpu=1
)
```
%% Cell type:code id: tags:
``` python
import MDAnalysis as mda
import nglview as nv
u = mda.Universe("sims/1-argon/3-bonds.psf","sims/1-argon/3-bonds.pdb", "sims/1-argon/output/3-bonds.dcd")
w = nv.show_mdanalysis(u)
w.clear_representations()
w.add_ball_and_stick("index < 10",radius=10)
w
```
%% Cell type:code id: tags:
``` python
## Why are there so many long bonds?
print(bond.range_) ## Distances are in Angstroms
r = np.linspace(*bond.range_, 20)
print(r)
print(bond.potential(r)) ## Energies are kcal/mol
## Outside of this range, the bond applies zero force
```
%% Cell type:markdown id: tags:
### ARBD is configured through plain text files
Use your favorite editor to examine the files named `sims/1-argon/*.bd`.
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.
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.
.ipynb_checkpoints/2-polymer-objects-checkpoint.ipynb deleted 100644 → 0
+ 0
− 597
View file @ dbc8129d
This diff is collapsed.
.ipynb_checkpoints/3-run-and-customized-checkpoint.ipynb deleted 100644 → 0
+ 0
− 168
View file @ dbc8129d
%% Cell type:markdown id:024005b4 tags:
# Learn to run simulation and implement customized coarsed-grained polymer models
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!
We will constructing a polymer model by coarse-graining the HpsModel.
## Step 1: Model Construction
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 id:04d59bef tags:
``` python
import numpy as np
from pathlib import Path
from arbdmodel import ArbdEngine
from arbdmodel.coords import readArbdCoords
from arbdmodel.polymer import PolymerSection
from arbdmodel.hps_polymer_model import _types
from arbdmodel.hps_polymer_model import HpsModel as NupModel
#from info import NupModel, create_arbd_polymer_objects, dimensions
temperature = 298.15
ion_concentration = 250 # in mM
gpu = 0
decomp_period = 50
skin_depth = 8
## 3.2
model_name = 'HPS'
step = 1
## units "sqrt(80 epsilon0 k K /(2 * (mM/particle) * e**2))" AA
debye_length = 10
polymers, sequences = create_arbd_polymer_objects()
model = NupModel( polymers, sequences,
debye_length = debye_length,
temperature=temperature,
damping_coefficient = 50000, # units of 1/ns
decomp_period = decomp_period,
pairlist_distance = 50+skin_depth,
dimensions = dimensions
)
engine = ArbdEngine( integrator = 'Langevin',
num_steps=1e7,
output_period=1e4,
gpu = 0
)
directory = f'1-one_bead_per_res-{model_name}'
outname = 'run'
if not Path(f'{directory}/output/{outname}.dcd').exists():
# """ "Minimization" """
# add_restraints(model,restraints)
# model.simulate(output_name='{}-min-{}'.format(model_name,step),
# num_steps=1e3, output_period=500,
# gpu=gpu,
# )
# coords = readArbdCoords('output/{}-min-{}.restart'.format(model_name,step))
# model.update_splines(coords)
""" Production """
model.set_damping_coefficient( 100 )
engine.simulate(model, output_name=outname,
directory = directory
)
```
%% Cell type:markdown id:5fb80332 tags:
Then we can run the simulation with
%% Cell type:code id:28a28e18 tags:
``` python
_seq = 'GLFG' * 8 # 10 repeats per polymer
## We're going to tile polymers along x and y in a square lattice
Nx = 5
Ny = 5
N_polymers = Nx*Ny
density = 55 # mg / mL
""" Calculated parameters """
polymer_mass = sum( [_types[aa].mass for aa in _seq] ) # in daltons
## units "mg/ml" "dalton/AA**3"
_conversion = 0.0006022142
dimensions = (Nx*Ny*polymer_mass/(density*_conversion))**(1/3) # in Angstroms, the unit of length in ARBD
dimensions = [dimensions]*3 # along x,y,z
xs = np.linspace( -dimensions[0]*0.5, dimensions[0]*0.5, Nx+1 ) # row edges (one extra value)
ys = np.linspace( -dimensions[1]*0.5, dimensions[1]*0.5, Ny+1 ) # column edges
xs = (xs[1:] + xs[:-1])*0.5 # row centers
ys = (ys[1:] + ys[:-1])*0.5 # column centers
z = np.arange(len(_seq))*3.8*0.5 # times coordinate for every amino acid in a polymer, compressed a bit
""" Build peptide list consisting of sequence, coordinates"""
peptides = []
for x in xs:
for y in ys:
r = np.empty( (len(_seq),3) ) # allocate array for coordinates of each amino acid in polymer
r[:,0] = x
r[:,1] = y
r[:,2] = z
peptides.append( (_seq, r ) )
```
%% Cell type:code id:b5aef887 tags:
``` python
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment