diff --git a/main.py b/main.py
index 9b94adededdffdb78e64378e70f6dab0d14e49da..a69a97d7a410e61a995aa830a86354e45ec872f1 100644
--- a/main.py
+++ b/main.py
@@ -12,6 +12,7 @@ DISCET=[6, 6, 6]
 NUM_FR=1
 TRAJ=SYSTEM_NAME+'/eq_traj.trr'
 GRO=SYSTEM_NAME+'/eq_final.gro'
+ITP_DIC={'NS':'CER_SOVOVA.itp','FFA':'FFA_CG.itp','CHOL':'CHOL_CG.itp'}
 ###################################################
 # {NAME:[QUEUE OF PROCESSES]}
 #
@@ -26,7 +27,7 @@ GRO=SYSTEM_NAME+'/eq_final.gro'
 #
 #                       if NAME is COMBINATION: [tilt] 
 ###################################################
-GROUPS={'HD_GROUP':['surf',['save', ['pkl','gmx_ndx'],'time_domain_c-normal-cg']],
+GROUPS={'HD_GROUP':['surf',['save', ['pkl','gmx_ndx', 'json'],'time_domain_c-normal-cg']],
         'TL_GROUP':['vector'],
         'COMBINE':['tilt']
 }
@@ -44,6 +45,14 @@ print('================')
 _,data_num,_,res_num,res_type,atom_type,atom_num,_ = tbf.read_gro(GRO)
 print(' ')
 ###################################################
+#Read .itp files
+#weights={}
+#for MOL in ITP_DIC.keys():
+#    weights_tmp = tbf.read_itp(SYSTEM_NAME+'/'+ITP_DIC[MOL])
+#    weights[MOL]=weights_tmp
+#    print(' ')
+#print(weights)
+###################################################
 if os.path.isfile('./'+SYSTEM_NAME+'/'+SYSTEM_NAME+'_data.pkl'):
     if LOG==1:
         print('WARNING: Preprocessing files exist.')
@@ -109,4 +118,5 @@ for i in GROUPS.keys():
     else:
         for j in GROUPS[i]:
             if j=='tilt':
-                tbf.SurfVector_angle(surf,vector)
\ No newline at end of file
+                tbf.SurfVector_angle(surf,vector)
+###################################################
\ No newline at end of file
diff --git a/tooba_f.py b/tooba_f.py
index 0b28dd3af8c7b7f9db9a83141089abe20b361704..77233ba9f186166a687cd0c82a6e837d3a3ca3c4 100644
--- a/tooba_f.py
+++ b/tooba_f.py
@@ -47,6 +47,7 @@ def angle_between3D(v1, v2):
         return 360 - angle
 
 def center_gravity(a):
+    #Todo somatidia diaforetikis mazas -> weighted me maza ????
     m=len(a)
     cg = np.sum(a)/m
     return cg
@@ -221,6 +222,69 @@ def read_gro(gro):
         #print(system,data_num,box_size)
     return system,data_num,box_size,res_num,res_type,atom_type,atom_num,rest_dt
 
+def read_itp(itp):
+    """
+    Read mass and atomtype from .itp file.
+
+    parameters:     itp = [.itp]
+
+    output:         dict {atomtype{mass}}
+
+    ATTENTION:  the format of itp is free but the current function
+                has some restrictions
+                1. sections starting points doesn't have spaces [ex. [atoms] ]
+                2. the sections have to be seperated by an empty line
+    """
+
+    ll=0
+    #mass=[]
+    #atomtype=[]
+    weights={}
+    num_lines = sum(1 for line in open(itp))
+    with open(itp, 'r') as F:
+        bar = Bar('Read .itp', max=num_lines)
+        for line in F:
+
+            if line.strip()=='[atoms]':
+                ll=1
+                #print(line.strip())
+                continue
+            elif line.strip()=='[atomtypes]':
+                ll=2
+                #print(line.strip())
+                continue
+            elif len(line.strip()) == 0:
+                #print(len(line.strip()))
+                ll=0
+                continue
+            elif line.strip()[0] == ';':
+                continue
+
+            if ll==1: #molecule structure itp
+                word=line.split()
+                try:
+                    if float(word[7])==0.0:
+                        print('WARNING: Atomtype mass of '+word[4]+' is 0.0')
+                except:
+                    print('Skip line: '+word[7])
+                    continue
+                #print(word[4],word[7])
+                weights[word[4]]=word[7]
+                #mass.append(word[1]) 
+                #atomtype.append(word[0])
+            if ll==2: #molecule forcefield itp
+                word=line.split()
+                try:
+                    if float(word[1])==0.0:
+                        print('WARNING: Atomtype mass of '+word[0]+' is 0.0')
+                except:
+                    print('Skip line: '+word[1])
+                    continue
+                weights[word[0]]=word[1]
+            bar.next()
+        bar.finish()
+    return weights
+
 def domain_decomposition(data,dx,dy,dz):
     """
     Identify subdomain for each frame of 
@@ -529,14 +593,24 @@ def coord2norm2cg(coord_vector,img=True):
         bar = Bar('Step: '+step, max=len(coord_vector[step].keys()))
         for subdomain in coord_vector[step].keys():
             surf[step][subdomain]={}
+            xx=[]
+            yy=[]
+            zz=[]
             #surf_points=[]
             res_cg={}
             cnt=0
             for resid in coord_vector[step][subdomain].keys():
                 cnt=cnt+1
-                cgx = center_gravity(coord_vector[step][subdomain][resid])
-                cgy = center_gravity(coord_vector[step][subdomain][resid])
-                cgz = center_gravity(coord_vector[step][subdomain][resid])
+                for ii in coord_vector[step][subdomain][resid]:
+                    #print(ii)
+                    #print(ii[0],ii[1],ii[2])
+                    #print(coord_vector[step][subdomain][resid][0])
+                    xx.append(ii[0])
+                    yy.append(ii[1])
+                    zz.append(ii[2])
+                cgx = center_gravity(xx)
+                cgy = center_gravity(yy)
+                cgz = center_gravity(zz)
                 res_cg[resid]=[cgx,cgy,cgz]
                 if cnt==1:
                     surf_points=coord_vector[step][subdomain][resid]