2
2
from PyCEC .cec_system import CECSystem
3
3
4
4
# MDAnalysis
5
- import MDAnalysis as mda # TODO: This needs to be fixed tbh, not efficient
5
+ import MDAnalysis as mda # TODO: This needs to be fixed tbh, not efficient
6
6
from MDAnalysis .analysis .hydrogenbonds .hbond_analysis import HydrogenBondAnalysis as HBA
7
7
8
8
from tqdm import tqdm
9
- import time
10
9
11
10
from matplotlib import pyplot as plt
12
11
import numpy as np
@@ -16,9 +15,9 @@ class CVAnalysis:
16
15
"""
17
16
Class to analyse the collective variable of the CEC system.
18
17
"""
19
-
20
- def __init__ ( self , universe , initial_resid , target_resid , other_resids = None ,
21
- ligand = None , cyzone_dim = [ 8 , 8 , - 8 ], frame_n = 0 ):
18
+ def __init__ ( self , universe , initial_resid , target_resid ,
19
+ other_resids = None , ligand = None , cyzone_dim = [ 8 , 8 , - 8 ] ,
20
+ frame_n = 0 ):
22
21
# Attributes
23
22
self .u = universe
24
23
self .initial_resid = initial_resid
@@ -40,31 +39,29 @@ def __init__(self, universe, initial_resid, target_resid, other_resids=None,
40
39
def get_traj_times (self ):
41
40
"""
42
41
Get the time of each frame in the trajectory.
43
-
42
+
44
43
"""
45
44
times = []
46
45
for ts in self .u .trajectory :
47
46
times .append (ts .time )
48
47
49
48
return times
50
49
51
-
52
50
def get_water_counts (self ):
53
51
"""
54
52
Iterate through the frames and get the water counts.
55
-
53
+
56
54
"""
57
55
water_counts = []
58
- for ts in tqdm (self .u .trajectory , desc = "Getting water counts..." ): # tqdm progress bar
56
+ for ts in tqdm (self .u .trajectory , desc = "Getting water counts..." ):
59
57
self .cv .set_frame (ts .frame )
60
58
water_counts .append (len (self .cv .waters .atoms ))
61
59
62
60
return water_counts
63
61
64
-
65
62
def get_water_counts_h (self ):
66
63
"""
67
- Iterate through the frames and get the water counts in a readable format.
64
+ Iterate through the frames and get water counts in a readable format.
68
65
69
66
"""
70
67
water_counts = self .get_water_counts ()
@@ -78,7 +75,7 @@ def get_water_counts_h(self):
78
75
for i , t , w in water_c :
79
76
print (f"Frame { i } at time { t } ps has { w } water molecules." )
80
77
81
- #return water_c
78
+ # return water_c
82
79
83
80
84
81
class WaterAnalysis :
@@ -94,7 +91,7 @@ def __init__(self, universe, initial_resid, target_resid, other_resids,
94
91
self .other_resids = other_resids
95
92
self .ligand = ligand
96
93
self .cyzone_dim = cyzone_dim
97
- self .frame_n = frame_n # TODO: is frame relevant here?
94
+ self .frame_n = frame_n # TODO: is frame relevant here?
98
95
99
96
# Trajectory information
100
97
self .n_frames = len (self .u .trajectory )
@@ -109,39 +106,37 @@ def __init__(self, universe, initial_resid, target_resid, other_resids,
109
106
self .waters = self .cv .waters
110
107
self .waters_sele_str = self .cv .waters_sele_str
111
108
112
-
113
109
def get_traj_times (self ):
114
110
"""
115
111
Get the time of each frame in the trajectory.
116
-
112
+
117
113
"""
118
114
times = []
119
115
for ts in self .u .trajectory :
120
116
times .append (ts .time )
121
117
122
118
return times
123
119
124
- def get_hydrogen_bonds (self , donors_sel = None , hydrogens_sel = None , acceptors_sel = None ,
125
- distance_cutoff = 3.5 , angle_cutoff = 120.0 ):
120
+ def get_hydrogen_bonds (self , donors_sel = None , hydrogens_sel = None ,
121
+ acceptors_sel = None , distance_cutoff = 3.5 ,
122
+ angle_cutoff = 120.0 ):
126
123
"""
127
124
Get the hydrogen bonds between the a pair of molecules.
128
-
125
+
129
126
"""
130
127
# Create the HBA object
131
128
hba = HBA (universe = self .u , donors_sel = donors_sel ,
132
129
hydrogens_sel = hydrogens_sel , acceptors_sel = acceptors_sel ,
133
130
d_h_cutoff = distance_cutoff , d_h_a_angle_cutoff = angle_cutoff )
134
- hba .run () # run the analysis
131
+ hba .run () # run the analysis
135
132
136
133
# Get counts
137
134
counts = hba .count_by_time ()
138
-
139
- return counts
140
-
135
+
141
136
def get_water_h_bonds (self , distance_cutoff = 3.5 , angle_cutoff = 120.0 ):
142
137
"""
143
138
Get the hydrogen bonds between the water molecules.
144
-
139
+
145
140
"""
146
141
# Get the hydrogen bonds
147
142
counts = self .get_hydrogen_bonds (donors_sel = self .waters_sele_str ,
@@ -150,30 +145,29 @@ def get_water_h_bonds(self, distance_cutoff=3.5, angle_cutoff=120.0):
150
145
distance_cutoff = distance_cutoff ,
151
146
angle_cutoff = angle_cutoff )
152
147
return counts
153
-
148
+
154
149
def plot_hbond_counts (self , counts , save = False ):
155
150
"""
156
151
Plot the hydrogen bond counts.
157
-
152
+
158
153
"""
159
154
fig , ax = plt .subplots ()
160
155
ax .plot (self .times , counts )
161
156
ax .set_ylabel ("Number of hydrogen bonds" )
162
157
ax .set_xlabel ("Frame number" )
163
158
ax .set_title ("Hydrogen bond counts over time" )
164
-
159
+
165
160
if save :
166
161
fig .savefig ("hbond_counts.png" , dpi = 300 )
167
-
168
- plt .show ()
169
162
163
+ plt .show ()
170
164
171
165
def get_max_hbond_counts (self , counts ):
172
166
"""
173
167
Get the frames and times of the 10 frames with the most hydrogen bonds.
174
168
"""
175
169
# Get the frames with the most hydrogen bonds
176
- max_counts = sorted (counts , reverse = True )[:10 ] # get the 10 highest counts
170
+ max_counts = sorted (counts , reverse = True )[:10 ] # get 10 highest counts
177
171
print (f"\n Max counts: { max_counts } " )
178
172
max_frames = []
179
173
for c in max_counts :
@@ -184,7 +178,7 @@ def get_max_hbond_counts(self, counts):
184
178
print (f"Times of the frames with the most hydrogen bonds: { max_times } " )
185
179
186
180
return max_frames , max_times
187
-
181
+
188
182
def find_consistent_waters (self ):
189
183
"""
190
184
Find waters consistently included in the waters list.
@@ -198,8 +192,6 @@ def get_residence_times(self):
198
192
pass
199
193
200
194
201
-
202
-
203
195
if __name__ == "__main__" :
204
196
205
197
# Directory and title
@@ -225,4 +217,4 @@ def get_residence_times(self):
225
217
226
218
counts = wa .get_water_h_bonds ()
227
219
228
- wa .plot_hbond_counts (counts , save = True )
220
+ wa .plot_hbond_counts (counts , save = True )
0 commit comments