Python bindings for Grid3D::points()

resolves #10

Author: mdimura

FlexLabel/include/FlexLabel/FlexLabel.h

@@ -112,8 +112,8 @@ class Grid3D
 		return xyz(tmp);
 	}
 
-	/// @brief returns a coordinate list with all coordinates having
-	/// values > 0.0
+	/// @brief returns a list of cartesian points and corresponding
+	/// weights for cells, that have density values > 0.0
 	Eigen::Matrix4Xf points() const
 	{
 		const int gridSize = grid.size();

FlexLabel/python/LabelLib_pymol.py

@@ -86,38 +86,16 @@ def makeAtom(index, xyz, vdw, name='AV'):
 
 def avToModel(av):
   m = chempy.models.Indexed()
-  
-  area = av.shape[0] * av.shape[1]
-  nx, ny, nz = av.shape
-  ox, oy, oz = av.originXYZ
-  dx = av.discStep
-  g = np.array(av.grid).reshape((nx, ny, nz),order='F')
-  
-  MP=np.array([0.0,0.0,0.0])
-  vol=0.0
-  iat = 0
-  for iz in range(nz):
-    for iy in range(ny):
-      for ix in range(nx):
-        val = g[ix, iy, iz]
-        if val <= 0.0:
-          continue
-        
-        iat += 1
-        resi = int(iat / 10)
-        
-        x = ix * dx + ox
-        y = iy * dx + oy
-        z = iz * dx + oz
-        
-        MP+=np.array([x,y,z])*val
-        vol+=val
-        
-        m.add_atom(makeAtom(iat,[x,y,z],dx * 0.5))
-  
-  MP=MP/vol
-  if iat>0:
-    m.add_atom(makeAtom(iat+1,list(MP),2.0,'AVmp'))
+  points = av.points()
+  r = av.discStep * 0.5
+
+  for i, p in enumerate(points.T):
+    x, y, z, w = p
+    m.add_atom(makeAtom(i+1, [x, y, z], r))
+
+  MP = np.average(points[:3,:],axis=1,weights=points[3])
+  if points.shape[1]>0:
+    m.add_atom(makeAtom(points.shape[1]+1,list(MP),2.0,'AVmp'))
 
   m.update_index()
   return m

FlexLabel/python/bindings.cpp

@@ -60,5 +60,6 @@ PYBIND11_MODULE(LabelLib, m)
 	        .def_readwrite("discStep", &Grid3D::discStep)
 	        .def_readwrite("originXYZ", &Grid3D::originXYZ)
 	        .def_readwrite("shape", &Grid3D::shape)
-	        .def_readwrite("grid", &Grid3D::grid);
+	        .def_readwrite("grid", &Grid3D::grid)
+	        .def("points",&Grid3D::points);
 }

FlexLabel/python/mdtraj_example.py

@@ -77,38 +77,5 @@ def genAV(traj, frameId, attachmentIdx, length, width, radius, resolution, cvFra
 	acv.grid=list(acvarr)
 	return acv
 
-def savePqr(fileName, grid):
-  with open(fileName, "w") as out:
-      area = grid.shape[0] * grid.shape[1]
-      nx, ny, nz = grid.shape
-      ox, oy, oz = grid.originXYZ
-      dx = grid.discStep
-      g = np.array(grid.grid).reshape((nx, ny, nz),order='F')
-      
-      iat = 0
-      for iz in range(nz):
-        for iy in range(ny):
-          for ix in range(nx):
-            val = g[ix, iy, iz]
-            if val <= 0.0:
-              continue
-            
-            iat += 1
-            resi = int(iat / 10)
-            
-            x = ix * dx + ox
-            y = iy * dx + oy
-            z = iz * dx + oz
-            
-            sz = 'ATOM{0: 7}   AV  AV{1: 6}{2:12.1f}{3:8.1f}{4:8.1f}{5:8.2f}{6:7.3f}\n'
-            sz = sz.format(iat, resi, x, y, z, val, dx * 0.5)
-            out.write(sz)
-
-def savePqrFromAtoms(fileName, atoms):
-  sz = 'ATOM{0: 7}    X   X{1: 6}{2:12.1f}{3:8.1f}{4:8.1f}{5:8.2f}{6:7.3f}\n'
-  with open(fileName, "w") as out:
-    for i,at in enumerate(atoms.T):
-      out.write(sz.format(i, i, at[0], at[1], at[2], 1.0, at[3]))
-
 if __name__ == '__main__':
     sys.exit(main(sys.argv[1:]))

FlexLabel/python/usage.py

@@ -2,31 +2,14 @@
 import numpy as np
 
 def savePqr(fileName, grid):
+  points = grid.points()
+  template = 'ATOM{0: 7}   AV  AV{1: 6}{2:12.1f}{3:8.1f}{4:8.1f}{5:8.2f}{6:7.3f}\n'
+  r = grid.discStep * 0.5
   with open(fileName, "w") as out:
-      area = grid.shape[0] * grid.shape[1]
-      nx, ny, nz = grid.shape
-      ox, oy, oz = grid.originXYZ
-      dx = grid.discStep
-      g = np.array(grid.grid).reshape((nx, ny, nz),order='F')
-      
-      iat = 0
-      for iz in range(nz):
-        for iy in range(ny):
-          for ix in range(nx):
-            val = g[ix, iy, iz]
-            if val <= 0.0:
-              continue
-            
-            iat += 1
-            resi = int(iat / 10)
-            
-            x = ix * dx + ox
-            y = iy * dx + oy
-            z = iz * dx + oz
-            
-            sz = 'ATOM{0: 7}   AV  AV{1: 6}{2:12.1f}{3:8.1f}{4:8.1f}{5:8.2f}{6:7.3f}\n'
-            sz = sz.format(iat, resi, x, y, z, val, dx * 0.5)
-            out.write(sz)
+    for i, p in enumerate(points.T):
+      x, y, z, w = p
+      sz = template.format(i + 1, 1, x, y, z, w, r)
+      out.write(sz)
 
 def savePqrFromAtoms(fileName, atoms):
   sz = 'ATOM{0: 7}    X   X{1: 6}{2:12.1f}{3:8.1f}{4:8.1f}{5:8.2f}{6:7.3f}\n'