buildH_calcOP_test module
This script builds hydrogens from a united-atom trajectory and calculate the order parameter for each C-H bond.
- It works in two modes :
A slow mode when an output trajectory (e.g. in xtc format) is requested by the user. In this case, the whole trajectory including newly built hydrogens are written to this trajectory.
A fast mode without any output trajectory.
For both modes, the order parameter is written to an output file in a format similar to the code of @jmelcr: https://github.com/NMRLipids/MATCH/blob/master/scripts/calcOrderParameters.py
This code has been checked against the one from @jmelcr. You might find minor differences due to rounding errors (in xtc, only 3 digits are written).
The way of building H is largely inspired from a code of Jon Kapla originally written in fortran : https://github.com/kaplajon/trajman/blob/master/module_trajop.f90#L242.
Note: that all coordinates in this script are handled using numpy 1D-arrays of 3 elements, e.g. atom_coor = np.array((x, y, z)). Note2: sometimes numpy is slow on small arrays, thus we wrote a few “in-house” functions for vectorial operations (e.g. cross product).
- buildH_calcOP_test.calc_OP(C, H)[source]
Returns the Order Parameter of a CH bond (OP).
OP is calculated according to equation:
S = 1/2 * (3*cos(theta)^2 -1)
theta is the angle between CH bond and the z(vertical) axis: z ^ H | / |/ C
Inspired from a function written by @jmelcr.
- Parameters:
C (numpy 1D-array) – Coordinates of C atom.
H (numpy 1D-array) – Coordinates of H atom.
- Returns:
The normalized vector.
- Return type:
float
- buildH_calcOP_test.normalize(vec)[source]
Normalizes a vector.
- Parameters:
vec (numpy 1D-array)
- Returns:
The normalized vector.
- Return type:
numpy 1D-array
- buildH_calcOP_test.norm(vec)[source]
Returns the norm of a vector.
- Parameters:
vec (numpy 1D-array)
- Returns:
The magniture of the vector.
- Return type:
float
- buildH_calcOP_test.calc_angle(atom1, atom2, atom3)[source]
Calculates the valence angle between atom1, atom2 and atom3.
Note: atom2 is the central atom.
- Parameters:
atom1 (numpy 1D-array.)
atom2 (numpy 1D-array.)
atom3 (numpy 1D-array.)
- Returns:
The calculated angle in radians.
- Return type:
float
- buildH_calcOP_test.vec2quaternion(vec, theta)[source]
Translates a vector of 3 elements and angle theta to a quaternion.
- Parameters:
vec (numpy 1D-array) – Vector of the quaternion.
theta (float) – Angle of the quaternion in radian.
- Returns:
The full quaternion (4 elements).
- Return type:
numpy 1D-array
- buildH_calcOP_test.calc_rotation_matrix(quaternion)[source]
Translates a quaternion to a rotation matrix.
- Parameters:
quaternion (numpy 1D-array of 4 elements.)
- Returns:
The rotation matrix.
- Return type:
numpy 2D-array (dimension [3, 3])
- buildH_calcOP_test.apply_rotation(vec_to_rotate, rotation_axis, rad_angle)[source]
Rotates a vector around an axis by a given angle.
Note: the rotation axis is a vector of 3 elements.
- Parameters:
vec_to_rotate (numpy 1D-array)
rotation_axis (numpy 1D-array)
rad_angle (float)
- Returns:
The final rotated (normalized) vector.
- Return type:
numpy 1D-array
- buildH_calcOP_test.pandasdf2pdb(df)[source]
Returns a string in PDB format from a pandas dataframe.
- Parameters:
df (pandas dataframe with columns "atnum", "atname", "resname", "resnum",) – “x”, “y”, “z”
- Returns:
A string representing the PDB.
- Return type:
str
- buildH_calcOP_test.cross_product(A, B)[source]
Returns the cross product between vectors A & B.
Source: http://hyperphysics.phy-astr.gsu.edu/hbase/vvec.html. Note: on small vectors (i.e. of 3 elements), computing cross products
with this functions is faster than np.cross().
- Parameters:
A (numpy 1D-array) – A vector of 3 elements.
B (numpy 1D-array) – Another vector of 3 elements.
- Returns:
Cross product of A^B.
- Return type:
numpy 1D-array
- buildH_calcOP_test.get_CH2(atom, helper1, helper2)[source]
Reconstructs the 2 hydrogens of a sp3 carbon (methylene group).
- Parameters:
atom (numpy 1D-array) – Central atom on which we want to reconstruct hydrogens.
helper1 (numpy 1D-array) – Heavy atom before central atom.
helper2 (numpy 1D-array) – Heavy atom after central atom.
- Returns:
Coordinates of the two hydrogens: ([x_H1, y_H1, z_H1], [x_H2, y_H2, z_H2]).
- Return type:
tuple of numpy 1D-arrays
- buildH_calcOP_test.get_CH(atom, helper1, helper2, helper3)[source]
Reconstructs the unique hydrogen of a sp3 carbon.
- Parameters:
atom (numpy 1D-array) – Central atom on which we want to reconstruct the hydrogen.
helper1 (numpy 1D-array) – First neighbor of central atom.
helper2 (numpy 1D-array) – Second neighbor of central atom.
helper3 (numpy 1D-array) – Third neighbor of central atom.
- Returns:
Coordinates of the rebuilt hydrogen: ([x_H, y_H, z_H]).
- Return type:
numpy 1D-array
- buildH_calcOP_test.get_CH_double_bond(atom, helper1, helper2)[source]
Reconstructs the hydrogen of a sp2 carbon.
- Parameters:
atom (numpy 1D-array) – Central atom on which we want to reconstruct the hydrogen.
helper1 (numpy 1D-array) – Heavy atom before central atom.
helper2 (numpy 1D-array) – Heavy atom after central atom.
- Returns:
Coordinates of the rebuilt hydrogen: ([x_H, y_H, z_H]).
- Return type:
tuple of numpy 1D-arrays
- buildH_calcOP_test.get_CH3(atom, helper1, helper2)[source]
Reconstructs the 3 hydrogens of a sp3 carbon (methyl group).
- Parameters:
atom (numpy 1D-array) – Central atom on which we want to reconstruct hydrogens.
helper1 (numpy 1D-array) – Heavy atom before central atom.
helper2 (numpy 1D-array) – Heavy atom before helper1 (two atoms away from central atom).
- Returns:
Coordinates of the 3 hydrogens: ([x_H1, y_H1, z_H1], [x_H2, y_H2, z_H2], [x_H3, y_H3, z_H3]).
- Return type:
tuple of numpy 1D-arrays
- buildH_calcOP_test.buildHs_on_1C(atom, dic_lipid)[source]
Builds 1, 2 or 3 H on a given carbon.
This function is a wrapper which gathers the coordinates of the helpers and call the function that builds 1, 2 or 3 H.
The name of the helpers as well as the type of H to build are described in a dictionnary stored in dic_lipids.py.
- Parameters:
atom (MDAnalysis Atom instance) – This instance contains the carbon on which we want to build Hs.
dic_lipid (dictionnary) – Comes from dic_lipids.py. Contains carbon names and helper names needed for reconstructing hydrogens.
- Returns:
Each element of the tuple is a numpy 1D-array containing 1, 2 or 3 reconstructed hydrogen(s). !!! IMPORTANT !!! This function should return a tuple even if there’s only one H that has been rebuilt.
- Return type:
tuple of numpy 1D-arrays
- buildH_calcOP_test.build_all_Hs_calc_OP(universe_woH, dic_lipid, dic_Cname2Hnames, universe_wH=None, dic_OP=None, dic_corresp_numres_index_dic_OP=None, return_coors=False)[source]
Main function that builds all hydrogens and calculates order parameters.
This function shall be used in two modes :
1) The first time this function is called, we have to construct a new universe with hydrogens. One shall call it like this :
new_data_frame = build_all_Hs_calc_OP(universe_woH, return_coors=True)
The boolean return_coors set to True indicates to the function to return a pandas dataframe. This latter will be used later to build a new universe with H.
2) For all the other frames, we just need to update the coordinates in the universe with hydrogens. One shall call it like this :
- build_all_Hs_calc_OP(universe_woH, dic_lipid, dic_Cname2Hnames,
universe_wH=universe_wH, dic_OP=dic_OP, dic_corresp_numres_index_dic_OP=dic_corresp_numres_index_dic_OP)
In this case, the function also calculates the order parameter and returns nothing. The coordinates of the universe with H are updated in place. The order parameter is also added in place (within dic_OP dictionnary).
NOTE: This function in mode 2 is slow, thus it shall be used when one wants to create a trajectory with H (such as .xtc or whatever format).
NOTE2: This function assumes all possible C-H pairs are present in the .def file (with -d option). They are needed since we want to build an xtc with the whole system. If one is interested in calculating only a subset of OPs, please use the function fast_build_all_Hs_calc_OP() instead.
- Parameters:
universe_woH (MDAnalysis universe instance) – This is the universe without hydrogen.
dic_lipid (dictionnary) – Comes from dic_lipids.py. Contains carbon names and helper names needed for reconstructing hydrogens.
dic_Cname2Hnames (dictionnary) –
This dict gives the correspondance Cname -> Hname. It is a dict of tuples. If there is more than 1 H for a given C, they need to be ordered like in the PDB. e.g. for CHARMM POPC : {‘C13’: (‘H13A’, ‘H13B’, ‘H13C’), …, ‘C33’: (‘H3X’, ‘H3Y’),
…, ‘C216’: (‘H16R’, ‘H16S’), …}
universe_wH (MDAnalysis universe instance (optional)) – This is the universe with hydrogens.
dic_OP (ordered dictionnary) – Each key of this dict is a couple carbon/H, and at the beginning it contains an empty list, e.g. OrderedDict([ (‘C1’, ‘H11): [], (‘C1’, ‘H12’): [], … ]) See function init_dic_OP() below to see how it is organized.
dic_corresp_numres_index_dic_OP (dictionnary) – This dict should contain the correspondance between the numres and the corresponding index in dic_OP. For example {…, 15: 14, …} means the residue numbered 15 in the PDB has an index of 14 in dic_OP.
return_coors (boolean (optional)) – If True, the function will return a pandas dataframe containing the system with hydrogens.
- Returns:
pandas dataframe (optional) – If parameter return_coors is True, this dataframe contains the system with hydrogens is returned.
None – If parameter return_coors is False.
- buildH_calcOP_test.fast_buildHs_on_1C(dic_lipids_with_indexes, ts, Cname, ix_first_atom_res)[source]
Builds 1, 2 or 3 H on a given carbon using fast indexing.
This function is a fast wrapper which gathers the coordinates of the helpers and call the function that builds 1, 2 or 3 H.
The name of the helpers as well as the type of H to build are described in a dictionnary stored in dic_lipids.py.
- Parameters:
dic_lipids_with_indexes (dictionnary) – The dictionnary made in function make_dic_lipids_with_indexes().
ts (MDAnalysis Timestep instance) – This object contains the actual frame under analysis.
Cname (str) – The carbon name on which we want to build H(s).
ix_first_atom_res (int) – The index of the first atom in the lipid under analysis.
- Returns:
Each element of the tuple is a numpy 1D-array containing 1, 2 or 3 reconstructed hydrogen(s). !!! IMPORTANT !!! This function should return a tuple even if there’s only one H that has been rebuilt.
- Return type:
tuple of numpy 1D-arrays
- buildH_calcOP_test.get_indexes(atom, universe_woH, dic_lipid)[source]
Returns the index of helpers for a given carbon.
- Parameters:
atom (MDAnalysis Atom instance) – This is an Atom instance of a carbon on which we want to build Hs.
universe_woH (MDAnalysis universe instance) – The universe without hydrogens.
dic_lipid (dictionnary) – Comes from dic_lipids.py. Contains carbon names and helper names needed for reconstructing hydrogens.
- Returns:
The tuple contains the index of the 2 (or 3) helpers for the atom that was passed as argument. (e.g. for atom C37 with index 99, the function returns a tuple containing 98 (index of C36 = helper 1) and 100 (index of C38=helper2).
- Return type:
tuple of 2 or 3 int
- buildH_calcOP_test.make_dic_lipids_with_indexes(universe_woH, dic_OP, dic_lipid)[source]
Expands dic_lipid and adds the index of each atom and helper.
- IMPORTANT: the index of each atom/helper is given with respect to the
first atom in that residue.
For example, if we have a POPC where C1 is the first atom, and C50 the last one, we want in the end: {‘C1’: (‘CH3’, ‘N4’, ‘C5’, 0, 3, 4), …,
‘C50’: (‘CH3’, ‘C49’, ‘C48’, 49, 48, 47)}
Where the 3 last int are the index (ix) of the atom, helper1, helper2 (possibly helper3) with respect to the first atom. Thus for C1 : 0 is index of C1, N4 is 3 atoms away from C1 and C5 is 4 atoms away from C1. For C50: C50 is 49 atoms away from C1, C49 is 48 atoms away from C1, C48 is 47 atoms away from C1.
- Parameters:
universe_woH (MDAnalysis Universe instance) – The universe without hydrogens.
dic_lipid (dictionnary) – Comes from dic_lipids.py. Contains carbon names and helper names needed for reconstructing hydrogens.
- Returns:
The returned dictionnary as described above in this docstring.
- Return type:
dictionnary
- buildH_calcOP_test.fast_build_all_Hs_calc_OP(universe_woH, dic_OP, dic_lipid, dic_Cname2Hnames)[source]
Build Hs and calc OP using fast indexing.
This function uses fast indexing to carbon atoms and helper atoms. It should be used when the user doesn’t want any output traj with hydrogens.
- Parameters:
universe_woH (MDAnalysis universe instance) – This is the universe without hydrogen.
dic_OP (ordered dictionnary) – Each key of this dict is a couple carbon/H, and at the beginning it contains an empty list, e.g. OrderedDict([ (‘C1’, ‘H11): [], (‘C1’, ‘H12’): [], … ]) See function init_dic_OP() below to see how it is organized.
dic_lipid (dictionnary) – Comes from dic_lipids.py. Contains carbon names and helper names needed for reconstructing hydrogens.
dic_Cname2Hnames (dictionnary) –
This dict gives the correspondance Cname -> Hname. It is a dict of tuples. If there is more than 1 H for a given C, they need to be ordered like in the PDB. e.g. for CHARMM POPC : {‘C13’: (‘H13A’, ‘H13B’, ‘H13C’), …, ‘C33’: (‘H3X’, ‘H3Y’),
…, ‘C216’: (‘H16R’, ‘H16S’), …}
- Returns:
This function returns nothing, dic_OP is changed in place.
- Return type:
None
- buildH_calcOP_test.make_dic_atname2genericname(filename)[source]
Make a dict of correspondance between generic H names and PDB names.
This dict will look like the following: {(‘C1’, ‘H11’): ‘gamma1_1’, …}. Useful for outputing OP with generic names (such as beta1, beta 2, etc.). Such files can be found on the NMRlipids MATCH repository: https://github.com/NMRLipids/MATCH/tree/master/scripts/orderParm_defs.
- Parameters:
filename (str) – Filename containing OP definition (e.g. order_parameter_definitions_MODEL_Berger_POPC.def).
- Returns:
Keys are tuples of (C, H) name, values generic name (as described above in this docstring). The use of an ordered dictionnary ensures we get always the same order in the output OP.
- Return type:
Ordered dictionnary