In [1]:
#initialize.http://127.0.0.1:8889/notebooks/Force_Analyze.ipynb#
%autosave 120
%pylab inline
%qtconsole
In [2]:
#change working directory
import os
working_directory = "/var/run/media/chen/Simulation/Simulation/100mM_1000f_vs_100f"
os.chdir(working_directory)
In [3]:
#Other initialize
import pandas as pd
import numpy as np
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: