Read gro file

CHANGELOG

@@ -12,6 +12,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - read all or specific part of .trr files
 - calculate surface from points
 - calculate angle between two vectors(3D) or lines(2D)
+- plot points and surface
+- read .gro for determine residue type and atom type 
 
 ### Changed
 - None

README.md

 # TooBBA
 
-# Tool for Bilayer in Blocks Analysis (TooBBA)
-It is an tool that reads .gro files, perform a domain decomposition (user defined) and analyzes each block in terms of lipids distribution. 
\ No newline at end of file
+## 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.

main.py

 import tooba_f as tbf
 
-P1=(0,1,0)
-P2=(1,1,0)
-print(tbf.angle_between3D(P1,P2))
+#P1=(0,1,0)
+#P2=(1,1,0)
+#print(tbf.angle_between3D(P1,P2))
 
 
-tbf.last_fr_export(trajfile='traj.trr')
-print('test')
\ No newline at end of file
+#tbf.last_fr_export(trajfile='traj.trr')
+tbf.read_gro_types('initial.gro')
\ No newline at end of file

requirements.txt

+numpy 
+scipy
+mpl_toolkits
+matplotlib
+pytrr
\ No newline at end of file

tooba_f.py

@@ -4,6 +4,12 @@ import scipy.optimize
 from mpl_toolkits.mplot3d import Axes3D
 import matplotlib.pyplot as plt
 
+from pytrr import (
+    read_trr_header,
+    read_trr_data,
+    skip_trr_data,
+    )
+
 def fitPlaneLTSQ(XYZ):
     (rows, cols) = XYZ.shape
     G = np.ones((rows, 3))
@@ -36,7 +42,7 @@ def angle_between3D(v1, v2):
     else:
         return 360 - angle
 
-def plot_surf(data, normal):
+def plot_surf(data, normal, c):
     #Plot surface
     fig = plt.figure()
     ax = fig.gca(projection='3d')
@@ -66,12 +72,6 @@ def plot_surf(data, normal):
     plt.show()
 
 def count_frames(trajfile='traj.trr'):
-    from pytrr import (
-        read_trr_header,
-        read_trr_data,
-        skip_trr_data,
-        )
-
     cnt_fr=0
     with open(trajfile, 'rb') as inputfile:
         for i in range(1000):
@@ -84,15 +84,10 @@ def count_frames(trajfile='traj.trr'):
                 pass
     return cnt_fr
 
-def last_fr_export(trajfile='traj.trr'):
-    from pytrr import (
-        read_trr_header,
-        read_trr_data,
-        skip_trr_data,
-        )
+def fr_export(trajfile='traj.trr', num_frames=1):
     cnt_fr=count_frames(trajfile)
     with open(trajfile, 'rb') as inputfile:
-        for i in range(cnt_fr-1):
+        for i in range(cnt_fr-num_frames):
             header = read_trr_header(inputfile)
             #print('Step: {step}, time: {time}'.format(**header))
             skip_trr_data(inputfile, header)
@@ -103,34 +98,24 @@ def last_fr_export(trajfile='traj.trr'):
         print(data['box'])
         print(data['x'][0])
 
-#Random data
-#data = np.random.randn(100, 3)/3
-#data[:, 2] /=10
-
-#Find surface of points-> for lipid plane
-#c, normal = fitPlaneLTSQ(data)
-#print('Normal of surface is: ',normal)
-
-#Angle between two lines-> for finding tilt
-#P1=(normal[0],normal[1], normal[2])
-#P2=(0,0,-0.99998042)
-
-#print('Line is: ',P2)
-#print('Angle between vectors is: ',angle_between3D(P1,P2))
-
-#Fit line to points->for lipid chain
-#from scipy import stats
-#slope, intercept, r_value, p_value, std_err = stats.linregress(data[:,0],data[:,1])
-#print('y=',slope,'x+',intercept)
-#print('r_value',r_value)
-
-#Angle between two lines-> for finding tilt
-#P1=(normal[0],normal[1])
-#P2=(1,1)
-
-#print('Line is: ',P2)
-#print('Angle between lines is: ',angle_between2D(P1,P2))
-
-#plot_surf(data=data, normal=normal)
-
-
+def read_gro_types(gro):
+    cnt=0
+    data_num=0
+    with open(gro, 'r') as F:
+        for line in F:
+            cnt=cnt+1
+            print(cnt)
+            if cnt>2:
+                res_num = line[:5]  
+                res_type = line[5:10]  
+                atom_type = line[10:15] 
+                atom_num  = line[15:20]
+                rest_dt = line[20:]
+                print(res_num,res_type,atom_type,atom_num,rest_dt)
+            elif cnt==1:
+                system=line[:10]
+            elif cnt==2:
+                data_num=int(line[:7])
+            if cnt>data_num:
+                box_size=line[:50]
+        print(system,data_num,box_size)
\ No newline at end of file