From f8332608bd434f7808a652d2192e137ce3b67bd6 Mon Sep 17 00:00:00 2001
From: skarozis <skarozis@ipta.demokritos.gr>
Date: Thu, 2 Jul 2020 17:04:57 +0300
Subject: [PATCH] Initial support for gmx commands

---
 CHANGELOG    | 11 +++++++
 main.py      | 90 ++++++++++++++++++++++++++++++++++++++++------------
 tooba_gmx.py | 72 +++++++++++++++++++++++++++++++++++++++++
 3 files changed, 153 insertions(+), 20 deletions(-)
 create mode 100644 tooba_gmx.py

diff --git a/CHANGELOG b/CHANGELOG
index 58a0195..73eb59c 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -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
diff --git a/main.py b/main.py
index a69a97d..2ecb5f4 100644
--- a/main.py
+++ b/main.py
@@ -1,5 +1,6 @@
 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
diff --git a/tooba_gmx.py b/tooba_gmx.py
new file mode 100644
index 0000000..325f91c
--- /dev/null
+++ b/tooba_gmx.py
@@ -0,0 +1,72 @@
+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
-- 
2.24.1