Skip to content

Commit ed99dba

Browse files
alsinmralsinmr
committed
Updates
1 parent 047e2df commit ed99dba

File tree

4 files changed

+281
-7
lines changed

4 files changed

+281
-7
lines changed

‎Frames/frames.py

Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -401,6 +401,121 @@ def sub():
401401

402402
return sub,frame_index,{'PPfun':'AvgGauss','sigma':sigma}
403403

404+
def lipid_chain_chi(molecule,n_bonds=0,resids=None,segids=None,filter_str=None,sigma=0):
405+
"""
406+
Returns a frame that accounts for motion around bonds in the lipid chain. This
407+
will step n_bonds up the chain to take a carbon-carbon bond as the reference
408+
frame.
409+
410+
Note that the initial selection (sel1) should only include carbons in the
411+
chains (names that start with C2 or C3)
412+
413+
414+
Parameters
415+
----------
416+
molecule : TYPE
417+
DESCRIPTION.
418+
n_bonds : TYPE, optional
419+
DESCRIPTION. The default is 1.
420+
resids : TYPE, optional
421+
DESCRIPTION. The default is None.
422+
segids : TYPE, optional
423+
DESCRIPTION. The default is None.
424+
filter_str : TYPE, optional
425+
DESCRIPTION. The default is None.
426+
sigma : TYPE, optional
427+
DESCRIPTION. The default is 0.
428+
429+
Returns
430+
-------
431+
None.
432+
433+
"""
434+
435+
sel0=selt.sel0_filter(molecule,resids=resids,segids=segids,filter_str=filter_str)
436+
437+
mdmode=molecule._mdmode
438+
molecule._mdmode=True
439+
440+
chain_list={**{f'C2{k}':f'C2{k-1}' for k in range(2,50)},
441+
**{f'C3{k}':f'C3{k-1}' for k in range(2,50)}}
442+
443+
replace={'C31':'O31','O31':'C3','C21':'O21','O21':'C1'}
444+
for name,value in replace.items():
445+
chain_list[name]=value
446+
447+
sel1,q=np.unique(molecule.sel1,return_inverse=True)
448+
for _ in range(n_bonds):
449+
s1=[]
450+
for k,s in enumerate(sel1):
451+
if s is not None:
452+
if s.name not in chain_list:
453+
s1.append(None)
454+
continue
455+
456+
i=np.logical_and(s.resid==sel0.resids,sel0.names==chain_list[s.name])
457+
if np.any(i):
458+
s1.append(sel0[i][0])
459+
else:
460+
s1.append(None)
461+
else:
462+
s1.append(None)
463+
464+
sel1=s1
465+
466+
s0=[sel1]
467+
for _ in range(2):
468+
s0.append([])
469+
for k,s in enumerate(s0[-2]):
470+
if s is not None:
471+
if s.name not in chain_list:
472+
s0[-1].append(None)
473+
continue
474+
475+
i=np.logical_and(s.resid==sel0.resids,sel0.names==chain_list[s.name])
476+
if np.any(i):
477+
s0[-1].append(sel0[i][0])
478+
else:
479+
s0[-1].append(None)
480+
else:
481+
s0[-1].append(None)
482+
483+
sel2,sel3=s0[1],s0[2]
484+
485+
486+
fi0=np.array(sel3,dtype=bool)
487+
fi=np.cumsum(fi0).astype(object)-1
488+
fi[fi0==False]=np.nan
489+
fi=fi.astype(float)
490+
491+
frame_index=fi[q]
492+
493+
494+
sel3=np.array(sel3,dtype=object)
495+
sel1=np.array(sel1,dtype=object)[sel3.astype(bool)]
496+
sel2=np.array(sel2,dtype=object)[sel3.astype(bool)]
497+
sel3=sel3[sel3.astype(bool)]
498+
499+
sel1=sel1.sum()
500+
sel2=sel2.sum()
501+
sel3=sel3.sum()
502+
503+
def sub():
504+
box=molecule.box
505+
v1=sel2.positions-sel1.positions
506+
v2=sel1.positions-sel3.positions
507+
v1=vft.pbc_corr(v1.T,box)
508+
v2=vft.pbc_corr(v2.T,box)
509+
return v1,v2
510+
511+
molecule._mdmode=mdmode
512+
513+
return sub,frame_index,{'PPfun':'AvgGauss','sigma':sigma}
514+
515+
516+
517+
518+
404519
def librations(molecule,sel1=None,sel2=None,Nuc=None,resids=None,segids=None,filter_str=None,full=True,sigma=0):
405520
"""
406521
Defines a frame for which librations are visible. That is, for a given bond,

‎Frames/special_frames.py

Lines changed: 113 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -82,15 +82,15 @@ def sub()
8282
from pyDR.Selection import select_tools as selt
8383
from copy import copy
8484

85-
def hop_setup(uni,sel1,sel2,sel3,sel4,ntest=1000):
85+
def hop_setup(mol,sel1,sel2,sel3,sel4,ntest=1000):
8686
"""
8787
Function that determines where the energy minima for a set of bonds can be
8888
found. Use for chi_hop and hop_3site.
8989
"""
9090
v12,v23,v34=list(),list(),list()
91-
box=uni.dimensions
92-
traj=uni.trajectory
93-
step=np.floor(traj.n_frames/ntest).astype(int)
91+
box=mol.box
92+
traj=mol.traj
93+
step=np.floor(len(traj)/ntest).astype(int)
9494

9595
for _ in traj[::step]:
9696
v12.append(vft.pbc_corr((sel1.positions-sel2.positions).T,box[:3]))
@@ -164,7 +164,7 @@ def chi_hop(molecule,n_bonds=1,Nuc=None,resids=None,segids=None,filter_str=None,
164164
"Next, we sample the trajectory to get an estimate of the energy minima of the hopping"
165165
#Note that we are assuming that minima are always separated by 120 degrees
166166

167-
vr=hop_setup(molecule.uni,sel1,sel2,sel3,sel4,ntest)
167+
vr=hop_setup(molecule,sel1,sel2,sel3,sel4,ntest)
168168

169169
box=molecule.box
170170
if sigma!=0:
@@ -186,6 +186,113 @@ def sub():
186186
return vft.R(v12s.T,*sc),v23s #Rotate back into original frame
187187

188188
return sub,frame_index
189+
190+
def lipid_hop(molecule,n_bonds=0,Nuc=None,resids=None,segids=None,filter_str=None,ntest=1000,sigma=0):
191+
"""
192+
Determines contributions to motion due to 120 degree hops across three sites
193+
for some bond within a side chain. Motion of the frame will be the three site
194+
hoping plus any outer motion (could be removed with additional frames), and
195+
motion within the frame will be all rotation around the bond excluding
196+
hopping.
197+
198+
One provides the same arguments as lipid_chain_chi, where we specify the
199+
nucleus of interest (ch3,ivl,ivla,ivlr,ivll, etc.), plus any other desired
200+
filters. We also provide n_bonds, which will determine how many bonds away
201+
from the methyl group (only methyl currently implemented) we want to observe
202+
the motion (usually 1 or 2).
203+
"""
204+
205+
"First we get the selections, and simultaneously determine the frame_index"
206+
207+
sel0=selt.sel0_filter(molecule,resids=resids,segids=segids,filter_str=filter_str)
208+
209+
mdmode=molecule._mdmode
210+
molecule._mdmode=True
211+
212+
chain_list={**{f'C2{k}':f'C2{k-1}' for k in range(2,50)},
213+
**{f'C3{k}':f'C3{k-1}' for k in range(2,50)}}
214+
replace={'C31':'O31','O31':'C3','C21':'O21','O21':'C1'}
215+
for name,value in replace.items():
216+
chain_list[name]=value
217+
218+
sel1,sel2=molecule.sel1,molecule.sel2
219+
220+
for _ in range(n_bonds):
221+
s2=sel1
222+
s1=[]
223+
for k,s in enumerate(sel1):
224+
if s is not None:
225+
if s.name not in chain_list:
226+
s1.append(None)
227+
continue
228+
229+
i=np.logical_and(s.resid==sel0.resids,sel0.names==chain_list[s.name])
230+
if np.any(i):
231+
s1.append(sel0[i][0])
232+
else:
233+
s1.append(None)
234+
else:
235+
s1.append(None)
236+
237+
sel1=s1
238+
sel2=s2
239+
240+
sel2,p,q=np.unique(sel2,return_index=True,return_inverse=True)
241+
sel1=np.array(sel1)[p]
242+
243+
s0=[sel2,sel1]
244+
for _ in range(2):
245+
s0.append([])
246+
for k,s in enumerate(s0[-2]):
247+
if s is not None:
248+
if s.name not in chain_list:
249+
s0[-1].append(None)
250+
continue
251+
252+
i=np.logical_and(s.resid==sel0.resids,sel0.names==chain_list[s.name])
253+
if np.any(i):
254+
s0[-1].append(sel0[i][0])
255+
else:
256+
s0[-1].append(None)
257+
else:
258+
s0[-1].append(None)
259+
260+
261+
262+
263+
fi0=np.array(s0[-1],dtype=bool)
264+
fi=np.cumsum(fi0).astype(object)-1
265+
fi[fi0==False]=np.nan
266+
fi=fi.astype(float)
267+
268+
frame_index=fi[q]
269+
270+
sel4=np.array(s0[-1])
271+
i=sel4.astype(bool)
272+
sel1,sel2,sel3,sel4=[np.array(s)[i].sum() for s in s0]
273+
274+
vr=hop_setup(molecule,sel1,sel2,sel3,sel4,ntest)
275+
276+
box=molecule.box
277+
if sigma!=0:
278+
def sub():
279+
return [vft.pbc_corr((s1.positions-s2.positions).T,box[:3]) \
280+
for s1,s2 in zip([sel1,sel2,sel3],[sel2,sel3,sel4])]
281+
return sub,frame_index,{'PPfun':'AvgHop','vr':vr,'sigma':sigma}
282+
else:
283+
284+
def sub():
285+
v12s,v23s,v34s=[vft.pbc_corr((s1.positions-s2.positions).T,box[:3]) \
286+
for s1,s2 in zip([sel1,sel2,sel3],[sel2,sel3,sel4])]
287+
v12s=vft.norm(v12s)
288+
sc=vft.getFrame(v23s,v34s)
289+
v12s=vft.R(v12s,*vft.pass2act(*sc)) #Into frame defined by v23,v34
290+
# print((v12s*vr).sum(axis=1)[:,0])
291+
i=np.argmax((v12s*vr).sum(axis=1),axis=0) #Index of best fit to reference vectors (product is cosine, which has max at nearest value)
292+
v12s=vr[i,:,np.arange(v12s.shape[1])] #Replace v12 with one of the three reference vectors
293+
return vft.R(v12s.T,*sc),v23s #Rotate back into original frame
294+
295+
return sub,frame_index
189296

190297
def hops_3site(molecule,sel1=None,sel2=None,sel3=None,sel4=None,\
191298
Nuc=None,resids=None,segids=None,filter_str=None,ntest=1000,sigma=0):
@@ -239,7 +346,7 @@ def hops_3site(molecule,sel1=None,sel2=None,sel3=None,sel4=None,\
239346
if not(sel4):
240347
sel4=selt.find_bonded(sel3,sel0,exclude=sel2,n=1,sort='cchain',d=1.65)[0]
241348

242-
vr=hop_setup(molecule.uni,sel1,sel2,sel3,sel4,ntest)
349+
vr=hop_setup(molecule,sel1,sel2,sel3,sel4,ntest)
243350

244351
box=molecule.box
245352
if sigma!=0:

‎Project/Project.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1812,7 +1812,7 @@ def modes2bonds(self,inclOverall:bool=False,calcCC='auto'):
18121812
for d in self:
18131813
if hasattr(d,'iRED') and 'Lambda' in d.iRED:
18141814
count+=1
1815-
d.modes2bonds(inclOverall=inclOverall)
1815+
d.modes2bonds(inclOverall=inclOverall,calcCC=calcCC)
18161816

18171817
# out=self[:0]
18181818
# out._index=np.setdiff1d(self.parent._index,index0)

‎Selection/select_tools.py

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@
1313
import MDAnalysis as mda
1414
import numpy as np
1515
import numbers
16+
import re
1617

1718

1819
def sel0_filter(mol,resids=None,segids=None,filter_str=None):
@@ -193,6 +194,7 @@ def sel_lists(mol,sel=None,resids=None,segids=None,filter_str=None):
193194
return sel
194195

195196
#%% Specific selections for proteins
197+
196198
def protein_defaults(Nuc:str,mol,resids:list=None,segids:list=None,filter_str:str=None)->tuple:
197199
"""
198200
Selects pre-defined pairs of atoms in a protein, where we use defaults based
@@ -323,12 +325,62 @@ def protein_defaults(Nuc:str,mol,resids:list=None,segids:list=None,filter_str:st
323325
sel2=sel0.select_atoms('name H{0}* and resname TYR PHE'.format(Nuc[-1].upper()))
324326
if Nuc[:4].lower()!='fy_z' and '1' in Nuc.lower():
325327
sel1,sel2=sel1[::2],sel2[::2]
328+
elif Nuc.lower() in ['sn','sn1','sn2','lipid']:
329+
return lipid_defaults(mol=mol,Nuc=Nuc,resids=resids,segids=segids,filter_str=filter_str)
326330
else:
327331
print('Unrecognized Nuc option')
328332
return
329333

330334
return sel1,sel2
331335

336+
def lipid_defaults(mol,Nuc:str,resids:list=None,segids:list=None,filter_str:str=None)->tuple:
337+
"""
338+
Returns some defaults for lipids. Nuc options are
339+
340+
'sn': 13C of both side chains
341+
'sn1': 13C of sn1 chain
342+
'sn2': 13C of sn2 chain
343+
'lipid': All 13C
344+
345+
Parameters
346+
----------
347+
mol : TYPE
348+
DESCRIPTION.
349+
mol : MolSelect object or AtomGroup
350+
resids : list/array/single element, optional
351+
Restrict selected residues. The default is None.
352+
segids : list/array/single element, optional
353+
Restrict selected segments. The default is None.
354+
filter_str : str, optional
355+
Restricts selection to atoms selected by the provided string. String
356+
is applied to the MDAnalysis select_atoms function. The default is None.
357+
358+
Returns
359+
-------
360+
tuple
361+
sel1: atom group
362+
sel2: atom group
363+
364+
"""
365+
sel0=sel0_filter(mol,resids,segids,filter_str)
366+
sel00=sel0_filter(mol,resids,segids)
367+
368+
sel1=sel0[sel0.types=='C']
369+
if Nuc.lower() in ['sn','sn1','sn2','sn_1','sn1_1','sn2_1']:
370+
match={'sn':'C[2/3].','sn1':'C2.','sn2':'C3.','sn_1':'C[2/3].','sn1_1':'C2.','sn2_1':'C3.'}
371+
r=re.compile(match[Nuc.lower()])
372+
i=np.array([bool(r.match(name)) for name in sel1.names])
373+
sel1=sel1[i]
374+
375+
bond=find_bonded(sel1,sel0=sel00)
376+
i=np.array([np.logical_or(b.types=='H',b.types=='O') for b in bond])
377+
sel2=np.sum(np.array(bond).T[i.T])
378+
sel2=sel2[np.logical_not(np.logical_or(sel2.names=='O21',sel2.names=='O31'))]
379+
sel1=find_bonded(sel2,sel0=sel0,n=1)[0]
380+
381+
return sel1,sel2
382+
383+
332384
def find_methyl(mol,resids=None,segids=None,filter_str=None,select=None):
333385
"""
334386
Finds methyl groups in a protein for a list of residues. Standard selection

0 commit comments

Comments
 (0)