fix cube file grid ordering bug, resize_grid point_tree bug & Bq atom validation

The meshgrid in the density code path used default 'xy' indexing which
produced (z-slow, x-mid, y-fast) ordering, mismatching the cube file's
(x-slow, y-mid, z-fast) density storage. This scrambled occupied grid
point positions, giving correct Mol_Vol but wrong %V_Bur. Fixed by
using indexing='ij' with .reshape(3,-1).T.

Also fixed resize_grid not returning an updated point_tree (sphere
volume error ~6% for cube files), and allowed Bq ghost atoms through
the VDW radii validation check so --sambvca works when atom1 is H.

Adds benzene and Ne test cube files at coarse/medium/fine resolutions.

Author: bobbypaton

dbstep/Dbstep.py

@@ -28,7 +28,6 @@ def read_input(molecule, ext, options):
 		DataParser object with parsed molecule data to be used by the rest of the program
 	"""
 	if ext == ".cube":
-		options.surface = "density"
 		mol = CubeParser(molecule, "cube")
 	else:
 		if ext in [".xyz", ".com", ".gjf"]:

dbstep/parse_data.py

@@ -143,17 +143,11 @@ def occupied(grid, coords, radii, origin, options):
 	for n in range(len(coords)):
 		center = coords[n] + origin
 		idx.append(point_tree.query_ball_point(center, radii[n], workers=-1))
-	# construct a list of indices of the grid array that are occupied / unoccupied
+	# construct a list of indices of the grid array that are occupied
 	jdx = [y for x in idx for y in x]
-	if options.qsar:
-		kdx = [i for i in range(len(grid)) if i not in jdx]
 
 	# removes duplicates since a voxel can only be occupied once
 	jdx = list(set(jdx))
-	if options.qsar:
-		kdx = list(set(kdx))
-		onehot = np.zeros(len(grid))
-		onehot[jdx] = 1.0
 
 	if options.verbose:
 		print("   There are {} occupied grid points.".format(len(jdx)))
@@ -167,13 +161,8 @@ def occupied(grid, coords, radii, origin, options):
 
 		u = pptk.viewer(grid)
 		v = pptk.viewer(grid[jdx])
-		if options.qsar:
-			w = pptk.viewer(grid[kdx])
 
-	if options.qsar:
-		return grid[jdx], grid[kdx], onehot, point_tree
-	else:
-		return grid[jdx], point_tree, occ_vol
+	return grid[jdx], point_tree, occ_vol
 
 
 def occupied_dens(grid, dens, options):
@@ -239,8 +228,9 @@ def resize_grid(x_max, y_max, z_max, x_min, y_min, z_min, options, mol):
 	y_vals = np.linspace(y_min, y_max, mol.ydim)
 	z_vals = np.linspace(z_min, z_max, mol.zdim)
 	grid = np.array(np.meshgrid(x_vals, y_vals, z_vals)).T.reshape(-1, 3)
+	point_tree = spatial.cKDTree(grid, balanced_tree=False, compact_nodes=False)
 
-	return grid
+	return grid, point_tree
 
 
 def get_classic_sterimol(coords, radii, atoms):

dbstep/sterics.py

@@ -0,0 +1,3 @@
+1
+Ne
+Ne   1.219400   -0.165200    2.160000 

tests/cube_files/Ne.xyz

@@ -0,0 +1,14 @@
+12
+benzene
+H  -2.426483  0.480591  -0.000006
+C  -1.364235  0.27023  0.000002
+C  -0.916128  -1.046358  0.000005
+H  -1.62943  -1.861106  0.000002
+C  0.448089  -1.31655  -0.000008
+H  0.797011  -2.341683  -0.00002
+C  1.364252  -0.270207  0.000012
+H  2.426484  -0.480624  -0.000005
+C  0.916127  1.046345  -0.000005
+H  1.629476  1.861062  0.000001
+C  -0.448107  1.316554  0.000001
+H  -0.797041  2.341678  -0.000018