In [1]:

    

import mdtraj as md
traj = md.load('./material/md.pdb')
top = traj.topology
print traj
print top


    

    




<mdtraj.Trajectory with 101 frames, 4018 atoms, 401 residues, and unitcells>
<mdtraj.Topology with 1 chains, 401 residues, 4018 atoms, 4088 bonds>

In [2]:

    

# traj.xyz[frame_idx, atom_idx,:]
traj.xyz[0,0]


    

    
Out[2]:





array([ 2.63599992,  1.25999999,  9.58300018], dtype=float32)

In [3]:

    

cas = top.select('name CA')
print cas
len(cas)


    

    




[   4   12   22   35   44   65   73   82   90  108  126  131  138  146  154
  175  188  198  204  213  222  231  240  257  264  270  278  287  293  306
  312  330  339  348  358  366  380  391  399  420  426  435  449  455  462
  469  477  486  494  502  512  520  532  542  550  558  567  577  588  596
  605  615  629  646  657  666  687  700  711  722  731  739  749  761  770
  782  792  801  810  819  827  836  857  866  878  886  895  907  915  922
  930  943  952  960  968  977  984  992  997 1003 1008 1016 1023 1032 1041
 1049 1057 1066 1075 1087 1093 1099 1107 1120 1129 1137 1154 1163 1171 1179
 1187 1196 1210 1228 1235 1240 1248 1254 1259 1276 1282 1291 1300 1313 1320
 1331 1340 1353 1362 1379 1390 1395 1408 1412 1420 1427 1440 1451 1459 1467
 1476 1484 1496 1503 1512 1527 1532 1541 1557 1565 1573 1581 1589 1598 1610
 1619 1628 1637 1648 1653 1661 1670 1676 1686 1696 1706 1714 1722 1731 1748
 1756 1765 1776 1793 1802 1813 1824 1830 1843 1852 1861 1870 1878 1890 1899
 1912 1922 1930 1938 1948 1957 1968 1975 1984 2000 2008 2019 2030 2041 2050
 2067 2080 2088 2097 2111 2120 2124 2132 2137 2154 2160 2177 2195 2204 2213
 2218 2228 2237 2246 2251 2260 2269 2286 2298 2304 2318 2325 2336 2345 2353
 2370 2376 2389 2410 2421 2430 2439 2448 2461 2473 2482 2490 2499 2512 2521
 2538 2548 2560 2577 2587 2598 2611 2620 2629 2637 2654 2665 2679 2687 2695
 2700 2705 2713 2721 2731 2740 2748 2757 2771 2779 2796 2807 2814 2819 2824
 2834 2851 2868 2886 2893 2904 2912 2921 2933 2942 2959 2970 2978 2987 3008
 3019 3030 3041 3050 3060 3065 3073 3084 3095 3104 3114 3119 3130 3139 3148
 3157 3165 3173 3180 3197 3206 3219 3231 3240 3249 3260 3269 3290 3302 3312
 3320 3325 3338 3344 3353 3371 3376 3383 3391 3400 3417 3422 3434 3443 3460
 3467 3475 3483 3494 3503 3512 3517 3526 3535 3544 3553 3570 3579 3584 3589
 3598 3609 3619 3630 3635 3644 3654 3663 3680 3696 3704 3709 3714 3719 3728
 3737 3754 3763 3774 3795 3812 3820 3830 3839 3857 3870 3888 3901 3909 3917
 3930 3939 3948 3956 3965 3970 3978 3983 3991 4000 4013]






    
Out[3]:





401

In [4]:

    

import numpy as np
cm = np.zeros((401,401))


    

In [5]:

    

for f in range(traj.n_frames):
    for i, cur_i in enumerate(cas):
        cur_ca = traj.xyz[f, cur_i]    
        for j, com_i in enumerate(cas):
            dis = np.sqrt(np.sum((traj.xyz[f,com_i]-cur_ca)**2))
            if dis < 6.0:
                cm[i][j] = cm[i][j] + 1


    

In [6]:

    

%matplotlib inline
import matplotlib.pyplot as plt
plt.imshow(cm, interpolation='nearest', cmap=plt.cm.ocean)
plt.colorbar()
plt.show()


    

In [9]:

    

cm


    

    
Out[9]:





array([[ 101.,   87.,   84., ...,   63.,   44.,   59.],
       [  87.,  101.,   72., ...,   75.,   50.,   63.],
       [  84.,   72.,  101., ...,   50.,   39.,   60.],
       ..., 
       [  63.,   75.,   50., ...,  101.,   76.,   73.],
       [  44.,   50.,   39., ...,   76.,  101.,   70.],
       [  59.,   63.,   60., ...,   73.,   70.,  101.]])

In [ ]: