Commit f8332608 authored by Stelios Karozis's avatar Stelios Karozis

Initial support for gmx commands

parent ab32f8b9
......@@ -6,6 +6,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
## [Unreleased]
## [0.0.6] - 2020-07-02
### Added
- Initial support for gmx commands
- Add gmx density profile function
### Changed
- more dynamic I/O of gmx_ndx
### Removed
- None
## [0.0.5] - 2020-05-21
### Added
- Add progress bar per function
......
import os
import tooba_f as tbf
import tooba_gmx as tbgmx
###################################################
#NOTICE: resids of head in each subdomain may differ in tail case
# keep all atoms of group in the first occurent subdomain
......@@ -12,6 +13,7 @@ DISCET=[6, 6, 6]
NUM_FR=1
TRAJ=SYSTEM_NAME+'/eq_traj.trr'
GRO=SYSTEM_NAME+'/eq_final.gro'
TPR=SYSTEM_NAME+'/eq_run.tpr'
ITP_DIC={'NS':'CER_SOVOVA.itp','FFA':'FFA_CG.itp','CHOL':'CHOL_CG.itp'}
###################################################
# {NAME:[QUEUE OF PROCESSES]}
......@@ -19,18 +21,26 @@ ITP_DIC={'NS':'CER_SOVOVA.itp','FFA':'FFA_CG.itp','CHOL':'CHOL_CG.itp'}
# NAME: It is user defined. A dictionary must follows with the same name.
# The dict structure has to be: {res_type:[atom_types]}
#
# if NAME is COMBINATION then it needs part or all tha info from aforementioned
# if NAME is COMBINE then it needs part or all the info from aforementioned
# groups to execute a process. You cannot use combination as first group.
#
# QUEUE OF PROCESSES: surf, vector, [save, [type], save_name]
# type: pkl, json, gmx_ndx (one ndx for every subdomain)
# QUEUE OF PROCESSES: surf, vector, tilt, index, density, gmx_ndx, [save, [type], save_name]
#
# surf: Determine surface from atoms (ex. Head of lipid)
# vector: Determine vector that fits atoms (ex. Tail of lipid)
# tilt: Use surf and vector result to calculate angle (if NAME is COMBINE)
# index: Creates unique code (md5) for every subdomain to use in data saving process
# density: Detrmine density profile of x,y,z and save peaks of directions with the least number
# gmx_ndx: Saves one ndx for every subdomain
# [save, [type], save_name]: Save function and properties aka, type: pkl, json, Name
#
# if NAME is COMBINATION: [tilt]
###################################################
GROUPS={'HD_GROUP':['surf',['save', ['pkl','gmx_ndx', 'json'],'time_domain_c-normal-cg']],
GROUPS={'ALL':['gmx_ndx','index',['save', ['pkl'],'index'],'density',['save', ['pkl'],'dens']],
'HD_GROUP':['surf',['save', ['pkl', 'json'],'time_domain_c-normal-cg']],
'TL_GROUP':['vector'],
'COMBINE':['tilt']
'COMBINE':[['HD_GROUP','surf'],['TL_GROUP','vector'],['COMB','tilt']]
}
ALL={'NS':['C6', 'Na', 'P4', 'P3', 'C7','C3', 'C4', 'C5', 'C8', 'C9', 'C10'], 'CHOL':['ROH','R1', 'R2', 'R3', 'R4', 'R5'], 'FFA':['AC','C1', 'C2', 'C3', 'C4']}
HD_GROUP={'NS':['C6', 'Na', 'P4', 'P3', 'C7'], 'CHOL':['ROH'], 'FFA':['AC']}
TL_GROUP={'NS':['C3', 'C4', 'C5', 'C8', 'C9', 'C10'], 'CHOL':['R1', 'R2', 'R3', 'R4', 'R5'], 'FFA':['C1', 'C2', 'C3', 'C4']}
###################################################
......@@ -74,7 +84,7 @@ for i in GROUPS.keys():
group_ndx=tbf.frompickle('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_ndx.pkl')
else:
#Find atom type index in lists created above
group_ndx=tbf.atomid_data(res_num, res_type, atom_type, atom_num, group=HD_GROUP)
group_ndx=tbf.atomid_data(res_num, res_type, atom_type, atom_num, group=locals()[i])
tbf.topickle(fl=group_ndx, sv_name='./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_ndx')
###################################################
if os.path.isfile('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_box.pkl'):
......@@ -98,25 +108,65 @@ for i in GROUPS.keys():
for j in GROUPS[i]:
if len(j) > 1:
if j=='surf':
if j not in locals():
surf={}
#Creates dictionary with c, normal per subdomain for each frame
surf=tbf.coord2norm2cg(coord_vector,img=False)
sv_data=surf
surf[locals()[i]]=tbf.coord2norm2cg(coord_vector,img=False)
sv_data=surf[locals()[i]]
if j=='vector':
vector=tbf.coord2vector(coord_vector)
sv_data=vector
if vector not in locals():
vector={}
vector[locals()[i]]=tbf.coord2vector(coord_vector)
sv_data=vector[locals()[i]]
if j=='index':
uniq_id=tbgmx.ndx_index(SYSTEM_NAME)
sv_data=uniq_id
if j=='density':
dens_df={}
for iidd in uniq_id.keys():
dens_df[iidd]={}
fl='./'+uniq_id[iidd]['system']+'/gmx_ndx/'+uniq_id[iidd]['domain']
cnt=-1
for mol in locals()[i].keys():
cnt=cnt+1
for d in ('x','y','z'):
peaks = tbgmx.density_picks(TRR=TRAJ,TPR=TPR,IND=fl,SLC=400,ST=1000,EN=-1,normal=d,fld='./'+uniq_id[iidd]['system'],arg=cnt,dist_pk=20)
if d=='x':
tmp=peaks
else:
print(len(tmp),len(peaks))
if len(tmp)<len(peaks):
peaks=tmp
tmp=peaks
dens_df[iidd][mol]=peaks
print(dens_df)
sv_data=dens_df
if j=='gmx_ndx':
tbf.togmxndx(box_res, fld='./'+SYSTEM_NAME, sv_name=SYSTEM_NAME+'_'+i)
if len(j)==3:
if j[0]=='save':
for k in j[1]:
if k=='pkl':
tbf.topickle(fl=sv_data, sv_name='./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+j[2])
if k=='json':
tbf.tojson(fl=sv_data, sv_name='./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+j[2])
if k=='gmx_ndx':
tbf.togmxndx(box_res, fld='./'+SYSTEM_NAME, sv_name=SYSTEM_NAME)
try:
sv_data
except NameError:
pass
else:
for k in j[1]:
if k=='pkl':
tbf.topickle(fl=sv_data, sv_name='./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_'+j[2])
if k=='json':
tbf.tojson(fl=sv_data, sv_name='./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_'+j[2])
###################################################
else:
for j in GROUPS[i]:
if j=='tilt':
tbf.SurfVector_angle(surf,vector)
for jj in j:
if jj[0]!='COMB':
if jj[1]=='surf':
surf=surf[locals()[jj[0]]]
elif jj[1]=='vector':
vector=vector[locals()[jj[0]]]
if jj[0]=='COMB':
if jj[1]=='tilt':
tbf.SurfVector_angle(surf,vector)
###################################################
\ No newline at end of file
import os
import subprocess
import numpy as np
from scipy.signal import find_peaks
from scipy.signal import savgol_filter
import hashlib
#cmd="gmx covar -s test/eq_run.tpr -f test/eq_traj.trr -n test/gmx_ndx/test_\(1\,\ 0\,\ 3\).ndx"
#os.system(cmd)
#cmd="gmx anaeig -v test/eigenvec.trr -s test/eq_run.tpr -f test/eq_traj.trr -n test/gmx_ndx/test_\(1\,\ 0\,\ 3\).ndx -proj -first 1 -last 20"
#os.system(cmd)
#cmd="gmx_clusterByFeatures cluster -s test/eq_run.tpr -f test/eq_traj.trr -n test/gmx_ndx/test_\(1\,\ 0\,\ 3\).ndx \
# -feat proj.xvg -method kmeans -nfeature 20 -cmetric ssr-sst -ncluster 5 \
# -fit2central -sort features -cpdb test/central.pdb \
# -fout test/cluster.xtc -plot test/pca_cluster.png"
#os.system(cmd)
def ndx_index(SYSTEM_NAME):
uniq_id={}
path='./'+SYSTEM_NAME+'/gmx_ndx/'
for f in os.listdir(path):
text=SYSTEM_NAME+'_'+f
iidd=hashlib.md5(text.encode('utf-8')).hexdigest()
uniq_id[iidd]={}
uniq_id[iidd]={'system':SYSTEM_NAME,'domain':f}
return uniq_id
def read_xvg(XVG):
x=[]
y=[]
with open(XVG) as f:
for line in f:
if line.startswith(("#","@")):
continue
else:
cols = line.split()
x.append(cols[0])
y.append(cols[1])
x = np.asarray(x)
y = np.asarray(y)
return x,y
def density_picks(TRR,TPR,IND,SLC,ST,EN,normal,fld,arg,dist_pk=1):
cmd=str(arg)+'\n'+str(arg)+'\n'
if EN==-1:
p = subprocess.Popen(['gmx', 'density', '-symm', '-b', str(ST), '-sl',str(SLC), '-d', normal,
'-f', TRR,
'-s',TPR,
'-n',IND,
'-o', fld+'/tmp.xvg'],
stdin=subprocess.PIPE)
else:
p = subprocess.Popen(['gmx', 'density', '-symm', '-b', str(ST), '-e', str(EN),'-sl',str(SLC), '-d', normal,
'-f', TRR,
'-s',TPR,
'-n',IND,
'-o', fld+'/tmp.xvg'],
stdin=subprocess.PIPE)
p.communicate(cmd.encode('UTF-8'))
p.wait()
x,y=read_xvg(XVG=fld+'/tmp.xvg')
yhat = savgol_filter(y, 15, 4) # window size 15, polynomial order 4
peaks, _ = find_peaks(yhat, distance=dist_pk)
pathname = os.path.abspath(os.path.join(fld, 'tmp.xvg'))
os.remove(pathname)
return peaks #res
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment