memory opt + debugging

CHANGELOG

@@ -6,6 +6,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.0.8] - 2020-07-16
+### Added
+- memory optimizations
+- Add check of molecule existance in .ndx file
+
+### Changed
+- Update I/O to functions
+
+### Removed
+- Remove log message option
+
 ## [0.0.7] - 2020-07-03
 ### Added
 - unified properties output
@@ -18,7 +29,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Removed
 - None
 
-
 ## [0.0.6] - 2020-07-02
 ### Added
 - Initial support for gmx commands

main.py

@@ -7,10 +7,9 @@ import tooba_gmx as tbgmx
 #           in case of tail is the one closest to the head, hence
 #           the code is a good approximation 
 ###################################################
-LOG=0
 SYSTEM_NAME='test'
-DISCET=[6, 6, 6]
-NUM_FR=1
+DISCET=[3.5, 3.5, 3.5]
+NUM_FR=100
 TRAJ=SYSTEM_NAME+'/eq_traj.trr'
 GRO=SYSTEM_NAME+'/eq_final.gro'
 TPR=SYSTEM_NAME+'/eq_run.tpr'
@@ -37,7 +36,7 @@ TPR=SYSTEM_NAME+'/eq_run.tpr'
 ###################################################
 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'],
+        'TL_GROUP':['vector',['save', ['pkl'],'vec']],
         'COMBINE':[['HD_GROUP','surf'],['TL_GROUP','vector'],['COMB','tilt'],['save', ['pkl'],'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']}
@@ -56,6 +55,14 @@ _,data_num,_,res_num,res_type,atom_type,atom_num,_ = tbf.read_gro(GRO)
 print(' ')
 ###################################################
 #--------------------------------------------------
+#Count frames & calculate time of calculations
+cnt_fr,time_index=tbf.count_frames(TRAJ,True)
+MAX_FR=list(time_index.keys())[-1]
+#print(MAX_FR,time_index[MAX_FR])
+ST_FR=MAX_FR-NUM_FR
+#print(ST_FR,time_index[ST_FR])
+###################################################
+#--------------------------------------------------
 #Read .itp files
 #weights={}
 #for MOL in ITP_DIC.keys():
@@ -66,15 +73,12 @@ print(' ')
 ###################################################
 #--------------------------------------------------
 if os.path.isfile('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_data.pkl'):
-    if LOG==1:
-        print('WARNING: Preprocessing files exist.')
-        print('         Erase data.pkl if the system is new.')
-        print('--------------------------------------------')
-    data_all=tbf.frompickle('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_data.pkl')
+    pass
 else:
     #Read .trr file
     data_all=tbf.fr_export(trajfile=TRAJ,num_frames=NUM_FR)
     tbf.topickle(fl=data_all, sv_name='./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_data')
+    del data_all
 ###################################################
 #--------------------------------------------------
 #Check save files if exist in order to skip functions
@@ -112,24 +116,18 @@ for i in GROUPS.keys():
 #not COMBINE section
     if i!='COMBINE': 
         if os.path.isfile('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_ndx.pkl'):
-            if LOG==1: 
-                print('WARNING: Preprocessing files exist.')
-                print('         Erase ndx_HEAD.pkl if the system is new.')
-                print('--------------------------------------------')
-            group_ndx=tbf.frompickle('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_ndx.pkl')
+            pass
         else:
             #Find atom type index in lists created above
             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')
+            del group_ndx
 #--------------------------------------------------
         if os.path.isfile('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_box.pkl'):
-            if LOG==1:
-                print('WARNING: Preprocessing files exist.')
-                print('         Erase box.pkl if the system is new')
-                print('         or new grid is applied !')
-                print('--------------------------------------------')
-            box_res=tbf.frompickle('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_box.pkl')
+            pass            
         else:
+            data_all=tbf.frompickle('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_data.pkl')
+            group_ndx=tbf.frompickle('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_ndx.pkl')
             #Create subdomains coordinates
             box_p=tbf.domain_decomposition(data=data_all,dx=DISCET[0],dy=DISCET[1],dz=DISCET[2])
             #Assign desired atoms (from above function) to subdomains
@@ -137,14 +135,32 @@ for i in GROUPS.keys():
             ##result2: {step:{res:{atom_type:{(subX,subYsubZ):[atom_num]}}}}
             _,box_res=tbf.atom2grid(data_all,box_p, group_ndx)
             tbf.topickle(fl=box_res, sv_name='./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_box')
+            del data_all
+            del group_ndx
+            del box_p
+            del box_res
+###################################################
+        if os.path.isfile('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_FR'+str(NUM_FR)+'.pkl'):
+            pass   
+        else:
+            #Creates dictionary with coordinates per subdomain for each frame
+            data_all=tbf.frompickle('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_data.pkl')
+            box_res=tbf.frompickle('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_box.pkl')
+            _,coord_vector=tbf.sub_coord(box=box_res, data=data_all, res_num=res_num)
+            tbf.topickle(fl=coord_vector, sv_name='./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_FR'+str(NUM_FR))
+            del data_all
+            del box_res
+            del coord_vector
 ###################################################
-        #Creates dictionary with coordinates per subdomain for each frame
-        _,coord_vector=tbf.sub_coord(box=box_res, data=data_all, res_num=res_num)
         for j in GROUPS[i]:
             if len(j) > 1:
                 if j=='surf' and sv_index[j]['status']=='not exist':
+                    if j not in locals():
+                        surf={}
                     #Creates dictionary with c, normal per subdomain for each frame
+                    coord_vector=tbf.frompickle('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_FR'+str(NUM_FR)+'.pkl')
                     surf[i]=tbf.coord2norm2cg(coord_vector,img=False)
+                    del coord_vector
                     sv_data=surf[i]
                 elif j=='surf' and sv_index[j]['status']=='exist':
                     if j not in locals():
@@ -154,7 +170,9 @@ for i in GROUPS.keys():
                 if j=='vector' and sv_index[j]['status']=='not exist':
                     if j not in locals():
                         vector={}
+                    coord_vector=tbf.frompickle('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_FR'+str(NUM_FR)+'.pkl')
                     vector[i]=tbf.coord2vector(coord_vector)
+                    del coord_vector
                     sv_data=vector[i]
                 elif j=='vector' and sv_index[j]['status']=='exist':
                     if j not in locals():
@@ -171,12 +189,15 @@ for i in GROUPS.keys():
                     dens_df={}
                     for iidd in uniq_id.keys():
                         dens_df[iidd]={}
-                        fl='./'+uniq_id[iidd]['system']+'/gmx_ndx/'+uniq_id[iidd]['ndx_file']
+                        fl='./'+uniq_id[iidd]['system']+'/gmx_ndx/'+uniq_id[iidd]['ndx_file']                  
                         cnt=-1
                         for mol in locals()[i].keys():
+                            ind=tbf.search_pattern(fl,mol)
+                            if ind=='not exist':
+                                break
                             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)
+                                peaks = tbgmx.density_picks(TRR=TRAJ,TPR=TPR,IND=fl,SLC=400,ST=time_index[ST_FR],EN=-1,normal=d,fld='./'+uniq_id[iidd]['system'],arg=cnt,dist_pk=20)
                                 if d=='x':
                                     tmp=peaks
                                 else:
@@ -191,7 +212,9 @@ for i in GROUPS.keys():
                     dens_df=tbf.frompickle(sv_index[j]['name'])
 #--------------------------------------------------
                 if j=='gmx_ndx':
+                    box_res=tbf.frompickle('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_'+i+'_box.pkl')
                     tbf.togmxndx(box_res, fld='./'+SYSTEM_NAME, sv_name=SYSTEM_NAME+'_'+i)
+                    del box_res
 #--------------------------------------------------
             # Save module
             if len(j)==3:
@@ -214,16 +237,20 @@ for i in GROUPS.keys():
             if j[0]!='COMB':
                 if j[1]=='surf':
                     surf=surf[j[0]]
+                    #del surf[j[0]]
                 if j[1]=='vector':
                     vector=vector[j[0]]
+                    #del vector[j[0]]
             #Calculate COMBINE property
             if j[0]=='COMB':
                 if j[1]=='tilt' and sv_index[str(j)]['status']=='not exist':
                     tilt=tbf.SurfVector_angle(surf,vector)
+                    del surf
+                    del vector
                     #Loop over timesteps and keep avgs tilts for each step
                     avg={}
-                    ss=[]
                     for step in tilt.keys():
+                        ss=[]
                         for sub in tilt[step].keys():
                             avgs=tilt[step][sub]['avg/frame']
                             if sub not in ss:
@@ -231,6 +258,7 @@ for i in GROUPS.keys():
                                 avg[sub]=avgs
                             else:
                                 avg[sub].append(avgs)
+                    del tilt
                     #Calculate total average 
                     tot_avg={}
                     for sub in avg.keys():

tooba_f.py

@@ -5,6 +5,7 @@ import matplotlib.pyplot as plt
 import scipy.optimize
 import pickle as pkl
 import json
+import re
 from progress.bar import Bar
 
 from pytrr import (
@@ -13,6 +14,18 @@ from pytrr import (
     skip_trr_data,
     )
 
+def search_pattern(fl,pattern):
+    textfile = open(fl, 'r')
+    filetext = textfile.read()
+    textfile.close()
+    matches = re.findall(pattern, filetext)
+    
+    if len(matches)==0:
+        ind='not exist'
+    else:
+        ind='exist'
+    return ind
+
 def fitPlaneLTSQ(XYZ):
     #Source: https://gist.github.com/RustingSword/e22a11e1d391f2ab1f2c
     (rows, cols) = XYZ.shape
@@ -123,23 +136,30 @@ def points2vector(data):
 
     return vv[0]
 
-def count_frames(trajfile='traj.trr'):
+def count_frames(trajfile='traj.trr', pr_time=False):
     """
     Count total frames of .trr file
 
     parameters:     trajfile = [.trr]
     """
+    time_index={}
     cnt_fr=0
     with open(trajfile, 'rb') as inputfile:
-        for i in range(1000):
+        for i in range(1000000):
             try: 
                 header = read_trr_header(inputfile)
-                #print('Step: {step}, time: {time}'.format(**header))
                 skip_trr_data(inputfile, header)
+                if pr_time==True:
+                    #print('Step: {step}, time: {time}'.format(**header))
+                    #print('Step: {step}, time: {time}'.format(**header))
+                    time_index[cnt_fr]='{time}'.format(**header)
                 cnt_fr=cnt_fr+1
             except EOFError:
                 pass
-    return cnt_fr
+    if pr_time==True:
+        return cnt_fr,time_index
+    else:
+        return cnt_fr
 
 def fr_export(trajfile='traj.trr', num_frames=1):
     """
@@ -616,10 +636,11 @@ def coord2norm2cg(coord_vector,img=True):
                     surf_points=coord_vector[step][subdomain][resid]
                 else:
                     surf_points=surf_points+coord_vector[step][subdomain][resid]
+
             c,normal = fitPlaneLTSQ(np.array(surf_points))
             surf[step][subdomain]={'c':c, 'normal':normal}
             surf[step][subdomain].update(res_cg)
-    
+            
             if img==True:
                 fl_save='png/'
                 save=fl_save+str(subdomain)
@@ -643,7 +664,7 @@ def coord2vector(coord_vector):
     output:             dictionary = {step:{subdomain:{resid:[x, y, z]}}
     """
     print(' ')
-    print('Finding surface progress:')
+    print('Finding vector progress:')
     print('==============================')
     vector={}
     for step in coord_vector.keys():
@@ -652,7 +673,10 @@ def coord2vector(coord_vector):
         for subdomain in coord_vector[step].keys():
             vector[step][subdomain]={}
             for resid in coord_vector[step][subdomain].keys():
-                vv = points2vector(np.array(coord_vector[step][subdomain][resid]))
+                grp=[]
+                for ii in coord_vector[step][subdomain][resid]:
+                    grp.append(ii)
+                vv = points2vector(np.array(grp))
                 vector[step][subdomain][resid]=vv
             bar.next()
         bar.finish()
@@ -676,6 +700,11 @@ def SurfVector_angle(surf,vector):
         angle[step]={}
         bar = Bar('Step: '+step, max=len(surf[step].keys()))
         for sudomain in surf[step].keys():
+            try:
+                vector[step][sudomain].keys()
+            except KeyError:
+                #print(step, sudomain)
+                continue
             angle[step][sudomain]={}
             tot=[]
             for resid in vector[step][sudomain].keys():