diff --git a/CHANGELOG b/CHANGELOG index f4f9c5c3c33d87ccb2d118be0c7a58b4a94f778b..9cbfa481b76d0c8383ee634b31f7d267f44cac28 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.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 diff --git a/main.py b/main.py index 315feaf1ba7a79ed588c4dbb653e1b4a0a50b8fb..27e1f2cbcfd98008238f3b5f6be8a73c054b336a 100644 --- a/main.py +++ b/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(): diff --git a/tooba_f.py b/tooba_f.py index c7b8e4f85394a731619bd5c686dd41f0c3fa3954..6cb6f51e72afaa5a61fcbbf9a5181732dff5a0bd 100644 --- a/tooba_f.py +++ b/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():