In [19]:
#initialize.http://127.0.0.1:8889/notebooks/Force_Analyze.ipynb#
%autosave 120
%pylab inline
%qtconsole
In [1]:
#Other initilize
from IPython.display import HTML, Image
import os
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 [3]:
Image(url="https://raw.github.com/maverickmw/pfdreport/master/movies/initial_perspect.gif")
Out[3]:
In [6]:
#Typical Flow Chart of MD Simulation steps in Gromacs.
Image("gromacs_flow_chart.png")
Out[6]:
In [4]:
#The nucleosome that we familiar with
Image(url = "http://upload.wikimedia.org/wikipedia/commons/d/d9/Nucleosome_1KX5_colour_coded.png", width=300, height=300)
Out[4]:
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 [10]:
#VMD as Conformation
Image("../pfdreport/conformation.png")
Out[10]:
In [9]:
#VMD as Fragment
Image("../pfdreport/fragment.png")
Out[9]:
In [11]:
#Expected~
Image(url="http://www.accessexcellence.org/RC/VL/GG/images/nucleosome.gif")
Out[11]:
In [8]:
# VMD resid
Image("../pfdreport/resid.png")
Out[8]:
In [11]:
#Nucleosome Boom!
Image(url="https://raw.github.com/maverickmw/pfdreport/master/movies/nucleosome_openmm_100ps_gpu_vacuum.gif")
Out[11]:
In [16]:
#Bonded force
Image("bonded_force.png")
Out[16]:
While (simulation loop):
do_force() //All on Atom Level
{
do_ns(); //cost :$O(N)$
do_force_lowlevel();
{
do_nonbonded(); //cost :$O(N^2)$ brute forcely; $O(N)$ after ns()
//calc bonds
}
//anything else;
}
In [16]:
#pull 1300ps
Image(url="https://raw.github.com/maverickmw/pfdreport/master/movies/pull_1300.gif")
Out[16]:
In [17]:
#relax 6000
Image(url="https://raw.github.com/maverickmw/pfdreport/master/movies/relax_6000.gif")
Out[17]:
In [20]:
#residue range:
H3_1 = arange(1,136) #0
H3_2 = arange(488, 623) #4
H4_1 = arange(136,238) #1
H4_2 = arange(623,725) #5
H2A_1 = arange(238,366) #2
H2A_2 = arange(725,853) #6
H2B_1 = arange(366,488) #3
H2B_2 = arange(853,975) #7
DNA = arange(975,1307)
DNA_1 = arange(975,1141) #8
DNA_2 = arange(1141,1307) #9
#print DNA_2
In [22]:
f_in = pd.read_csv("../pfdreport/nucleosome_30nm_md300_1300ps_novs_skip_10/force_external.csv")
f_bond = pd.read_csv("../pfdreport/nucleosome_30nm_md300_1300ps_novs_skip_10/force_bonded.csv")
In [25]:
cl=f_in.columns
print cl
In [27]:
dt_63 = f_in[cl[63]]
dt_83 = f_in[cl[83]]
time = f_in[cl[0]]
In [30]:
plot(time,dt_63)
plot(time,dt_83)
legend([cl[63],cl[83]])
show()
In [32]:
Image("time_scale.png")
Out[32]:
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 [5]:
#Image(url="http://ww3.sinaimg.cn/bmiddle/61e8a1fdjw1ebcpe8vqghg203z03zq3h.gif")
In [ ]: