diff --git a/core/PhysiCell_cell.cpp b/core/PhysiCell_cell.cpp
index 1f3e2646b..2755d7d5d 100644
--- a/core/PhysiCell_cell.cpp
+++ b/core/PhysiCell_cell.cpp
@@ -609,22 +609,25 @@ Cell* Cell::divide( )
 	rand_vec *= radius; // multiply direction times the displacement 
 	*/
 
+	// direction to displace daughter cell
 	std::vector<double> rand_vec; 
 	if( this->functions.cell_division_direction_function )
 	{ 
 		rand_vec = this->functions.cell_division_direction_function( this ); 
         if( default_microenvironment_options.simulate_2D == true )  // ensure vec in XY plane
 	        { rand_vec[2] = 0.0; }
-        normalize( &rand_vec );  // ensure normalized
 	}
 	else
 	{
 		rand_vec = cell_division_orientation(); 
+		// make it orthogonal to the cell's orientation
 		rand_vec = rand_vec- phenotype.geometry.polarity*(rand_vec[0]*state.orientation[0]+ 
 			rand_vec[1]*state.orientation[1]+rand_vec[2]*state.orientation[2])*state.orientation;	
 	}
-
-	rand_vec *= phenotype.geometry.radius;
+	// make sure it is a unit vector
+	normalize( &rand_vec );
+	// push the cells far enough apart to avoid strong overlap, but still maintain some overlap
+	rand_vec *= 0.5 * phenotype.geometry.radius;
 
 	child->assign_position(position[0] + rand_vec[0],
 						   position[1] + rand_vec[1],
@@ -632,8 +635,8 @@ Cell* Cell::divide( )
 						 
 	//change my position to keep the center of mass intact 
 	// and then see if I need to update my voxel index
-	static double negative_one_half = -0.5; 
-	axpy( &position, negative_one_half , rand_vec ); // position = position - 0.5*rand_vec; 
+	static double negative_one = -1.0; 
+	axpy( &position, negative_one , rand_vec ); // position = position - rand_vec; 
 
 	//If this cell has been moved outside of the boundaries, mark it as such.
 	//(If the child cell is outside of the boundaries, that has been taken care of in the assign_position function.)