Commit 17b12c3f authored by lvs's avatar lvs
Browse files

corredted fire

parent 20cf44f1
......@@ -207,10 +207,10 @@ NRUN := 200
POPSIZE := 500
TESTNAME := test
TESTLS := fire
NUM_LSIT := 300
test: odock
$(BIN_DIR)/$(TARGET) -ffile ./input/$(PDB)/derived/$(PDB)_protein.maps.fld -lfile ./input/$(PDB)/derived/$(PDB)_ligand.pdbqt -nrun $(NRUN) -psize $(POPSIZE) -resnam $(TESTNAME) -gfpop 1 -lsmet $(TESTLS) -lsit $(NUM_LSIT) -smooth 0.5
$(BIN_DIR)/$(TARGET) -ffile ./input/$(PDB)/derived/$(PDB)_protein.maps.fld -lfile ./input/$(PDB)/derived/$(PDB)_ligand.pdbqt -nrun $(NRUN) -psize $(POPSIZE) -resnam $(TESTNAME) -gfpop 0 -lsmet $(TESTLS)
ASTEX_PDB := 2bsm
ASTEX_NRUN:= 10
......
/*
OCLADock, an OpenCL implementation of AutoDock 4.2 running a Lamarckian Genetic Algorithm
Copyright (C) 2017 TU Darmstadt, Embedded Systems and Applications Group, Germany. All rights reserved.
AutoDock is a Trade Mark of the Scripps Research Institute.
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
//#define DEBUG_ENERGY_KERNEL
/*
#include "calcenergy_basic.h"
*/
// All related pragmas are in defines.h (accesible by host and device code)
// The GPU device function calculates the energy's gradient (forces or derivatives)
// of the entity described by genotype, dockpars and the ligand-data
// arrays in constant memory and returns it in the "gradient_genotype" parameter.
// The parameter "run_id" has to be equal to the ID of the run
// whose population includes the current entity (which can be determined with get_group_id(0)),
// since this determines which reference orientation should be used.
//#define DEBUG_GRAD_TRANSLATION_GENES
//#define DEBUG_GRAD_ROTATION_GENES
//#define DEBUG_GRAD_TORSION_GENES
// Atomic operations used in gradients of intra contributors.
// Only atomic_cmpxchg() works on floats.
// So for atomic add on floats, this link was used:
// https://streamhpc.com/blog/2016-02-09/atomic-operations-for-floats-in-opencl-improved/
void atomicAdd_g_f(volatile __local float *addr, float val)
{
union{
unsigned int u32;
float f32;
} next, expected, current;
current.f32 = *addr;
do{
expected.f32 = current.f32;
next.f32 = expected.f32 + val;
current.u32 = atomic_cmpxchg( (volatile __local unsigned int *)addr, expected.u32, next.u32);
} while( current.u32 != expected.u32 );
}
void atomicSub_g_f(volatile __local float *addr, float val)
{
union{
unsigned int u32;
float f32;
} next, expected, current;
current.f32 = *addr;
do{
expected.f32 = current.f32;
next.f32 = expected.f32 - val;
current.u32 = atomic_cmpxchg( (volatile __local unsigned int *)addr, expected.u32, next.u32);
} while( current.u32 != expected.u32 );
}
void gpu_calc_EG(
int dockpars_rotbondlist_length,
char dockpars_num_of_atoms,
char dockpars_gridsize_x,
char dockpars_gridsize_y,
char dockpars_gridsize_z,
// g1 = gridsize_x
uint dockpars_gridsize_x_times_y, // g2 = gridsize_x * gridsize_y
uint dockpars_gridsize_x_times_y_times_z, // g3 = gridsize_x * gridsize_y * 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,
// 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* genotype,
__local float* energy,
__local int* run_id,
__local float* calc_coords_x,
__local float* calc_coords_y,
__local float* calc_coords_z,
__local float* partial_energies,
#if defined (DEBUG_ENERGY_KERNEL)
__local float* partial_interE,
__local float* partial_intraE,
#endif
__constant float* atom_charges_const,
__constant char* atom_types_const,
__constant char* intraE_contributors_const,
#if 0
bool debug,
#endif
float dockpars_smooth,
__constant float* reqm,
__constant float* reqm_hbond,
__constant uint* atom1_types_reqm,
__constant uint* atom2_types_reqm,
__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
__constant int* rotbonds_const,
__constant int* rotbonds_atoms_const,
__constant int* num_rotating_atoms_per_rotbond_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
,
int dockpars_num_of_genes,
__local float* gradient_inter_x,
__local float* gradient_inter_y,
__local float* gradient_inter_z,
__local float* gradient_intra_x,
__local float* gradient_intra_y,
__local float* gradient_intra_z,
__local float* gradient_genotype
)
//The GPU device function calculates the energy of the entity described by genotype, dockpars and the liganddata
//arrays in constant memory and returns it in the energy parameter. The parameter run_id has to be equal to the ID
//of the run whose population includes the current entity (which can be determined with blockIdx.x), since this
//determines which reference orientation should be used.
{
partial_energies[get_local_id(0)] = 0.0f;
#if defined (DEBUG_ENERGY_KERNEL)
partial_interE[get_local_id(0)] = 0.0f;
partial_intraE[get_local_id(0)] = 0.0f;
#endif
// Initializing gradients (forces)
// Derived from autodockdev/maps.py
for (uint atom_id = get_local_id(0);
atom_id < dockpars_num_of_atoms;
atom_id+= NUM_OF_THREADS_PER_BLOCK) {
// Intermolecular gradients
gradient_inter_x[atom_id] = 0.0f;
gradient_inter_y[atom_id] = 0.0f;
gradient_inter_z[atom_id] = 0.0f;
// Intramolecular gradients
gradient_intra_x[atom_id] = 0.0f;
gradient_intra_y[atom_id] = 0.0f;
gradient_intra_z[atom_id] = 0.0f;
}
// Initializing gradient genotypes
for (uint gene_cnt = get_local_id(0);
gene_cnt < dockpars_num_of_genes;
gene_cnt+= NUM_OF_THREADS_PER_BLOCK) {
gradient_genotype[gene_cnt] = 0.0f;
}
barrier(CLK_LOCAL_MEM_FENCE);
uchar g1 = dockpars_gridsize_x;
uint g2 = dockpars_gridsize_x_times_y;
uint g3 = dockpars_gridsize_x_times_y_times_z;
// ================================================
// CALCULATING ATOMIC POSITIONS AFTER ROTATIONS
// ================================================
for (uint rotation_counter = get_local_id(0);
rotation_counter < dockpars_rotbondlist_length;
rotation_counter+=NUM_OF_THREADS_PER_BLOCK)
{
int rotation_list_element = rotlist_const[rotation_counter];
if ((rotation_list_element & RLIST_DUMMY_MASK) == 0) // If not dummy rotation
{
uint atom_id = rotation_list_element & RLIST_ATOMID_MASK;
// Capturing atom coordinates
float atom_to_rotate[3];
if ((rotation_list_element & RLIST_FIRSTROT_MASK) != 0) // If first rotation of this atom
{
atom_to_rotate[0] = ref_coords_x_const[atom_id];
atom_to_rotate[1] = ref_coords_y_const[atom_id];
atom_to_rotate[2] = ref_coords_z_const[atom_id];
}
else
{
atom_to_rotate[0] = calc_coords_x[atom_id];
atom_to_rotate[1] = calc_coords_y[atom_id];
atom_to_rotate[2] = calc_coords_z[atom_id];
}
// Capturing rotation vectors and angle
float rotation_movingvec[3];
float quatrot_left_x, quatrot_left_y, quatrot_left_z, quatrot_left_q;
float quatrot_temp_x, quatrot_temp_y, quatrot_temp_z, quatrot_temp_q;
if ((rotation_list_element & RLIST_GENROT_MASK) != 0) // If general rotation
{
// Rotational genes in the Shoemake space expressed in radians
float u1 = genotype[3];
float u2 = genotype[4];
float u3 = genotype[5];
// 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);
quatrot_left_x = native_sqrt(1 - u1) * native_cos(PI_TIMES_2*u2);
quatrot_left_y = native_sqrt(u1) * native_sin(PI_TIMES_2*u3);
quatrot_left_z = native_sqrt(u1) * native_cos(PI_TIMES_2*u3);
rotation_movingvec[0] = genotype[0];
rotation_movingvec[1] = genotype[1];
rotation_movingvec[2] = genotype[2];
}
else // If rotating around rotatable bond
{
uint rotbond_id = (rotation_list_element & RLIST_RBONDID_MASK) >> RLIST_RBONDID_SHIFT;
float rotation_unitvec[3];
rotation_unitvec[0] = rotbonds_unit_vectors_const[3*rotbond_id];
rotation_unitvec[1] = rotbonds_unit_vectors_const[3*rotbond_id+1];
rotation_unitvec[2] = rotbonds_unit_vectors_const[3*rotbond_id+2];
float rotation_angle = genotype[6+rotbond_id]*DEG_TO_RAD;
rotation_movingvec[0] = rotbonds_moving_vectors_const[3*rotbond_id];
rotation_movingvec[1] = rotbonds_moving_vectors_const[3*rotbond_id+1];
rotation_movingvec[2] = rotbonds_moving_vectors_const[3*rotbond_id+2];
// Performing additionally the first movement which
// is needed only if rotating around rotatable bond
atom_to_rotate[0] -= rotation_movingvec[0];
atom_to_rotate[1] -= rotation_movingvec[1];
atom_to_rotate[2] -= rotation_movingvec[2];
// Transforming torsion angles into quaternions
rotation_angle = rotation_angle * 0.5f;
float sin_angle = native_sin(rotation_angle);
quatrot_left_q = native_cos(rotation_angle);
quatrot_left_x = sin_angle*rotation_unitvec[0];
quatrot_left_y = sin_angle*rotation_unitvec[1];
quatrot_left_z = sin_angle*rotation_unitvec[2];
}
// Performing rotation
if ((rotation_list_element & RLIST_GENROT_MASK) != 0) // If general rotation,
// two rotations should be performed
// (multiplying the quaternions)
{
// Calculating quatrot_left*ref_orientation_quats_const,
// which means that reference orientation rotation is the first
quatrot_temp_q = quatrot_left_q;
quatrot_temp_x = quatrot_left_x;
quatrot_temp_y = quatrot_left_y;
quatrot_temp_z = quatrot_left_z;
quatrot_left_q = quatrot_temp_q*ref_orientation_quats_const[4*(*run_id)]-
quatrot_temp_x*ref_orientation_quats_const[4*(*run_id)+1]-
quatrot_temp_y*ref_orientation_quats_const[4*(*run_id)+2]-
quatrot_temp_z*ref_orientation_quats_const[4*(*run_id)+3];
quatrot_left_x = quatrot_temp_q*ref_orientation_quats_const[4*(*run_id)+1]+
ref_orientation_quats_const[4*(*run_id)]*quatrot_temp_x+
quatrot_temp_y*ref_orientation_quats_const[4*(*run_id)+3]-
ref_orientation_quats_const[4*(*run_id)+2]*quatrot_temp_z;
quatrot_left_y = quatrot_temp_q*ref_orientation_quats_const[4*(*run_id)+2]+
ref_orientation_quats_const[4*(*run_id)]*quatrot_temp_y+
ref_orientation_quats_const[4*(*run_id)+1]*quatrot_temp_z-
quatrot_temp_x*ref_orientation_quats_const[4*(*run_id)+3];
quatrot_left_z = quatrot_temp_q*ref_orientation_quats_const[4*(*run_id)+3]+
ref_orientation_quats_const[4*(*run_id)]*quatrot_temp_z+
quatrot_temp_x*ref_orientation_quats_const[4*(*run_id)+2]-
ref_orientation_quats_const[4*(*run_id)+1]*quatrot_temp_y;
}
quatrot_temp_q = 0 -
quatrot_left_x*atom_to_rotate [0] -
quatrot_left_y*atom_to_rotate [1] -
quatrot_left_z*atom_to_rotate [2];
quatrot_temp_x = quatrot_left_q*atom_to_rotate [0] +
quatrot_left_y*atom_to_rotate [2] -
quatrot_left_z*atom_to_rotate [1];
quatrot_temp_y = quatrot_left_q*atom_to_rotate [1] -
quatrot_left_x*atom_to_rotate [2] +
quatrot_left_z*atom_to_rotate [0];
quatrot_temp_z = quatrot_left_q*atom_to_rotate [2] +
quatrot_left_x*atom_to_rotate [1] -
quatrot_left_y*atom_to_rotate [0];
atom_to_rotate [0] = 0 -
quatrot_temp_q*quatrot_left_x +
quatrot_temp_x*quatrot_left_q -
quatrot_temp_y*quatrot_left_z +
quatrot_temp_z*quatrot_left_y;
atom_to_rotate [1] = 0 -
quatrot_temp_q*quatrot_left_y +
quatrot_temp_x*quatrot_left_z +
quatrot_temp_y*quatrot_left_q -
quatrot_temp_z*quatrot_left_x;
atom_to_rotate [2] = 0 -
quatrot_temp_q*quatrot_left_z -
quatrot_temp_x*quatrot_left_y +
quatrot_temp_y*quatrot_left_x +
quatrot_temp_z*quatrot_left_q;
// Performing final movement and storing values
calc_coords_x[atom_id] = atom_to_rotate [0] + rotation_movingvec[0];
calc_coords_y[atom_id] = atom_to_rotate [1] + rotation_movingvec[1];
calc_coords_z[atom_id] = atom_to_rotate [2] + rotation_movingvec[2];
} // End if-statement not dummy rotation
barrier(CLK_LOCAL_MEM_FENCE);
} // End rotation_counter for-loop
// ================================================
// CALCULATING INTERMOLECULAR ENERGY
// ================================================
for (uint atom_id = get_local_id(0);
atom_id < dockpars_num_of_atoms;
atom_id+= NUM_OF_THREADS_PER_BLOCK)
{
uint atom_typeid = atom_types_const[atom_id];
float x = calc_coords_x[atom_id];
float y = calc_coords_y[atom_id];
float z = calc_coords_z[atom_id];
float q = atom_charges_const[atom_id];
if ((x < 0) || (y < 0) || (z < 0) || (x >= dockpars_gridsize_x-1)
|| (y >= dockpars_gridsize_y-1)
|| (z >= dockpars_gridsize_z-1)){
partial_energies[get_local_id(0)] += 16777216.0f; //100000.0f;
// Setting gradients (forces) penalties.
// These 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;
#if defined (DEBUG_ENERGY_KERNEL)
partial_interE[get_local_id(0)] += 16777216.0f;
#endif
}
else
{
// Getting coordinates
int x_low = (int)floor(x);
int y_low = (int)floor(y);
int z_low = (int)floor(z);
int x_high = (int)ceil(x);
int y_high = (int)ceil(y);
int z_high = (int)ceil(z);
float dx = x - x_low;
float dy = y - y_low;
float dz = z - z_low;
// Calculating interpolation weights
float weights[2][2][2];
weights [0][0][0] = (1-dx)*(1-dy)*(1-dz);
weights [1][0][0] = dx*(1-dy)*(1-dz);
weights [0][1][0] = (1-dx)*dy*(1-dz);
weights [1][1][0] = dx*dy*(1-dz);
weights [0][0][1] = (1-dx)*(1-dy)*dz;
weights [1][0][1] = dx*(1-dy)*dz;
weights [0][1][1] = (1-dx)*dy*dz;
weights [1][1][1] = dx*dy*dz;
// Capturing affinity values
uint ylow_times_g1 = y_low*g1;
uint yhigh_times_g1 = y_high*g1;
uint zlow_times_g2 = z_low*g2;
uint zhigh_times_g2 = z_high*g2;
// Grid offset
uint offset_cube_000 = x_low + ylow_times_g1 + zlow_times_g2;
uint offset_cube_100 = x_high + ylow_times_g1 + zlow_times_g2;
uint offset_cube_010 = x_low + yhigh_times_g1 + zlow_times_g2;
uint offset_cube_110 = x_high + yhigh_times_g1 + zlow_times_g2;
uint offset_cube_001 = x_low + ylow_times_g1 + zhigh_times_g2;
uint offset_cube_101 = x_high + ylow_times_g1 + zhigh_times_g2;
uint offset_cube_011 = x_low + yhigh_times_g1 + zhigh_times_g2;
uint offset_cube_111 = x_high + yhigh_times_g1 + zhigh_times_g2;
uint mul_tmp = atom_typeid*g3;
float cube[2][2][2];
cube [0][0][0] = *(dockpars_fgrids + offset_cube_000 + mul_tmp);
cube [1][0][0] = *(dockpars_fgrids + offset_cube_100 + mul_tmp);
cube [0][1][0] = *(dockpars_fgrids + offset_cube_010 + mul_tmp);
cube [1][1][0] = *(dockpars_fgrids + offset_cube_110 + mul_tmp);
cube [0][0][1] = *(dockpars_fgrids + offset_cube_001 + mul_tmp);
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);
// Calculating affinity energy
partial_energies[get_local_id(0)] += TRILININTERPOL(cube, weights);
#if defined (DEBUG_ENERGY_KERNEL)
partial_interE[get_local_id(0)] += TRILININTERPOL(cube, weights);
#endif
// -------------------------------------------------------------------
// 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;
// -------------------------------------------------------------------
// Calculating gradients (forces) corresponding to
// "atype" intermolecular energy
// Derived from autodockdev/maps.py
// -------------------------------------------------------------------
// 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;
//printf("%-15s %-3u %-10.8f %-10.8f %-10.8f %-10.8f %-10.8f %-10.8f\n", "atom aff", atom_id, vx_z0, vx_z1, vy_z0, vy_z1, vz_y0, vz_y1);
// Capturing electrostatic values
atom_typeid = dockpars_num_of_atypes;
mul_tmp = atom_typeid*g3;
cube [0][0][0] = *(dockpars_fgrids + offset_cube_000 + mul_tmp);
cube [1][0][0] = *(dockpars_fgrids + offset_cube_100 + mul_tmp);
cube [0][1][0] = *(dockpars_fgrids + offset_cube_010 + mul_tmp);
cube [1][1][0] = *(dockpars_fgrids + offset_cube_110 + mul_tmp);
cube [0][0][1] = *(dockpars_fgrids + offset_cube_001 + mul_tmp);
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);
// Calculating electrostatic energy
partial_energies[get_local_id(0)] += q * TRILININTERPOL(cube, weights);
#if defined (DEBUG_ENERGY_KERNEL)
partial_interE[get_local_id(0)] += q * TRILININTERPOL(cube, weights);
#endif
// -------------------------------------------------------------------
// Calculating gradients (forces) corresponding to
// "elec" intermolecular energy
// Derived from autodockdev/maps.py
// -------------------------------------------------------------------
/*
// Capturing electrostatic values
atom_typeid = dockpars_num_of_atypes;
mul_tmp = atom_typeid*g3;
cube [0][0][0] = *(dockpars_fgrids + offset_cube_000 + mul_tmp);
cube [1][0][0] = *(dockpars_fgrids + offset_cube_100 + mul_tmp);
cube [0][1][0] = *(dockpars_fgrids + offset_cube_010 + mul_tmp);
cube [1][1][0] = *(dockpars_fgrids + offset_cube_110 + mul_tmp);
cube [0][0][1] = *(dockpars_fgrids + offset_cube_001 + mul_tmp);
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);
*/
// 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] += q * ((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] += q *((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] += q *((1 - dy) * vz_y0 + dy * vz_y1);
//printf("%-15s %-3u %-10.8f %-10.8f %-10.8f %-10.8f %-10.8f %-10.8f\n", "elec", atom_id, vx_z0, vx_z1, vy_z0, vy_z1, vz_y0, vz_y1);
// Capturing desolvation values
atom_typeid = dockpars_num_of_atypes+1;
mul_tmp = atom_typeid*g3;
cube [0][0][0] = *(dockpars_fgrids + offset_cube_000 + mul_tmp);
cube [1][0][0] = *(dockpars_fgrids + offset_cube_100 + mul_tmp);
cube [0][1][0] = *(dockpars_fgrids + offset_cube_010 + mul_tmp);
cube [1][1][0] = *(dockpars_fgrids + offset_cube_110 + mul_tmp);
cube [0][0][1] = *(dockpars_fgrids + offset_cube_001 + mul_tmp);
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);
// Calculating desolvation energy
partial_energies[get_local_id(0)] += fabs(q) * TRILININTERPOL(cube, weights);