-
Notifications
You must be signed in to change notification settings - Fork 1
/
mdlib.py
157 lines (128 loc) · 4.99 KB
/
mdlib.py
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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
# Author: Samuel Genheden samuel.genheden@gmail.com
"""
Classes to help with the processing of MD trajectories
"""
import argparse
import MDAnalysis as md
import numpy as np
import elbalib
from sgenlib import moldyn
from sgenlib import mdactions
class _Residue(object):
"""
Class to store residue information, used by the AaCgTrajectoryProcessor class
Attributes
----------
ffmol : elbalib.ElbaMolecule
the force field molecule corresponding to this residue
mdsel : MDAnalysis.AtomSelection
the MD universe selection corresponding to this residue
resid : string
the identifier of this string
transmat : numpy.ndarray
transformation from AA to CG coordinates for this residue
xyz : numpy.ndarray
the current coordinates of this residue
"""
def __init__(self, id, ffmol, mduniverse):
self.resid = "resid %d" % id
self.mdsel = mduniverse.selectAtoms(self.resid)
self.ffmol = ffmol
self.transmat = None
self.xyz = None
def make_transmat(self,namtrans):
"""
Creates and store a transformation matrix from AA to CG for this residue
"""
if namtrans is None:
atomnames = [a.name.lower() for a in self.mdsel]
else:
atomnames = [namtrans[a.name].lower() for a in self.mdsel]
self.transmat = self.ffmol.transmat(atomnames)
def update_coordinates(self):
"""
Updates the xyz property by using the coordinates from mdsel,
if necessary transforms the coordinates
"""
if self.transmat is None:
self.xyz = self.mdsel.get_positions()
else:
self.xyz = np.dot(self.transmat,self.mdsel.get_positions())
class AaCgTrajectoryProcessor(moldyn.TrajectoryProcessor):
"""
A processor for a trajectory that can be either AA or CG
Attributes
----------
ff : elbalib.Elba
the force field definitions
namtrans : dictionary of strings
name translations
iscgtraj : boolean
if this is CG trajectory
refuniverse : MDAnalysis.Universe
the reference CG universe
residues : list of _Residue
residues to keep track of
"""
def __init__(self, description):
super(AaCgTrajectoryProcessor,self).__init__(description)
self.argparser.add_argument('-r', '--ref', help="a CG reference file", default="ref.pdb")
self.argparser.add_argument('-x', '--xml', help="an XML file with force field definitions")
self.argparser.add_argument('-n', '--naming', nargs=2, help="naming convention change")
self.argparser.add_argument('-m','--mol',help="the name of the molecule")
self._argnames.extend(["ref", "xml", "naming","mol"])
self.namtrans = None
self.ff = None
self.refuniverse = None
self.residues = []
self.iscgtraj = False
def setup(self,printargs=False):
super(AaCgTrajectoryProcessor,self).setup(printargs)
self.ff = elbalib.Elba()
self.ff.load(self.args.xml)
self.refuniverse = md.Universe(self.args.ref)
# Check if we need to translate the atom names
self.namtrans = None
if self.args.naming is not None:
namref = [line.strip() for line in open(self.args.naming[0], 'r').readlines()]
nammob = [line.strip() for line in open(self.args.naming[1], 'r').readlines()]
self.namtrans = {}
for nr, nm in zip(namref, nammob):
if nm != "*":
self.namtrans[nm] = nr
# Setup of a list of residues defined in the force field that can be processed
self.residues = []
for r in self.refuniverse.selectAtoms("all").residues:
if self.args.mol is not None and r.name.lower() != self.args.mol.lower():
continue
ffmol = self.ff.find_molecule(r.name)
if ffmol is None or not ffmol.beadnames:
continue
self.residues.append(_Residue(r.id, ffmol, self.universe))
if not self.iscgtraj : self.residues[-1].make_transmat(self.namtrans)
for action in self.actions:
action.setup_residues()
def setup_universe(self):
self.iscgtraj = self.args.struct is None
if self.iscgtraj:
self.args.struct = self.args.ref
super(AaCgTrajectoryProcessor,self).setup_universe()
def call_actions(self):
for action in self.actions:
action.preprocess()
for res in self.residues:
res.update_coordinates()
for action in self.actions:
action.resprocess(res)
for action in self.actions:
action.process()
if self.dosubsample and action.dosubsample and \
self.currtime % self.freq == 0:
action.subsample()
class AaCgAction(mdactions.TrajectoryAction):
def setup_residues():
pass
def preprocess(self):
pass
def resprocess(self,res):
pass