This notebook shows how to compute pairwise ditances, angles and dihedrals between CG beads given a CG mapping. The CG mapping used in this example is generated from DSGPM model. You must need MDAnalysis and NetworkX in your working environment to run this example.
import hoomd
import hoomd.htf as htf
import tensorflow as tf
import MDAnalysis as mda
import numpy as np
import os
os.environ['CUDA_VISIBLE_DEVICES'] = '-1'
import warnings
warnings.filterwarnings('ignore')
2 Physical GPUs, 2 Logical GPUs
In this example we read from a trajectory file to compute coarse-grained (CG) bond distances, angles and dihedrals. Here, we use two MDAnalysis universes with and without hydrogens as the mapping model used to compute the CG mappings (DSGPM model) only maps heavy atoms of a given molecule. Hence, we have to add the missing hydrogen atoms to the corresponding CG beads.
universe = mda.Universe('ex9_segA.pdb','ex9_segA.trr')
avg_cgr = tf.keras.metrics.MeanTensor()
avg_cga = tf.keras.metrics.MeanTensor()
avg_cgd = tf.keras.metrics.MeanTensor()
directory = os.getcwd()
jfile = os.path.join(directory,'ex9_cgmap_segA.json')
#mda universe without H's
uxh = mda.Universe(os.path.join(directory,'ex9_segA_xH.pdb'))
#mda universe with H's
uh = mda.Universe(os.path.join(directory,'ex9_segA.pdb'))
for inputs, ts in htf.iter_from_trajectory(16, universe, r_cut=10, start=400, end=700):
cg_fts = []
r_tensor = []
a_tensor = []
d_tensor = []
box = inputs[2]
#get CG bead indices that make bonds, angles, dihedrals and
#CG coordinates
cg_fts = htf.compute_cg_graph(DSGPM=True,infile=jfile,group_atoms=True,
u_no_H=uxh, u_H=uh)
for i in range(len(cg_fts[0])):
cg_r = htf.mol_bond_distance(CG = True, cg_positions = cg_fts[-1],
b1=cg_fts[0][i][0],b2=cg_fts[0][i][1],box=box)
r_tensor.append(cg_r)
for j in range(len(cg_fts[1])):
cg_a = htf.mol_angle(CG= True, cg_positions=cg_fts[-1],
b1=cg_fts[1][j][0],b2=cg_fts[1][j][1],b3=cg_fts[1][j][2],box=box)
a_tensor.append(cg_a)
for k in range(len(cg_fts[2])):
cg_d = htf.mol_dihedral(CG=True,cg_positions=cg_fts[-1],b1=cg_fts[2][k][0],
b2=cg_fts[2][k][1],b3=cg_fts[2][k][2],b4=cg_fts[2][k][3],box=box)
d_tensor.append(cg_d)
avg_cgr.update_state(r_tensor)
avg_cga.update_state(a_tensor)
avg_cgd.update_state(d_tensor)
100%|██████████| 801/801 [01:30<00:00, 8.83it/s]
cgR = avg_cgr.result().numpy()
cgD = avg_cgd.result().numpy()*180./np.pi
cgA = avg_cga.result().numpy()*180./np.pi
print('CG pairwise distances:',cgR,'\n')
print('CG angles:',cgA,'\n')
print('CG dihedral angles:',cgD)
CG pairwise distances: [ 5.4447026 1.1312671 6.8375177 2.9382892 2.4656532 4.4416947 3.199694 4.2150507 3.5845404 2.153627 7.9029765 3.8829455 6.7589035 6.4774413 2.255304 4.924929 15.143286 ] CG angles: [ 57.06865 75.22357 83.657074 113.90926 30.8918 61.174572 40.556293 27.594091 50.535973 149.74725 46.7441 91.21376 44.42922 157.15317 45.61479 121.53312 140.93109 90.67879 51.733078 156.72841 ] CG dihedral angles: [ 61.196575 177.25443 4.7860584 111.41965 176.07312 133.15497 84.99461 135.7665 147.13869 4.834345 168.7402 124.28182 175.61597 21.146255 163.78894 32.634514 9.021241 175.17809 10.565324 7.1954145]
Note that the above computation is only applied to one molecule in the system. If a user has multiple copies of a molecule, the calculation of indices of CG beads making bonds, angles and dihedrals must be applied to all molecules.
Let's assume there are 2 of the above molecules are available in the system. Each molecule has 18 CG beads. We can obtain the indices as follows. Here we use the outputs from compute_cg_graph
to both molecules.
r_ids,a_ids,d_ids = htf.mol_features_multiple(bnd_indices=cg_fts[0],ang_indices=cg_fts[1],
dih_indices=cg_fts[2],molecules=2,beads=18)
# For example here are the CG bead indices involved in making angles
print('angles in molecule 1: ',a_ids[:20])
print('\n angles in molecule 2:',a_ids[20:])
# Now the same calculation with mol_bond_distance,mol_angle and mol_dihedral can be used
# to calculate CG bond distances, angles and dihedrals
angles in molecule 1: [[ 0 1 2] [ 1 2 3] [ 2 3 5] [ 3 5 4] [ 3 5 6] [ 4 5 6] [ 5 6 7] [ 5 6 8] [ 6 8 10] [ 7 6 8] [ 8 10 9] [ 8 10 11] [ 9 10 11] [10 11 12] [11 12 13] [12 13 15] [13 15 14] [13 15 16] [14 15 16] [15 16 17]] angles in molecule 2: [[18 19 20] [19 20 21] [20 21 23] [21 23 22] [21 23 24] [22 23 24] [23 24 25] [23 24 26] [24 26 28] [25 24 26] [26 28 27] [26 28 29] [27 28 29] [28 29 30] [29 30 31] [30 31 33] [31 33 32] [31 33 34] [32 33 34] [33 34 35]]