In [4]:
#initialize.http://127.0.0.1:8889/notebooks/Force_Analyze.ipynb#
%autosave 120
%pylab inline
%qtconsole
In [2]:
#Other initilize
from IPython.display import HTML, Image
In [2]:
#change working directory
import os
working_directory = "/var/run/media/chen/Simulation/Simulation/100mM_1000f_vs_100f"
os.chdir(working_directory)
In [2]:
#Other initialize
import pandas as pd
import numpy as np
In [6]:
#Typical Flow Chart of MD Simulation steps in Gromacs.
Image("gromacs_flow_chart.png")
Out[6]:
In [8]:
#The nucleosome that we familiar with
Image(url = "http://upload.wikimedia.org/wikipedia/commons/d/d9/Nucleosome_1KX5_colour_coded.png")
Out[8]:
The crystal structure of the nucleosome core particle consisting of H2A , H2B , H3 and H4 core histones, and DNA. The view is from the top through the superhelical axis.
In [ ]:
#VMD as Conformation
In [8]:
#VMD as Fragment
In [11]:
#Expected~
Image(url="http://www.accessexcellence.org/RC/VL/GG/images/nucleosome.gif")
Out[11]:
In [12]:
# VMD resid
In [14]:
#Nucleosome Boom!
In [16]:
#Bonded force
Image("bonded_force.png")
Out[16]:
In [ ]:
Force is the main reason that drives the molecular dynamic simulation. It was calculated every simuatlion step, in the use of updating coordinate of simulation structure.
For a standard simulation procedure, from a certain pdb structure to simulation result, including the following steps:
Choose a certain Force Field, including the Atom discription(**mass, charge**), bonded force(**bonds, pairs for nb14, angles, and dihedrals X2**) and nonbonded force (**vdw and coulomb**)
Add water and salts.
for implicit solvent add other artificial effects
To avoid force between particle and particle that larger than 1000 force unit ($kJ\cdot mol^{-1}\cdot nm^{-1}$)
As for:
$1kJ = 1000 N\cdot m$
$1nm = 10 ^ {-9} m$
$1mol = 6.022\cdot 10 ^ {-23}$ uintless
So that, for 1 force unit:
while (steps < Set_to_simulation_steps)
{
dd_partition_system();
do_force();
write_traj();
update();
steps++;
}
inside the loop, the steps do_force() calculate the forces and update the trajectory:
do_force()
{
dd_move_x();
do_ns();
do_force_lowlevel();
dd_move_f();
calc_viral();
anything else;
}
do_force_lowlevel(): The core scheme to calculate force, include:
do_force_lowlevel()
{
do_walls(); //calc forces and viral for walls
do_nonbonded(); //calc forces and viral from nonbonded interaction
calc_bonds(); //calc bonds
}
Basic Consideration:
F(total) = F(bonded) + F(nonbonded)
For a single molecule, sum(F(bonded)) = 0; sum(F(nonbonded)) = 0
then: F(nonbonded) = F(total) - F(bonded)
1 If only need to know the force that loaded on molecules (globaly), then only bonded force is needed to be calculated out.
2 Pairwise force??
ToDo: explain detail or redifine //As: F(total) = F(internal) + F(external)
sum(F(internal)) == 0
F(external) = ??(try: Should over write ns()? //or use an OO search type) ==> modified nblists afterwardsInclude:
Bonded (cals_bonds(): run sequencially, not so difficult)
As the calculation is atomic based, first calculate the bonded force for each atom, then sequncially binned to resiedue.
(To check: A-B-C F(A) + F(C) == F(B)? )
Nonbonded (do_nonbonded.c and in kernel. Different kernels, merged, seperatable?)
elec (short and long, no PME)
vdw (should exclude nb internal. How?) //==> use two nblist to store internal and external nblists.
Method:
If only bonded force included, the world will be happy. Simply used the build in function clac_bonds.
To eveluate the nonbonded force, modification of kernel is not possible.
Improvement will be using the ns() (with different cut off?)
However, severial problems will occurs:
1 nb14
Not included in nb_kernel, dumped into bounded. (mostly zeros)
2 internal v.s external? (for DNA internal should define as Double Helix.)
F = F(internal_bond) + F(internal_nonbonded) + F(external_nonbonded) ~ 0. (Can not be achieved since PME, water and ions are gone, the md external will never be the same with analyze external)
In [2]:
Pairwise force calculation, (only) support 1 slection or 2 selections (1 selection can be one molecule or some atoms, or some residues, whatever), and calculate forces acording to particle(atom) and residue map
* For 2 selections:
* Bonded_force (most time meaningless) //to check for the parameter setting
* Nonbonded_force
* nb14 //to check for the parameter settingLooping through index and sum F_coul and F_vdw
Can use the nbkernal All_vs_all
Rely on nblist (cutoff setting),
input:
output:
Principle:
Steps:
In [5]:
os.listdir(working_directory)
Out[5]:
In [9]:
fm = pd.read_csv('force_external_)
cl=fm.columns
In [1]:
dt = fm[cl[0]]
In [29]:
df = fm[cl[491]]
a = plot(dt, df)
In [31]:
max(df)
cl[491]
Out[31]:
The most difficult part!!!!!!!
The structure of nblist:
save each frame residue based data as csv format.
!!!Take care of the virtual site force
Data Structure of t_nblist:
t_excl* excl; //Exclustions, only with enltypeCG
void* kernelptr_vf;
int simd_padding_width;
In [ ]:
When do a pulling simulation, the frequency of unit that has been loaded force, is going to increase.
Consider a simple Harmonic Oscillator, with out external force or friction:
In [13]:
Image(url="http://ww3.sinaimg.cn/bmiddle/61e8a1fdjw1ebcpe8vqghg203z03zq3h.gif")
Out[13]:
In [ ]: