Skip to content

Commit

Permalink
updates
Browse files Browse the repository at this point in the history
  • Loading branch information
@alsinmr
alsinmr committed Jan 25, 2024
1 parent 3a658fa commit d13f355
Show file tree
Hide file tree
Showing 4 changed files with 145 additions and 38 deletions.
Binary file added .temp.xtc_offsets.npz
View file
Binary file not shown.
166 changes: 136 additions & 30 deletions PCA/PCAmovies.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,8 +32,8 @@ def __init__(self,pca):
"""
self.pca=pca
self._n=None
self._data=None
self._direct_data=None
self._PCARef=None
self._BondRef=None
self._zrange=None
self._timescale=None
self._xtc=None
Expand All @@ -60,52 +60,44 @@ def n(self,n):
self._n=n

@property
def data(self):
def PCARef(self):
"""
Processed PCA detector analysis. Can be slow at the first call
Processed PCA detector analysis.
Returns
-------
data
data object containing a detector analysis of the PCA.
"""
if self._data is None:
d=self.pca.Data.PCARef
d.detect.r_no_opt(self.n)
noo=d.fit()
noo.detect.r_auto(self.n)
self._data=noo.fit()
if self._PCARef is None:
return self.pca.Data.PCARef

return self._data
return self._PCARef


@data.setter
@PCARef.setter
def data(self,data):
self._data=data
self._PCARef=data
self._zrange=None

@property
def direct_data(self):
def BondRef(self):
"""
Data object containing direct analysis of the bond motion (no PCA).
Data object containing analysis of the bond motion (no PCA).
Returns
-------
None.
"""
if self._direct_data is None:
d=self.pca.direct2data()
d.detect.r_no_opt(self.n)
noo=d.fit()
noo.detect.r_auto(self.n)
self._direct_data=noo.fit()
return self._direct_data
if self._BondRef is None:
return self.pca.Data.BondRef
return self._BondRef

@direct_data.setter
def direct_data(self,data):
self._direct_data=data
@BondRef.setter
def BondRef(self,data):
self._BondRef=data

@property
def sens(self):
Expand All @@ -117,7 +109,7 @@ def sens(self):
sens
"""
return self.data.sens
return self.PCARef.sens

@property
def select(self):
Expand Down Expand Up @@ -259,6 +251,101 @@ def weight(self,timescale):
return self.data.R@wt0


def xtc_bond(self,index:int,rho_index:int,frac:float=0.75,filename:str='temp.xtc',
nframes:int=150,framerate:int=15,scaling:float=1):
"""
Parameters
----------
index : int
DESCRIPTION.
rho_index : int, optional
DESCRIPTION. The default is None.
frac : float, optional
DESCRIPTION. The default is 0.75.
filename : str, optional
DESCRIPTION. The default is 'temp.xtc'.
nframes : int, optional
DESCRIPTION. The default is 300.
framerate : int, optional
DESCRIPTION. The default is 15.
Returns
-------
None.
"""
timescale=self.PCARef.sens.info['z0'][rho_index]
step=np.round(10**timescale*1e12/self.PCARef.source.select.traj.dt/framerate).astype(int)
if step==0:step=1
if step*nframes>len(self.pca.traj):
step=np.round(len(self.pca.traj)/nframes).astype(int)

wt=self.pca.Weighting.bond(index=index,rho_index=rho_index,PCAfit=self.PCARef,frac=frac)

Delta=np.concatenate([np.zeros([wt.sum(),1]),
np.diff(self.pca.PCamp[wt,:step*nframes:step],axis=1)],axis=1)
pos=self.pca.pos[0]+scaling*np.cumsum(self.pca.PCxyz[:,:,wt]@Delta,axis=-1).T


ag=copy(self.pca.atoms)
with Writer(filename, n_atoms=len(ag)) as w:
for pos0 in pos:
# ag.positions=self.pca.mean+((wt[:,k]*self.pca.PCamp[:,i[k]])*self.pca.PCxyz).sum(-1).T
ag.positions=pos0
w.write(ag)

self._xtc=filename

return self

def xtc_mode(self,mode_index,filename:str='temp.xtc',nframes:int=150,
framerate:int=15,scaling:float=1):
"""
Displays a specific mode
Parameters
----------
mode_index : TYPE
DESCRIPTION.
filename : str, optional
DESCRIPTION. The default is 'temp.xtc'.
nframes : int, optional
DESCRIPTION. The default is 150.
framerate : int, optional
DESCRIPTION. The default is 15.
Returns
-------
None.
"""

rho_index=np.argmax(self.PCARef.R[mode_index])
timescale=self.sens.info['z0'][rho_index]
step=np.round(10**timescale*1e12/self.PCARef.source.select.traj.dt/framerate).astype(int)
if step==0:step=1
if step*nframes>len(self.pca.traj):
step=np.round(len(self.pca.traj)/nframes).astype(int)

Delta=np.concatenate([np.zeros(1),
np.diff(self.pca.PCamp[mode_index,:step*nframes:step],axis=0)],axis=0)
pos=self.pca.pos[0]+scaling*np.cumsum(np.atleast_3d(self.pca.PCxyz[:,:,mode_index])@\
np.atleast_2d(Delta),axis=-1).T
ag=copy(self.pca.atoms)
with Writer(filename, n_atoms=len(ag)) as w:
for pos0 in pos:
# ag.positions=self.pca.mean+((wt[:,k]*self.pca.PCamp[:,i[k]])*self.pca.PCxyz).sum(-1).T
ag.positions=pos0
w.write(ag)

self._xtc=filename

return self



def xtc_log_swp(self,filename:str='temp.xtc',zrange=None,nframes:int=300,framerate:int=15):
"""
Create an xtc file for a log sweep of the timescales. Setup is the output
Expand Down Expand Up @@ -302,7 +389,10 @@ def xtc_log_swp(self,filename:str='temp.xtc',zrange=None,nframes:int=300,framera

PCamp=np.cumsum(DelPC,axis=-1)+np.atleast_2d(PCcorr).T@np.atleast_2d(i)/i[-1]

pos=self.pca.pos[0].T+np.array([(self.pca.PCxyz*a).sum(-1) for a in PCamp.T])
pos=self.pca.pos[0]+(self.pca.PCxyz@PCamp).T
# pos=self.pca.pos[0].T+np.array([(self.pca.PCxyz*a).sum(-1) for a in PCamp.T])



# pos=self.pca.pos[0]+(np.cumsum(DelPC,axis=-1)*self.pca.PCxyz).sum(-1).T+\
# +(np.atleast_2d(PCcorr).T@np.atleast_2d(i)/i[-1])
Expand All @@ -315,7 +405,7 @@ def xtc_log_swp(self,filename:str='temp.xtc',zrange=None,nframes:int=300,framera
with Writer(filename, n_atoms=len(ag)) as w:
for pos0 in pos:
# ag.positions=self.pca.mean+((wt[:,k]*self.pca.PCamp[:,i[k]])*self.pca.PCxyz).sum(-1).T
ag.positions=pos0.T
ag.positions=pos0
w.write(ag)

self._xtc=filename
Expand All @@ -340,6 +430,8 @@ def play_xtc(self):

return self



class Options():
def __init__(self,pca_movie):
self.dict={}
Expand Down Expand Up @@ -392,6 +484,20 @@ def TimescaleIndicator(self,tau=None,remove:bool=False):
self.dict['TimescaleIndicator']=lambda tau=tau:self.add_event('TimescaleIndicator', tau)
return self

def Detectors(self,index=None,rho_index=None,remove:bool=False):
"""
Returns
-------
None.
"""

out=dict(R=R,rho_index=rho_index,ids=ids)
# CMXRemote.remove_event(ID,'Detectors')
CMXRemote.add_event(ID,'Detectors',out)

def DetFader(self,tau=None,index=None,rho_index=None,remove:bool=False):
"""
Add a Detector Fader to the trajectory
Expand Down Expand Up @@ -419,7 +525,7 @@ def DetFader(self,tau=None,index=None,rho_index=None,remove:bool=False):
if __name__ in self.dict:self.dict.pop(__name__)
return

data=self.pca_movie.direct_data
data=self.pca_movie.BondRef
rhoz=data.sens.rhoz

if index is None:index=np.ones(data.R.shape[0],dtype=bool)
Expand All @@ -429,7 +535,7 @@ def DetFader(self,tau=None,index=None,rho_index=None,remove:bool=False):
if rho_index is None:
rho_index=np.ones(data.R.shape[1],dtype=bool)
if rhoz[-1][-1]>0.9:rho_index[-1]=False
print(rho_index)


x=x[:,rho_index]
rhoz=rhoz[rho_index]
Expand Down
14 changes: 7 additions & 7 deletions PCA/PCAsubs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -651,7 +651,7 @@ def __init__(self,pca):
self._n=None
self._n_noopt=12
self._locs={}
self._DirectRef=None
self._BondRef=None
self._PCARef=None

def clear(self):
Expand All @@ -669,7 +669,7 @@ def clear(self):
self._locs={}

@property
def DirectRef(self):
def BondRef(self):
"""
Default data set when a directly-calculated data object is required.
Expand All @@ -681,14 +681,14 @@ def DirectRef(self):
data
"""
if self._DirectRef is not None:
return self._DirectRef
if self._BondRef is not None:
return self._BondRef
return self.direct()

@DirectRef.setter
def DirectRef(self,data):
@BondRef.setter
def BondRef(self,data):
assert data.source.status in ['proc','opt_fit'],'Data needs to have status proc or opt_fit'
self._DirectRef=data
self._BondRef=data

@property
def PCARef(self):
Expand Down
3 changes: 2 additions & 1 deletion PCA/testPCA.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,9 +11,10 @@
from pyDR.PCA.PCAclean import PCA

select=pyDR.MolSelect(topo='/Users/albertsmith/Documents/Dynamics/MDsims/HETs/backboneB.pdb',
traj_files='/Users/albertsmith/Documents/Dynamics/MDsims/HETs/backboneB.xtc',step=100)
traj_files='/Users/albertsmith/Documents/Dynamics/MDsims/HETs/backboneB.xtc',step=1,tf=10000)

pca=PCA(select)
pca.select_bond('15N',resids=np.concatenate([np.arange(225,248),np.arange(261,281)]))
pca.select_bond('15N')
pca.select_atoms('name N C CA')

0 comments on commit d13f355

Please sign in to comment.