Commit f79a7104 authored by Leonardo Solis's avatar Leonardo Solis
Browse files

corrections + code-restructure in grad

parent d266192b
......@@ -198,8 +198,8 @@ odock: check-env-all stringify $(SRC)
# Example
PDB := 3ce3
NRUN := 100
POPSIZE := 500
NRUN := 200
POPSIZE := 2048
TESTNAME:= test
test: odock
......
......@@ -172,7 +172,6 @@ void gpu_perform_elitist_selection(
for (gene_counter = get_local_id(0);
gene_counter < dockpars_num_of_genes;
gene_counter+= NUM_OF_THREADS_PER_BLOCK) {
dockpars_conformations_next[GENOTYPE_LENGTH_IN_GLOBMEM*get_group_id(0)+gene_counter] = dockpars_conformations_current[GENOTYPE_LENGTH_IN_GLOBMEM*get_group_id(0) +
GENOTYPE_LENGTH_IN_GLOBMEM*best_ID[0]+gene_counter];
dockpars_conformations_next[GENOTYPE_LENGTH_IN_GLOBMEM*get_group_id(0)+gene_counter] = dockpars_conformations_current[GENOTYPE_LENGTH_IN_GLOBMEM*get_group_id(0) + GENOTYPE_LENGTH_IN_GLOBMEM*best_ID[0]+gene_counter];
}
}
......@@ -2,11 +2,7 @@
// Implementation of auxiliary functions
// for the gradient-based minimizer
bool is_gradDescent_enabled(
__local bool* is_gNorm_gt_gMin,
__local bool* is_nIter_lt_maxIter,
__local bool* is_perturb_gt_gene_min,
__local bool* is_perturb_gt_genotype,
__local float* local_gNorm,
float gradMin_tol,
__local uint* local_nIter,
......@@ -16,10 +12,14 @@ bool is_gradDescent_enabled(
__local bool* is_gradDescentEn,
uint gradMin_numElements)
{
bool is_gNorm_gt_gMin;
bool is_nIter_lt_maxIter;
bool is_perturb_gt_genotype;
if (get_local_id(0) == 0) {
*is_gNorm_gt_gMin = (local_gNorm[0] >= gradMin_tol);
*is_nIter_lt_maxIter = (local_nIter[0] <= gradMin_maxiter);
*is_perturb_gt_genotype = true;
is_gNorm_gt_gMin = (local_gNorm[0] >= gradMin_tol);
is_nIter_lt_maxIter = (local_nIter[0] <= gradMin_maxiter);
is_perturb_gt_genotype = true;
}
// For every gene, let's determine
......@@ -35,153 +35,35 @@ bool is_gradDescent_enabled(
if (get_local_id(0) == 0) {
// Reduce all is_perturb_gt_gene_min's
// into their corresponding genotype
for(uint i = 0;
i < gradMin_numElements;
i++) {
*is_perturb_gt_genotype = *is_perturb_gt_genotype && is_perturb_gt_gene_min[i];
is_perturb_gt_genotype = is_perturb_gt_genotype && is_perturb_gt_gene_min[i];
/*
printf("is_perturb_gt_gene_min[%u]?: %s\n", i, (is_perturb_gt_gene_min[i] == true)?"yes":"no");
*/
}
// Reduce all three previous
// partial evaluations (gNorm, nIter, perturb) into a final one
is_gradDescentEn[0] = *is_gNorm_gt_gMin && *is_nIter_lt_maxIter && *is_perturb_gt_genotype;
*is_gradDescentEn = is_gNorm_gt_gMin && is_nIter_lt_maxIter && is_perturb_gt_genotype;
/*
if (get_local_id(0) == 0) {
printf("is_gNorm_gt_gMin?: %s\n", (is_gNorm_gt_gMin == true)?"yes":"no");
printf("is_nIter_lt_maxIter?: %s\n", (is_nIter_lt_maxIter == true)?"yes":"no");
printf("is_perturb_gt_genotype?: %s\n", (is_perturb_gt_genotype == true)?"yes":"no");
printf("Continue gradient iteration?: %s\n", (*is_gradDescentEn == true)?"yes":"no");
}
*/
}
barrier(CLK_LOCAL_MEM_FENCE);
return is_gradDescentEn[0];
}
void stepGPU (// Args for minimization
__local float* local_genotype, // originally as "d_x"
__local float* local_genotype_new, // originally as "d_xnew"
__local float* local_genotype_diff, // originally as "d_xdiff"
__local float* local_gradient, // originally as "d_g"
float gradMin_alpha, // originally as "alpha"
float gradMin_h, // originally as "h"
uint gradMin_inputSize, // originally as "M". initially labelled as "gradMin_M"
// Args for energy and gradient calculation
int dockpars_rotbondlist_length,
char dockpars_num_of_atoms,
char dockpars_gridsize_x,
char dockpars_gridsize_y,
char dockpars_gridsize_z,
__global const float* restrict dockpars_fgrids, // This is too large to be allocated in __constant
char dockpars_num_of_atypes,
int dockpars_num_of_intraE_contributors,
float dockpars_grid_spacing,
float dockpars_coeff_elec,
float dockpars_qasp,
float dockpars_coeff_desolv,
__local float* genotype,
__local float* energy,
__local int* run_id,
// Some OpenCL compilers don't allow declaring
// local variables within non-kernel functions.
// These local variables must be declared in a kernel,
// and then passed to non-kernel functions.
__local float* calc_coords_x,
__local float* calc_coords_y,
__local float* calc_coords_z,
__local float* partial_energies,
__constant float* atom_charges_const,
__constant char* atom_types_const,
__constant char* intraE_contributors_const,
__constant float* VWpars_AC_const,
__constant float* VWpars_BD_const,
__constant float* dspars_S_const,
__constant float* dspars_V_const,
__constant int* rotlist_const,
__constant float* ref_coords_x_const,
__constant float* ref_coords_y_const,
__constant float* ref_coords_z_const,
__constant float* rotbonds_moving_vectors_const,
__constant float* rotbonds_unit_vectors_const,
__constant float* ref_orientation_quats_const
// Gradient-related arguments
// Calculate gradients (forces) for intermolecular energy
// Derived from autodockdev/maps.py
// "is_enabled_gradient_calc": enables gradient calculation.
// In Genetic-Generation: no need for gradients
// In Gradient-Minimizer: must calculate gradients
,
__local bool* is_enabled_gradient_calc,
__local float* gradient_inter_x,
__local float* gradient_inter_y,
__local float* gradient_inter_z,
__local float* gradient_genotype
)
{
// Calculate gradient
// =============================================================
gpu_calc_energy(dockpars_rotbondlist_length,
dockpars_num_of_atoms,
dockpars_gridsize_x,
dockpars_gridsize_y,
dockpars_gridsize_z,
dockpars_fgrids,
dockpars_num_of_atypes,
dockpars_num_of_intraE_contributors,
dockpars_grid_spacing,
dockpars_coeff_elec,
dockpars_qasp,
dockpars_coeff_desolv,
genotype,
energy,
run_id,
// Some OpenCL compilers don't allow declaring
// local variables within non-kernel functions.
// These local variables must be declared in a kernel,
// and then passed to non-kernel functions.
calc_coords_x,
calc_coords_y,
calc_coords_z,
partial_energies,
atom_charges_const,
atom_types_const,
intraE_contributors_const,
VWpars_AC_const,
VWpars_BD_const,
dspars_S_const,
dspars_V_const,
rotlist_const,
ref_coords_x_const,
ref_coords_y_const,
ref_coords_z_const,
rotbonds_moving_vectors_const,
rotbonds_unit_vectors_const,
ref_orientation_quats_const
// Gradient-related arguments
// Calculate gradients (forces) for intermolecular energy
// Derived from autodockdev/maps.py
,
is_enabled_gradient_calc,
gradient_inter_x,
gradient_inter_y,
gradient_inter_z,
gradient_genotype
);
// =============================================================
for(uint i = get_local_id(0);
i < gradMin_inputSize;
i+= NUM_OF_THREADS_PER_BLOCK) {
// Taking step
local_genotype_new[i] = local_genotype[i] - gradMin_alpha * local_gradient[i];
// Updating termination metrics
local_genotype_diff[i] = local_genotype_new[i] - local_genotype[i];
// Updating current solution
local_genotype[i] = local_genotype_new[i];
}
return *is_gradDescentEn;
}
float inner_product(__local float* vector1,
......
......@@ -66,20 +66,7 @@ void gpu_calc_energy(
__constant float* ref_coords_z_const,
__constant float* rotbonds_moving_vectors_const,
__constant float* rotbonds_unit_vectors_const,
__constant float* ref_orientation_quats_const,
// Gradient-related arguments
// Calculate gradients (forces) for intermolecular energy
// Derived from autodockdev/maps.py
// "is_enabled_gradient_calc": enables gradient calculation.
// In Genetic-Generation: no need for gradients
// In Gradient-Minimizer: must calculate gradients
__local bool* is_enabled_gradient_calc,
__local float* gradient_inter_x,
__local float* gradient_inter_y,
__local float* gradient_inter_z,
__local float* gradient_genotype
__constant float* ref_orientation_quats_const
)
//The GPU device function calculates the energy of the entity described by genotype, dockpars and the liganddata
......@@ -89,20 +76,6 @@ void gpu_calc_energy(
{
partial_energies[get_local_id(0)] = 0.0f;
// -------------------------------------------------------------------
// Calculate gradients (forces) for intermolecular energy
// Derived from autodockdev/maps.py
// -------------------------------------------------------------------
if (*is_enabled_gradient_calc) {
for (uint atom_id = get_local_id(0);
atom_id < dockpars_num_of_atoms;
atom_id+= NUM_OF_THREADS_PER_BLOCK) {
gradient_inter_x[atom_id] = 0.0f;
gradient_inter_y[atom_id] = 0.0f;
gradient_inter_z[atom_id] = 0.0f;
}
}
#if 0
// Rotational genes in the Shoemake space expressed in radians
float u1, u2, u3;
......@@ -180,11 +153,9 @@ void gpu_calc_energy(
// NATIVE_PRECISION, HALF_PRECISION, Full precision
// Rotational genes in the Shoemake space expressed in radians
float u1, u2, u3;
u1 = genotype[3];
u2 = genotype[4]/**DEG_TO_RAD*/;
u3 = genotype[5]/**DEG_TO_RAD*/;
float u1 = genotype[3];
float u2 = genotype[4]/**DEG_TO_RAD*/;
float u3 = genotype[5]/**DEG_TO_RAD*/;
// u1, u2, u3 should be within their valid range of [0,1]
quatrot_left_q = native_sqrt(1 - u1) * native_sin(PI_TIMES_2*u2);
......@@ -298,57 +269,6 @@ void gpu_calc_energy(
} // End rotation_counter for-loop
// -------------------------------------------------------------------
// Calculate gradients (forces) for intermolecular energy
// Derived from autodockdev/maps.py
// -------------------------------------------------------------------
// Variables to store gradient of
// the intermolecular energy per each ligand atom
// Some OpenCL compilers don't allow declaring
// local variables within non-kernel functions.
// These local variables must be declared in a kernel,
// and then passed to non-kernel functions.
/*
__local float gradient_inter_x[MAX_NUM_OF_ATOMS];
__local float gradient_inter_y[MAX_NUM_OF_ATOMS];
__local float gradient_inter_z[MAX_NUM_OF_ATOMS];
*/
// Deltas dx, dy, dz are already normalized
// (by host/src/getparameters.cpp) in OCLaDock.
// The correspondance between vertices in xyz axes is:
// 0, 1, 2, 3, 4, 5, 6, 7 and 000, 100, 010, 001, 101, 110, 011, 111
/*
deltas: (x-x0)/(x1-x0), (y-y0...
vertices: (000, 100, 010, 001, 101, 110, 011, 111)
Z
'
3 - - - - 6
/. /|
4 - - - - 7 |
| ' | |
| 0 - - - + 2 -- Y
'/ |/
1 - - - - 5
/
X
*/
// Intermediate values for vectors in x-direction
float x10, x52, x43, x76;
float vx_z0, vx_z1;
// Intermediate values for vectors in y-direction
float y20, y51, y63, y74;
float vy_z0, vy_z1;
// Intermediate values for vectors in z-direction
float z30, z41, z62, z75;
float vz_y0, vz_y1;
// -------------------------------------------------------------------
// ================================================
// CALCULATE INTERMOLECULAR ENERGY
// ================================================
......@@ -366,19 +286,6 @@ void gpu_calc_energy(
|| (y >= dockpars_gridsize_y-1)
|| (z >= dockpars_gridsize_z-1)){
partial_energies[get_local_id(0)] += 16777216.0f; //100000.0f;
// -------------------------------------------------------------------
// Calculate gradients (forces) for intermolecular energy
// Derived from autodockdev/maps.py
// -------------------------------------------------------------------
if (*is_enabled_gradient_calc) {
// Penalty values are valid as long as they are high
gradient_inter_x[atom_id] += 16777216.0f;
gradient_inter_y[atom_id] += 16777216.0f;
gradient_inter_z[atom_id] += 16777216.0f;
}
// -------------------------------------------------------------------
}
else
{
......@@ -432,75 +339,6 @@ void gpu_calc_energy(
cube [1][0][1] = *(dockpars_fgrids + offset_cube_101 + mul_tmp);
cube [0][1][1] = *(dockpars_fgrids + offset_cube_011 + mul_tmp);
cube [1][1][1] = *(dockpars_fgrids + offset_cube_111 + mul_tmp);
// -------------------------------------------------------------------
// Calculate gradients (forces) corresponding to
// "atype" intermolecular energy
// Derived from autodockdev/maps.py
// -------------------------------------------------------------------
if (*is_enabled_gradient_calc) {
// vector in x-direction
/*
x10 = grid[int(vertices[1])] - grid[int(vertices[0])] # z = 0
x52 = grid[int(vertices[5])] - grid[int(vertices[2])] # z = 0
x43 = grid[int(vertices[4])] - grid[int(vertices[3])] # z = 1
x76 = grid[int(vertices[7])] - grid[int(vertices[6])] # z = 1
vx_z0 = (1-yd) * x10 + yd * x52 # z = 0
vx_z1 = (1-yd) * x43 + yd * x76 # z = 1
gradient[0] = (1-zd) * vx_z0 + zd * vx_z1
*/
x10 = cube [1][0][0] - cube [0][0][0]; // z = 0
x52 = cube [1][1][0] - cube [0][1][0]; // z = 0
x43 = cube [1][0][1] - cube [0][0][1]; // z = 1
x76 = cube [1][1][1] - cube [0][1][1]; // z = 1
vx_z0 = (1 - dy) * x10 + dy * x52; // z = 0
vx_z1 = (1 - dy) * x43 + dy * x76; // z = 1
gradient_inter_x[atom_id] += (1 - dz) * vx_z0 + dz * vx_z1;
// vector in y-direction
/*
y20 = grid[int(vertices[2])] - grid[int(vertices[0])] # z = 0
y51 = grid[int(vertices[5])] - grid[int(vertices[1])] # z = 0
y63 = grid[int(vertices[6])] - grid[int(vertices[3])] # z = 1
for (uint rotation_counter = get_local_id(0);
rotation_counter < dockpars_rotbondlist_length;
rotation_counter+=NUM_OF_THREADS_PER_BLOCK)
{ y74 = grid[int(vertices[7])] - grid[int(vertices[4])] # z = 1
vy_z0 = (1-xd) * y20 + xd * y51 # z = 0
vy_z1 = (1-xd) * y63 + xd * y74 # z = 1
gradient[1] = (1-zd) * vy_z0 + zd * vy_z1
*/
y20 = cube[0][1][0] - cube [0][0][0]; // z = 0
y51 = cube[1][1][0] - cube [1][0][0]; // z = 0
y63 = cube[0][1][1] - cube [0][0][1]; // z = 1
y74 = cube[1][1][1] - cube [1][0][1]; // z = 1
vy_z0 = (1 - dx) * y20 + dx * y51; // z = 0
vy_z1 = (1 - dx) * y63 + dx * y74; // z = 1
gradient_inter_y[atom_id] += (1 - dz) * vy_z0 + dz * vy_z1;
// vectors in z-direction
/*
z30 = grid[int(vertices[3])] - grid[int(vertices[0])] # y = 0
z41 = grid[int(vertices[4])] - grid[int(vertices[1])] # y = 0
z62 = grid[int(vertices[6])] - grid[int(vertices[2])] # y = 1
z75 = grid[int(vertices[7])] - grid[int(vertices[5])] # y = 1
vz_y0 = (1-xd) * z30 + xd * z41 # y = 0
vz_y1 = (1-xd) * z62 + xd * z75 # y = 1
gradient[2] = (1-yd) * vz_y0 + yd * vz_y1
*/
z30 = cube [0][0][1] - cube [0][0][0]; // y = 0
z41 = cube [1][0][1] - cube [1][0][0]; // y = 0
z62 = cube [0][1][1] - cube [0][1][0]; // y = 1
z75 = cube [1][1][1] - cube [1][1][0]; // y = 1
vz_y0 = (1 - dx) * z30 + dx * z41; // y = 0
vz_y1 = (1 - dx) * z62 + dx * z75; // y = 1
gradient_inter_z[atom_id] += (1 - dy) * vz_y0 + dy * vz_y1;
}
// -------------------------------------------------------------------
#else
// -------------------------------------------------------------------
// FIXME: this block within the "#else" preprocessor directive
......@@ -565,42 +403,6 @@ void gpu_calc_energy(
cube [1][0][1] = *(dockpars_fgrids + offset_cube_101 + mul_tmp);
cube [0][1][1] = *(dockpars_fgrids + offset_cube_011 + mul_tmp);
cube [1][1][1] = *(dockpars_fgrids + offset_cube_111 + mul_tmp);
// -------------------------------------------------------------------
// Calculate gradients (forces) corresponding to
// "elec" intermolecular energy
// Derived from autodockdev/maps.py
// -------------------------------------------------------------------
if (*is_enabled_gradient_calc) {
// vector in x-direction
x10 = cube [1][0][0] - cube [0][0][0]; // z = 0
x52 = cube [1][1][0] - cube [0][1][0]; // z = 0
x43 = cube [1][0][1] - cube [0][0][1]; // z = 1
x76 = cube [1][1][1] - cube [0][1][1]; // z = 1
vx_z0 = (1 - dy) * x10 + dy * x52; // z = 0
vx_z1 = (1 - dy) * x43 + dy * x76; // z = 1
gradient_inter_x[atom_id] += (1 - dz) * vx_z0 + dz * vx_z1;
// vector in y-direction
y20 = cube[0][1][0] - cube [0][0][0]; // z = 0
y51 = cube[1][1][0] - cube [1][0][0]; // z = 0
y63 = cube[0][1][1] - cube [0][0][1]; // z = 1
y74 = cube[1][1][1] - cube [1][0][1]; // z = 1
vy_z0 = (1 - dx) * y20 + dx * y51; // z = 0
vy_z1 = (1 - dx) * y63 + dx * y74; // z = 1
gradient_inter_y[atom_id] += (1 - dz) * vy_z0 + dz * vy_z1;
// vectors in z-direction
z30 = cube [0][0][1] - cube [0][0][0]; // y = 0
z41 = cube [1][0][1] - cube [1][0][0]; // y = 0
z62 = cube [0][1][1] - cube [0][1][0]; // y = 1
z75 = cube [1][1][1] - cube [1][1][0]; // y = 1
vz_y0 = (1 - dx) * z30 + dx * z41; // y = 0
vz_y1 = (1 - dx) * z62 + dx * z75; // y = 1
gradient_inter_z[atom_id] += (1 - dy) * vz_y0 + dy * vz_y1;
}
// -------------------------------------------------------------------
#else
// -------------------------------------------------------------------
// FIXME: this block within the "#else" preprocessor directive
......@@ -665,42 +467,6 @@ void gpu_calc_energy(
cube [1][0][1] = *(dockpars_fgrids + offset_cube_101 + mul_tmp);
cube [0][1][1] = *(dockpars_fgrids + offset_cube_011 + mul_tmp);
cube [1][1][1] = *(dockpars_fgrids + offset_cube_111 + mul_tmp);
// -------------------------------------------------------------------
// Calculate gradients (forces) corresponding to
// "dsol" intermolecular energy
// Derived from autodockdev/maps.py
// -------------------------------------------------------------------
if (*is_enabled_gradient_calc) {
// vector in x-direction
x10 = cube [1][0][0] - cube [0][0][0]; // z = 0
x52 = cube [1][1][0] - cube [0][1][0]; // z = 0
x43 = cube [1][0][1] - cube [0][0][1]; // z = 1
x76 = cube [1][1][1] - cube [0][1][1]; // z = 1
vx_z0 = (1 - dy) * x10 + dy * x52; // z = 0
vx_z1 = (1 - dy) * x43 + dy * x76; // z = 1
gradient_inter_x[atom_id] += (1 - dz) * vx_z0 + dz * vx_z1;
// vector in y-direction
y20 = cube[0][1][0] - cube [0][0][0]; // z = 0
y51 = cube[1][1][0] - cube [1][0][0]; // z = 0
y63 = cube[0][1][1] - cube [0][0][1]; // z = 1
y74 = cube[1][1][1] - cube [1][0][1]; // z = 1
vy_z0 = (1 - dx) * y20 + dx * y51; // z = 0
vy_z1 = (1 - dx) * y63 + dx * y74; // z = 1
gradient_inter_y[atom_id] += (1 - dz) * vy_z0 + dz * vy_z1;
// vectors in z-direction
z30 = cube [0][0][1] - cube [0][0][0]; // y = 0
z41 = cube [1][0][1] - cube [1][0][0]; // y = 0
z62 = cube [0][1][1] - cube [0][1][0]; // y = 1
z75 = cube [1][1][1] - cube [1][1][0]; // y = 1
vz_y0 = (1 - dx) * z30 + dx * z41; // y = 0
vz_y1 = (1 - dx) * z62 + dx * z75; // y = 1
gradient_inter_z[atom_id] += (1 - dy) * vz_y0 + dy * vz_y1;
}
// -------------------------------------------------------------------
#else
// -------------------------------------------------------------------
// FIXME: this block within the "#else" preprocessor directive
......@@ -870,14 +636,6 @@ void gpu_calc_energy(
}
} // End contributor_counter for-loop (INTRAMOLECULAR ENERGY)
barrier(CLK_LOCAL_MEM_FENCE);
if (get_local_id(0) == 0)
......@@ -892,232 +650,6 @@ void gpu_calc_energy(
}
}
// -------------------------------------------------------------------
// Calculate gradients (forces) corresponding to (interE + intraE)
// Derived from autodockdev/motions.py/forces_to_delta()
// -------------------------------------------------------------------
// Could be barrier removed if another work-item is used?
// (e.g. get_locla_id(0) == 1)
barrier(CLK_LOCAL_MEM_FENCE);
// FIXME: done so far only for interE
if (get_local_id(0) == 0) {
if (*is_enabled_gradient_calc) {
gradient_genotype [0] = 0.0f;
gradient_genotype [1] = 0.0f;
gradient_genotype [2] = 0.0f;
// ------------------------------------------
// translation-related gradients
// ------------------------------------------
for (uint lig_atom_id = 0;
lig_atom_id<dockpars_num_of_atoms;
lig_atom_id++) {
gradient_genotype [0] += gradient_inter_x[lig_atom_id]; // gradient for gene 0: gene x
gradient_genotype [1] += gradient_inter_y[lig_atom_id]; // gradient for gene 1: gene y
gradient_genotype [2] += gradient_inter_z[lig_atom_id]; // gradient for gene 2: gene z
}
// ------------------------------------------
// rotation-related gradients
// Transform gradients_inter_{x|y|z}
// into local_gradients[i] (with four quaternion genes)
// Derived from autodockdev/motions.py/forces_to_delta_genes()
// Transform local_gradients[i] (with four quaternion genes)
// into local_gradients[i] (with three Shoemake genes)
// Derived from autodockdev/motions.py/_get_cube3_gradient()
// ------------------------------------------
float3 torque_rot = (float3)(0.0f, 0.0f, 0.0f);
// center of rotation
// In getparameters.cpp, it indicates
// translation genes are in grid spacing (instead of Angstroms)
float about[3];
about[0] = genotype[0];
about[1] = genotype[1];
about[2] = genotype[2];
// Temporal variable to calculate translation differences.
// They are converted back to Angstroms here
float3 r;
for (uint lig_atom_id = 0;
lig_atom_id<dockpars_num_of_atoms;
lig_atom_id++) {
r.x = (calc_coords_x[lig_atom_id] - about[0]) * dockpars_grid_spacing;
r.y = (calc_coords_y[lig_atom_id] - about[1]) * dockpars_grid_spacing;
r.z = (calc_coords_z[lig_atom_id] - about[2]) * dockpars_grid_spacing;
torque_rot += cross(r, torque_rot);
}
const float rad = 1E-8;
const float rad_div_2 = native_divide(rad, 2);
float quat_w, quat_x, quat_y, quat_z;
// Derived from rotation.py/axisangle_to_q()
// genes[3:7] = rotation.axisangle_to_q(torque, rad)
torque_rot = fast_normalize(torque_rot);
quat_x = torque_rot.x;
quat_y = torque_rot.y;
quat_z = torque_rot.z;
// rotation-related gradients are expressed here in quaternions
quat_w = native_cos(rad_div_2);
quat_x = quat_x * native_sin(rad_div_2);
quat_y = quat_y * native_sin(rad_div_2);
quat_z = quat_z * native_sin(rad_div_2);
// convert quaternion gradients into Shoemake gradients
// Derived from autodockdev/motion.py/_get_cube3_gradient