Add order parameter calculation, correction in whole workflow etc

CHANGELOG

@@ -6,6 +6,20 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.0.9] - 2020-07-30
+### Added
+- calculate rdf profile for all combination of molecules
+- create ndx for order parameter calculation
+- calculate order parameter properties
+- debug save function
+
+### Changed
+- text that is converted to unique hashed id, is now independent of group name and filename
+- combine index and gmx ndx function
+
+### Removed
+- None
+
 ## [0.0.8] - 2020-07-17
 ### Added
 - memory optimizations

README.md

@@ -2,4 +2,4 @@
 
 ## Tool for Bilayer in Blocks Analysis (TooBBA)
 
-It is an tool that reads .trr files, perform a domain decomposition (user defined) and analyzes each block in terms of lipids distribution.
+It is an tool that reads .trr files, perform a domain decomposition (user defined) and calculates structural properties for each block (domain). The final output constitute a domain dataset, aiming to feed Machine Learning (ML) algorithms and identify coherent groups of observations, sharing the same characteristics and properties.

tooba_f.py

@@ -721,10 +721,10 @@ def SurfVector_angle(surf,vector):
 
 def togmxndx(box_res, fld, sv_name):
     print(' ')
-    print('Save to gmx_ndx progress:')
+    print('Save to ndx progress:')
     print('==============================')
     cnt=0
-    fl_save=fld+'/gmx_ndx/'
+    fl_save=fld+'/'
     if not os.path.exists(fl_save):
         os.makedirs(fl_save)
     else:

tooba_gmx.py

@@ -4,6 +4,8 @@ import numpy as np
 from scipy.signal import find_peaks
 from scipy.signal import savgol_filter
 import hashlib
+from progress.bar import Bar
+
 #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)
 
@@ -17,15 +19,18 @@ import hashlib
 
 
 #os.system(cmd)
-def ndx_index(SYSTEM_NAME):
+def ndx_index(SYSTEM_NAME, pth):
     uniq_id={}
-    path='./'+SYSTEM_NAME+'/gmx_ndx/'
+    path='./'+SYSTEM_NAME+'/'+pth+'/'
     for f in os.listdir(path):
-        text=SYSTEM_NAME+'_'+f
+        dom=f.split('(')[1].split(')')[0]
+        dom='('+dom+')'
+        #dom=f.split('_')[2].split('.')[0]
+        text=SYSTEM_NAME+'_'+dom
         iidd=hashlib.md5(text.encode('utf-8')).hexdigest()
-        dom=f.split('_')[2].split('.')[0]
+        #dom=f.split('_')[2].split('.')[0]
         uniq_id[iidd]={}
-        uniq_id[iidd]={'system':SYSTEM_NAME,'ndx_file':f,'domain':dom}
+        uniq_id[iidd]={'system':SYSTEM_NAME,'fld':pth,'ndx_file':f,'domain':dom}
     return uniq_id
 
 def read_xvg(XVG):
@@ -43,7 +48,7 @@ def read_xvg(XVG):
     y = np.asarray(y)
     return x,y
 
-def density_picks(TRR,TPR,IND,SLC,ST,EN,normal,fld,arg,dist_pk=1):
+def density_peaks(TRR,TPR,IND,SLC,ST,EN,normal,fld,arg,dist_pk=1):
 
     cmd=str(arg)+'\n'+str(arg)+'\n'
     if EN==-1:
@@ -54,7 +59,7 @@ def density_picks(TRR,TPR,IND,SLC,ST,EN,normal,fld,arg,dist_pk=1):
                             '-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,
+        p = subprocess.Popen(['gmx', 'density', '-symm', '-b', str(ST), '-e', str(EN), '-sl',str(SLC), '-d', normal,
                             '-f', TRR,
                             '-s',TPR,
                             '-n',IND,
@@ -70,4 +75,119 @@ def density_picks(TRR,TPR,IND,SLC,ST,EN,normal,fld,arg,dist_pk=1):
     pathname = os.path.abspath(os.path.join(fld, 'tmp.xvg'))
     os.remove(pathname)
 
-    return peaks #res
\ No newline at end of file
+    return peaks
+
+def rdf_peaks(TRR,TPR,IND,ST,EN,fld,arg1,arg2,dist_pk):
+
+    cmd=str(arg1)+'\n'+str(arg2)+'\n'
+    if EN==-1:
+        p = subprocess.Popen(['gmx', 'rdf', '-b', str(ST), '-selrpos', 'mol_com',
+                            '-f', TRR,
+                            '-s',TPR,
+                            '-n',IND,
+                            '-o', fld+'/tmp.xvg'],
+                        stdin=subprocess.PIPE)
+    else:
+        p = subprocess.Popen(['gmx', 'rdf', '-b', str(ST), '-e', str(EN), '-selrpos', 'mol_com',
+                            '-f', TRR,
+                            '-s',TPR,
+                            '-n',IND,
+                            '-o', fld+'/tmp.xvg'],
+                        stdin=subprocess.PIPE)
+    p.communicate(cmd.encode('UTF-8'))
+    p.send_signal('<Ctrl>-D')
+    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
+
+def order_ndx(box_res,fld,atoms,sv_name):
+#def togmxndx(box_res, fld, sv_name):
+    print(' ')
+    print('Save to order_ndx progress:')
+    print('==============================')
+    cnt=0
+    fl_save=fld+'/'
+    if not os.path.exists(fl_save):
+        os.makedirs(fl_save)
+    else:
+        print('WARNING: .ndx files exists. Nothing to do!')
+        return
+ 
+    ndx={}
+    for step in box_res.keys():
+        if cnt==1:
+            break
+        for resid in box_res[step].keys():
+            for res in box_res[step][resid].keys():
+                for subdomain in box_res[step][resid][res].keys():
+                    ndx[subdomain]={}
+
+        rs=[]
+        bar = Bar('Create format: ', max=len(box_res[step].keys()))
+        for resid in box_res[step].keys():
+            for res in box_res[step][resid].keys():
+                for subdomain in box_res[step][resid][res]:
+                    posl=-1
+                    for at in atoms:
+                        posl=posl+1
+                        if (at,subdomain) not in rs:
+                            rs.append((at,subdomain))
+                            ndx[subdomain][at]=[]
+                        posf=-1
+                        for atomnum in box_res[step][resid][res][subdomain]:    
+                            posf=posf+1
+                            if posl==posf:
+                                ndx[subdomain][at].append(atomnum)
+            bar.next()
+        bar.finish()                 
+        cnt=cnt+1
+
+    bar = Bar('Save file: ', max=len(ndx.keys()))
+    for subdomain in ndx.keys():
+        save=fl_save+sv_name+'_'+str(subdomain)+'.ndx'
+        with open(save, "a") as myfile:
+            for res in ndx[subdomain].keys():
+                myfile.write("["+res+"]\n")
+                for atomnum in ndx[subdomain][res]:
+                    myfile.write(atomnum+'\n')
+
+        bar.next()
+    bar.finish()
+
+def order(TRR,TPR,IND,ST,EN,normal,fld,dist_pk=1):
+
+    if EN==-1:
+        p = subprocess.Popen(['gmx', 'order', '-b', str(ST), '-d', normal,
+                            '-f', TRR,
+                            '-s',TPR,
+                            '-n',IND,
+                            '-o', fld+'/tmp1.xvg',
+                            '-od', fld+'/tmp2.xvg'],
+                        stdin=subprocess.PIPE)
+    else:
+        p = subprocess.Popen(['gmx', 'order', '-b', str(ST), '-e', str(EN), '-d', normal,
+                            '-f', TRR,
+                            '-s',TPR,
+                            '-n',IND,
+                            '-o', fld+'/tmp1.xvg',
+                            '-od', fld+'/tmp2.xvg'],
+                        stdin=subprocess.PIPE)
+    p.wait()
+
+    x,y=read_xvg(XVG=fld+'/tmp1.xvg')
+    #yhat = savgol_filter(y, 15, 4) # window size 15, polynomial order 4
+    #peaks, _ = find_peaks(yhat, distance=dist_pk)
+    #exit()
+    pathname = os.path.abspath(os.path.join(fld, 'tmp1.xvg'))
+    os.remove(pathname)
+    pathname = os.path.abspath(os.path.join(fld, 'tmp2.xvg'))
+    os.remove(pathname)
+
+    return y
\ No newline at end of file