1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
|
# Calculate atomic fluctuations from a trajectory.
#
# The trajectory file contains a protein. The fluctuations are calculated
# for all atoms, and then the average over the C-alpha atoms is determined.
# Note that calculating the fluctuations for only the C-alpha atoms is
# more complicated and no faster.
#
# This example illustrates:
# 1) Reading trajectory files
# 2) Selecting parts of a system
# 3) Calculating trajectory averages
#
from MMTK import *
from MMTK.Trajectory import Trajectory
from MMTK.Proteins import Protein
# Open the trajectory, use every tenth step
trajectory = Trajectory(None, 'lysozyme.nc')[::10]
universe = trajectory.universe
# Calculate the average conformation
average = ParticleVector(universe)
for step in trajectory:
average += step['configuration']
average /= len(trajectory)
# Calculate the fluctiations for all atoms
fluctuations = ParticleScalar(universe)
for step in trajectory:
d = step['configuration']-average
fluctuations += d*d
fluctuations /= len(trajectory)
# Collect the C-alpha atoms
calphas = Collection()
for protein in universe.objectList(Protein):
for chain in protein:
for residue in chain:
calphas.addObject(residue.peptide.C_alpha)
# Average over the C-alpha atoms
calpha_mask = calphas.booleanMask()
print (fluctuations*calpha_mask).sumOverParticles()/calphas.numberOfAtoms()
|
