kernel_sd.cl 25.7 KB
Newer Older
1
2
// Gradient-based steepest descent minimizer
// Alternative to Solis-Wetts
3

lvs's avatar
lvs committed
4
5
6
7
8
//#define DEBUG_ENERGY_KERNEL5
	//#define PRINT_ENERGIES
	//#define PRINT_GENES_AND_GRADS
	//#define PRINT_ATOMIC_COORDS

lvs's avatar
lvs committed
9
10
11
// Enable DEBUG_MINIMIZER for a seeing a detailed SD evolution
// If only PRINT_MINIMIZER_ENERGY_EVOLUTION is enabled,
// then a only a simplified SD evolution will be shown
12
//#define DEBUG_MINIMIZER
lvs's avatar
lvs committed
13
//	#define PRINT_MINIMIZER_ENERGY_EVOLUTION
14

15
// Enable this for debugging SD from a defined initial genotype
lvs's avatar
lvs committed
16
//#define DEBUG_INITIAL_2BRT
17

18
__kernel void __attribute__ ((reqd_work_group_size(NUM_OF_THREADS_PER_BLOCK,1,1)))
lvs's avatar
lvs committed
19
gradient_minSD(	
lvs's avatar
lvs committed
20
21
22
23
24
25
			    char   dockpars_num_of_atoms,
			    char   dockpars_num_of_atypes,
			    int    dockpars_num_of_intraE_contributors,
			    char   dockpars_gridsize_x,
			    char   dockpars_gridsize_y,
			    char   dockpars_gridsize_z,
26
							    		// g1 = gridsize_x
lvs's avatar
lvs committed
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
			    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
			    float  dockpars_grid_spacing,
	 __global const     float* restrict dockpars_fgrids, 		// This is too large to be allocated in __constant 
			    int    dockpars_rotbondlist_length,
			    float  dockpars_coeff_elec,
			    float  dockpars_coeff_desolv,
	  __global          float* restrict dockpars_conformations_next,
	  __global          float* restrict dockpars_energies_next,
  	  __global 	    int*   restrict dockpars_evals_of_new_entities,
	  __global          uint*  restrict dockpars_prng_states,
			    int    dockpars_pop_size,
			    int    dockpars_num_of_genes,
			    float  dockpars_lsearch_rate,
			    uint   dockpars_num_of_lsentities,
			    uint   dockpars_max_num_of_iters,
			    float  dockpars_qasp,
			    float  dockpars_smooth,

	  __constant        kernelconstant_interintra* 	 kerconst_interintra,
	  __global const    kernelconstant_intracontrib* kerconst_intracontrib,
	  __constant        kernelconstant_intra*	 kerconst_intra,
	  __constant        kernelconstant_rotlist*   	 kerconst_rotlist,
	  __constant        kernelconstant_conform*	 kerconst_conform
51
			,
lvs's avatar
lvs committed
52
53
54
55
56
57
58
	  __constant int*   	  rotbonds_const,
	  __global   const int*   rotbonds_atoms_const,
	  __constant int*   	  num_rotating_atoms_per_rotbond_const
			,
	  __global   const float* angle_const,
	  __constant       float* dependence_on_theta_const,
	  __constant       float* dependence_on_rotangle_const
59
60
61
62
63
64
65
66
67
)
//The GPU global function performs gradient-based minimization on (some) entities of conformations_next.
//The number of OpenCL compute units (CU) which should be started equals to num_of_minEntities*num_of_runs.
//This way the first num_of_lsentities entity of each population will be subjected to local search
//(and each CU carries out the algorithm for one entity).
//Since the first entity is always the best one in the current population,
//it is always tested according to the ls probability, and if it not to be
//subjected to local search, the entity with ID num_of_lsentities is selected instead of the first one (with ID 0).
{
68
	// -----------------------------------------------------------------------------
69
70
71
	// -----------------------------------------------------------------------------
	// -----------------------------------------------------------------------------

72
	// Determining entity, and its run, energy, and genotype
73
74
	__local int   entity_id;
	__local int   run_id;
Leonardo Solis's avatar
Leonardo Solis committed
75
76
77
  	__local float energy;
	__local float genotype[ACTUAL_GENOTYPE_LENGTH];

lvs's avatar
lvs committed
78
79
	// Iteration counter fot the minimizer
  	__local uint iteration_cnt;  	
80

Leonardo Solis's avatar
Leonardo Solis committed
81
82
83
	// Stepsize for the minimizer
	__local float stepsize;

84
85
	if (get_local_id(0) == 0)
	{
86
		// Choosing a random entity out of the entire population
lvs's avatar
lvs committed
87
		/*
88
		run_id = get_group_id(0);
Leonardo Solis's avatar
Leonardo Solis committed
89
90
		//entity_id = (uint)(dockpars_pop_size * gpu_randf(dockpars_prng_states));
		entity_id = 0;
lvs's avatar
lvs committed
91
		*/
92

93
94
		run_id = get_group_id(0) / dockpars_num_of_lsentities;
		entity_id = get_group_id(0) % dockpars_num_of_lsentities;
95
96
97
98
99
100
101
102
103
104

		// Since entity 0 is the best one due to elitism,
		// it should be subjected to random selection
		if (entity_id == 0) {
			// If entity 0 is not selected according to LS-rate,
			// choosing an other entity
			if (100.0f*gpu_randf(dockpars_prng_states) > dockpars_lsearch_rate) {
				entity_id = dockpars_num_of_lsentities;					
			}
		}
105
		
Leonardo Solis's avatar
Leonardo Solis committed
106
		energy = dockpars_energies_next[run_id*dockpars_pop_size+entity_id];
107

Leonardo Solis's avatar
Leonardo Solis committed
108
109
		// Initializing gradient-minimizer counters and flags
    		iteration_cnt  = 0;
Leonardo Solis's avatar
Leonardo Solis committed
110
		stepsize       = STEP_START;
111
112

		#if defined (DEBUG_MINIMIZER) || defined (PRINT_MINIMIZER_ENERGY_EVOLUTION)
113
114
115
116
117
118
119
120
		printf("\n");
		printf("-------> Start of SD minimization cycle\n");
		printf("%20s %6u\n", "run_id: ", run_id);
		printf("%20s %6u\n", "entity_id: ", entity_id);
		printf("\n");
		printf("%20s \n", "LGA genotype: ");
		printf("%20s %.6f\n", "initial energy: ", energy);
		printf("%20s %.6f\n", "initial stepsize: ", stepsize);
121
		#endif
122
123
	}

124
125
	barrier(CLK_LOCAL_MEM_FENCE);

lvs's avatar
lvs committed
126
127
128
  	event_t ev = async_work_group_copy(genotype,
  			      		   dockpars_conformations_next+(run_id*dockpars_pop_size+entity_id)*GENOTYPE_LENGTH_IN_GLOBMEM,
                              		   dockpars_num_of_genes, 0);
129

130
	// Asynchronous copy should be finished by here
131
	wait_group_events(1, &ev);
132

133
  	// -----------------------------------------------------------------------------
134
135
	// -----------------------------------------------------------------------------
	// -----------------------------------------------------------------------------
Leonardo Solis's avatar
Leonardo Solis committed
136
	           
137
	// Partial results of the gradient step
Leonardo Solis's avatar
Leonardo Solis committed
138
139
140
	__local float gradient          [ACTUAL_GENOTYPE_LENGTH];
	__local float candidate_energy;
	__local float candidate_genotype[ACTUAL_GENOTYPE_LENGTH];
141
142
143
144
145

	// -------------------------------------------------------------------
	// Calculate gradients (forces) for intermolecular energy
	// Derived from autodockdev/maps.py
	// -------------------------------------------------------------------
146
	// Gradient of the intermolecular energy per each ligand atom
lvs's avatar
lvs committed
147
	// Also used to store the accummulated gradient per each ligand atom
148
149
150
	__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];
151

152
	// Gradient of the intramolecular energy per each ligand atom
153
154
155
	__local float gradient_intra_x[MAX_NUM_OF_ATOMS];
	__local float gradient_intra_y[MAX_NUM_OF_ATOMS];
	__local float gradient_intra_z[MAX_NUM_OF_ATOMS];
156

157
	// Ligand-atom position and partial energies
158
159
160
161
162
	__local float calc_coords_x[MAX_NUM_OF_ATOMS];
	__local float calc_coords_y[MAX_NUM_OF_ATOMS];
	__local float calc_coords_z[MAX_NUM_OF_ATOMS];
	__local float partial_energies[NUM_OF_THREADS_PER_BLOCK];

lvs's avatar
lvs committed
163
	#if defined (DEBUG_ENERGY_KERNEL)
164
165
	__local float partial_interE[NUM_OF_THREADS_PER_BLOCK];
	__local float partial_intraE[NUM_OF_THREADS_PER_BLOCK];
166
	#endif
167

168
	// Enable this for debugging SD from a defined initial genotype
169
	#if defined (DEBUG_INITIAL_2BRT)
lvs's avatar
lvs committed
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
	if (get_local_id(0) == 0) {
		// 2brt
		genotype[0] = 24.093334;
		genotype[1] = 24.658667;
		genotype[2] = 24.210667;
		genotype[3] = 50.0;
		genotype[4] = 50.0;
		genotype[5] = 50.0;
		genotype[6] = 0.0f;
		genotype[7] = 0.0f;
		genotype[8] = 0.0f;
		genotype[9] = 0.0f;
		genotype[10] = 0.0f;
		genotype[11] = 0.0f;
		genotype[12] = 0.0f;
		genotype[13] = 0.0f;
		genotype[14] = 0.0f;
		genotype[15] = 0.0f;
		genotype[16] = 0.0f;
		genotype[17] = 0.0f;
		genotype[18] = 0.0f;
		genotype[19] = 0.0f;
		genotype[20] = 0.0f;
	}
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
	// Evaluating candidate
	barrier(CLK_LOCAL_MEM_FENCE);

	// =============================================================
	gpu_calc_energy(dockpars_rotbondlist_length,
			dockpars_num_of_atoms,
			dockpars_gridsize_x,
			dockpars_gridsize_y,
			dockpars_gridsize_z,
							    	// g1 = gridsize_x
			dockpars_gridsize_x_times_y, 		// g2 = gridsize_x * gridsize_y
			dockpars_gridsize_x_times_y_times_z,	// g3 = gridsize_x * gridsize_y * gridsize_z
			dockpars_fgrids,
			dockpars_num_of_atypes,
			dockpars_num_of_intraE_contributors,
			dockpars_grid_spacing,
			dockpars_coeff_elec,
			dockpars_qasp,
			dockpars_coeff_desolv,
lvs's avatar
lvs committed
213
214
			dockpars_smooth,

lvs's avatar
lvs committed
215
			genotype, /*WARNING: calculating the energy of the hardcoded genotype*/
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
			&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,
			#if defined (DEBUG_ENERGY_KERNEL)
			partial_interE,
			partial_intraE,
			#endif
#if 0
			true,
#endif
lvs's avatar
lvs committed
233
234
235
236
237
			kerconst_interintra,
			kerconst_intracontrib,
			kerconst_intra,
			kerconst_rotlist,
			kerconst_conform			
238
239
240
			);
	// =============================================================

241
242
243
	// WARNING: hardcoded has priority over LGA genotype.
	// That means, if DEBUG_INITIAL_2BRT is defined, then
	// LGA genotype is not used (only for debugging purposes)
244
245
	if (get_local_id(0) == 0)
	{
246
247
248
249
		printf("\n");
		printf("%20s \n", "hardcoded genotype: ");
		printf("%20s %.6f\n", "initial energy: ", energy);
		printf("%20s %.6f\n\n", "initial stepsize: ", stepsize);		
250
251
252
253
	}
	barrier(CLK_LOCAL_MEM_FENCE);
	#endif

254
	// Perform gradient-descent iterations
255

256
	#if 0
Leonardo Solis's avatar
Leonardo Solis committed
257
258
259
260
261
262
263
264
265
266
	// 7cpa
	float grid_center_x = 49.836;
	float grid_center_y = 17.609;
	float grid_center_z = 36.272;
	float ligand_center_x = 49.2216976744186;
	float ligand_center_y = 17.793953488372097;
	float ligand_center_z = 36.503837209302326;
	float shoemake_gene_u1 = 0.02;
	float shoemake_gene_u2 = 0.23;
	float shoemake_gene_u3 = 0.95;
267
	#endif
Leonardo Solis's avatar
Leonardo Solis committed
268

269
270
271
272
273
274
275
276
277
	#if 0
	// 3tmn
	float grid_center_x = 52.340;
	float grid_center_y = 15.029;
	float grid_center_z = -2.932;
	float ligand_center_x = 52.22740741;
	float ligand_center_y = 15.51751852;
	float ligand_center_z = -2.40896296;
	#endif
278

lvs's avatar
lvs committed
279
	// Calculating maximum possible stepsize (alpha)
lvs's avatar
lvs committed
280
	__local float max_trans_grad, max_rota_grad, max_tors_grad;
lvs's avatar
lvs committed
281
282
283
	__local float max_trans_stepsize, max_rota_stepsize, max_tors_stepsize;
	__local float max_stepsize;

lvs's avatar
lvs committed
284
285
	// Storing torsion gradients here
	__local float torsions_gradient[ACTUAL_GENOTYPE_LENGTH];
lvs's avatar
lvs committed
286

Leonardo Solis's avatar
Leonardo Solis committed
287
288
289
	// The termination criteria is based on 
	// a maximum number of iterations, and
	// the minimum step size allowed for single-floating point numbers 
lvs's avatar
lvs committed
290
	// (IEEE-754 single float has a precision of about 6 decimal digits)
291
	do {
292
		#if 0
293
		// Specific input genotypes for a ligand with no rotatable bonds (1ac8).
Leonardo Solis's avatar
Leonardo Solis committed
294
		// Translation genes must be expressed in grids in OCLADock (genotype [0|1|2]).
295
296
		// However, for testing purposes, 
		// we start using translation values in real space (Angstrom): {31.79575, 93.743875, 47.699875}
Leonardo Solis's avatar
Leonardo Solis committed
297
		// Rotation genes are expresed in the Shoemake space: genotype [3|4|5]
298
		// xyz_gene_gridspace = gridcenter_gridspace + (input_gene_realspace - gridcenter_realspace)/gridsize
299
300

		// 1ac8				
Leonardo Solis's avatar
Leonardo Solis committed
301
302
303
304
305
306
		genotype[0] = 30 + (31.79575  - 31.924) / 0.375;
		genotype[1] = 30 + (93.743875 - 93.444) / 0.375;
		genotype[2] = 30 + (47.699875 - 47.924) / 0.375;
		genotype[3] = 0.1f;
		genotype[4] = 0.5f;
		genotype[5] = 0.9f;
307
		#endif
308

309
		#if 0
310
		// 3tmn
Leonardo Solis's avatar
Leonardo Solis committed
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
		genotype[0] = 30 + (ligand_center_x - grid_center_x) / 0.375;
		genotype[1] = 30 + (ligand_center_y - grid_center_y) / 0.375;
		genotype[2] = 30 + (ligand_center_z - grid_center_z) / 0.375;
		genotype[3] = shoemake_gene_u1;
		genotype[4] = shoemake_gene_u2;
		genotype[5] = shoemake_gene_u3;
		genotype[6] = 0.0f;
		genotype[7] = 0.0f;
		genotype[8] = 0.0f;
		genotype[9] = 0.0f;
		genotype[10] = 0.0f;
		genotype[11] = 0.0f;
		genotype[12] = 0.0f;
		genotype[13] = 0.0f;
		genotype[14] = 0.0f;
		genotype[15] = 0.0f;
		genotype[16] = 0.0f;
		genotype[17] = 0.0f;
		genotype[18] = 0.0f;
		genotype[19] = 0.0f;
		genotype[20] = 0.0f;
332
		#endif
333

334
335
		#if 0
		// 2j5s
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
		genotype[0] = 28.464;
		genotype[1] = 25.792762;
		genotype[2] = 23.740571;
		genotype[3] = 50.0;
		genotype[4] = 50.0;
		genotype[5] = 50.0;
		genotype[6] = 0.0f;
		genotype[7] = 0.0f;
		genotype[8] = 0.0f;
		genotype[9] = 0.0f;
		genotype[10] = 0.0f;
		genotype[11] = 0.0f;
		genotype[12] = 0.0f;
		genotype[13] = 0.0f;
		genotype[14] = 0.0f;
		genotype[15] = 0.0f;
		genotype[16] = 0.0f;
		genotype[17] = 0.0f;
		genotype[18] = 0.0f;
		genotype[19] = 0.0f;
		genotype[20] = 0.0f;
357
358
		#endif

lvs's avatar
lvs committed
359
		#if 0
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
		// 2brt
		genotype[0] = 24.093334;
		genotype[1] = 24.658667;
		genotype[2] = 24.210667;
		genotype[3] = 50.0;
		genotype[4] = 50.0;
		genotype[5] = 50.0;
		genotype[6] = 0.0f;
		genotype[7] = 0.0f;
		genotype[8] = 0.0f;
		genotype[9] = 0.0f;
		genotype[10] = 0.0f;
		genotype[11] = 0.0f;
		genotype[12] = 0.0f;
		genotype[13] = 0.0f;
		genotype[14] = 0.0f;
		genotype[15] = 0.0f;
		genotype[16] = 0.0f;
		genotype[17] = 0.0f;
		genotype[18] = 0.0f;
		genotype[19] = 0.0f;
		genotype[20] = 0.0f;
		#endif
383

384
385
386
387
388
389
390
391
392
393
394
395
396
		// Printing number of stepest-descent iterations
		#if defined (DEBUG_MINIMIZER) 
		if (get_local_id(0) == 0) {
			printf("%s\n", "----------------------------------------------------------");	
		}
		#endif
		
		#if defined (DEBUG_MINIMIZER) || defined (PRINT_MINIMIZER_ENERGY_EVOLUTION)
		if (get_local_id(0) == 0) {
			printf("%-15s %-3u ", "# SD iteration: ", iteration_cnt);
		}
		#endif

397
		// Calculating gradient
Leonardo Solis's avatar
Leonardo Solis committed
398
399
		barrier(CLK_LOCAL_MEM_FENCE);

400
401
402
403
404
405
406
		// =============================================================
		gpu_calc_gradient(
				dockpars_rotbondlist_length,
				dockpars_num_of_atoms,
				dockpars_gridsize_x,
				dockpars_gridsize_y,
				dockpars_gridsize_z,
407
408
409
								    	// g1 = gridsize_x
				dockpars_gridsize_x_times_y, 		// g2 = gridsize_x * gridsize_y
				dockpars_gridsize_x_times_y_times_z,	// g3 = gridsize_x * gridsize_y * gridsize_z
410
411
412
413
414
415
416
				dockpars_fgrids,
				dockpars_num_of_atypes,
				dockpars_num_of_intraE_contributors,
				dockpars_grid_spacing,
				dockpars_coeff_elec,
				dockpars_qasp,
				dockpars_coeff_desolv,
lvs's avatar
lvs committed
417
418
				dockpars_smooth,

419
420
421
422
				// 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.
Leonardo Solis's avatar
Leonardo Solis committed
423
424
				genotype,
				&energy,
Leonardo Solis's avatar
Leonardo Solis committed
425
426
				&run_id,

427
428
429
430
				calc_coords_x,
				calc_coords_y,
				calc_coords_z,

lvs's avatar
lvs committed
431
432
433
434
435
436
				kerconst_interintra,
				kerconst_intracontrib,
				kerconst_intra,
				kerconst_rotlist,
				kerconst_conform
				,
437
438
439
				rotbonds_const,
				rotbonds_atoms_const,
				num_rotating_atoms_per_rotbond_const
440
441
442
443
				,
	     			angle_const,
	     			dependence_on_theta_const,
	     			dependence_on_rotangle_const
444
445
446
447
448
449
450
451
			 	// Gradient-related arguments
			 	// Calculate gradients (forces) for intermolecular energy
			 	// Derived from autodockdev/maps.py
				,
				dockpars_num_of_genes,
				gradient_inter_x,
				gradient_inter_y,
				gradient_inter_z,
452
453
454
				gradient_intra_x,
				gradient_intra_y,
				gradient_intra_z,
Leonardo Solis's avatar
Leonardo Solis committed
455
				gradient
456
457
				);
		// =============================================================
458

459
460
461
		// This could be enabled back for double checking
		#if 0
		#if defined (DEBUG_ENERGY_KERNEL5)	
lvs's avatar
lvs committed
462
		if (/*(get_group_id(0) == 0) &&*/ (get_local_id(0) == 0)) {
463
464
		
			#if defined (PRINT_GENES_AND_GRADS)
lvs's avatar
lvs committed
465
			for(uint i = 0; i < dockpars_num_of_genes; i++) {
466
467
468
469
470
				if (i == 0) {
					printf("\n%s\n", "----------------------------------------------------------");
					printf("%13s %13s %5s %15s %15s\n", "gene_id", "gene.value", "|", "gene.grad", "(autodockdevpy units)");
				}
				printf("%13u %13.6f %5s %15.6f %15.6f\n", i, genotype[i], "|", gradient[i], (i<3)? (gradient[i]/0.375f):(gradient[i]*180.0f/PI_FLOAT));
lvs's avatar
lvs committed
471
			}
472
			#endif
473

474
			#if defined (PRINT_ATOMIC_COORDS)
475
476
477
			for(uint i = 0; i < dockpars_num_of_atoms; i++) {
				if (i == 0) {
					printf("\n%s\n", "----------------------------------------------------------");
478
					printf("%s\n", "Coordinates calculated by calcgradient.cl");
lvs's avatar
lvs committed
479
					printf("%12s %12s %12s %12s\n", "atom_id", "coords.x", "coords.y", "coords.z");
480
				}
lvs's avatar
lvs committed
481
				printf("%12u %12.6f %12.6f %12.6f\n", i, calc_coords_x[i], calc_coords_y[i], calc_coords_z[i]);
482
483
			}
			printf("\n");
484
			#endif
lvs's avatar
lvs committed
485
		}
486
487
		#endif
		#endif
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517

		if (get_local_id(0) == 0) {
			// Finding maximum of the absolute value for the three translation gradients
			max_trans_grad = fmax(fabs(gradient[0]), fabs(gradient[1]));
			max_trans_grad = fmax(max_trans_grad, fabs(gradient[2]));

			// MAX_DEV_TRANSLATION needs to be expressed in grid size first
			max_trans_stepsize = native_divide(native_divide(MAX_DEV_TRANSLATION, dockpars_grid_spacing), max_trans_grad);

			// Finding maximum of the absolute value for the three rotation gradients
			max_rota_grad = fmax(fabs(gradient[3]), fabs(gradient[4]));	
			max_rota_grad = fmax(max_rota_grad, fabs(gradient[5]));	

			// Note that MAX_DEV_ROTATION is already expressed approprietly
			max_rota_stepsize = native_divide(MAX_DEV_ROTATION, max_rota_grad);
		}

		// Copying torsions genes
		for(uint i = get_local_id(0); 
			 i < dockpars_num_of_genes-6; 
			 i+= NUM_OF_THREADS_PER_BLOCK) {
			torsions_gradient[i] = fabs(gradient[i+6]);
		}
		barrier(CLK_LOCAL_MEM_FENCE);

		// Calculating maximum absolute torsional gene
		// https://stackoverflow.com/questions/36465581/opencl-find-max-in-array
		for (uint i=(dockpars_num_of_genes-6)/2; i>=1; i/=2){
			if (get_local_id(0) < i) {

518
			// This could be enabled back for details
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
			#if 0
			#if defined (DEBUG_MINIMIZER)
			printf("---====--- %u %u %10.10f %-0.10f\n", i, get_local_id(0), torsions_gradient[get_local_id(0)], torsions_gradient[get_local_id(0) + i]);
			#endif
			#endif

				if (torsions_gradient[get_local_id(0)] < torsions_gradient[get_local_id(0) + i]) {
					torsions_gradient[get_local_id(0)] = torsions_gradient[get_local_id(0) + i];
				}
			}
			barrier(CLK_LOCAL_MEM_FENCE);
		}
		if (get_local_id(0) == 0) {
			max_tors_grad = torsions_gradient[get_local_id(0)];
			max_tors_stepsize = native_divide(MAX_DEV_TORSION, max_tors_grad);
		}
		barrier(CLK_LOCAL_MEM_FENCE);

		if (get_local_id(0) == 0) {
			// Calculating the maximum stepsize using previous three
			max_stepsize = fmin(max_trans_stepsize, max_rota_stepsize);
			max_stepsize = fmin(max_stepsize, max_tors_stepsize);

			// Capping the stepsize
			stepsize = fmin(stepsize, max_stepsize);
		}
545
		barrier(CLK_LOCAL_MEM_FENCE);
546
547
		
		for(uint i = get_local_id(0); i < dockpars_num_of_genes; i+= NUM_OF_THREADS_PER_BLOCK) {
548
	     		// Taking step
Leonardo Solis's avatar
Leonardo Solis committed
549
			candidate_genotype[i] = genotype[i] - stepsize * gradient[i];	
550
551

			#if defined (DEBUG_MINIMIZER)
lvs's avatar
lvs committed
552
553
			if (i == 0) {
				printf("\n%s\n", "After calculating gradients:");
554
				printf("%13s %13s %5s %15s %5s %20s\n", "gene_id", "gene", "|", "grad", "|", "cand.gene");
lvs's avatar
lvs committed
555
556
			}
			printf("%13u %13.6f %5s %15.6f %5s %18.6f\n", i, genotype[i], "|", gradient[i], "|", candidate_genotype[i]);
557

558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
			if (i == 0) {
				// This could be enabled back for double checking
				#if 0
				for(uint i = 0; i < dockpars_num_of_genes; i++) {
					if (i == 0) {
						printf("\n%s\n", "----------------------------------------------------------");
						printf("\n%s\n", "After calculating gradients:");
						printf("%13s %13s %5s %15s %20s\n", "gene_id", "gene", "|", "grad", " grad (devpy units)");
					}
					printf("%13u %13.6f %5s %15.6f %18.6f\n", i, genotype[i], "|", gradient[i], (i<3)? (gradient[i]/0.375f):(gradient[i]*180.0f/PI_FLOAT));
				}
				#endif

				printf("\n");
				printf("%20s %10.6f\n", "max_trans_grad: ", max_trans_grad);
				printf("%20s %10.6f\n", "max_rota_grad: ", max_rota_grad);
				printf("%20s %10.6f\n", "max_tors_grad: ", max_tors_grad);

				printf("\n");
				printf("%20s %10.6f\n", "max_trans_stepsize: ", max_trans_stepsize);
				printf("%20s %10.6f\n", "max_rota_stepsize: " , max_rota_stepsize);
				printf("%20s %10.6f\n", "max_tors_stepsize: " , max_tors_stepsize);

				printf("\n");
				printf("%20s %10.6f\n", "max_stepsize: ", max_stepsize);
				printf("%20s %10.6f\n", "stepsize: ", stepsize);
			}
lvs's avatar
lvs committed
585
			#endif
586
	   	}
587
		
588
		// Evaluating candidate
589
		barrier(CLK_LOCAL_MEM_FENCE);
590

Leonardo Solis's avatar
Leonardo Solis committed
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
		// =============================================================
		gpu_calc_energy(dockpars_rotbondlist_length,
				dockpars_num_of_atoms,
				dockpars_gridsize_x,
				dockpars_gridsize_y,
				dockpars_gridsize_z,
								    	// g1 = gridsize_x
				dockpars_gridsize_x_times_y, 		// g2 = gridsize_x * gridsize_y
				dockpars_gridsize_x_times_y_times_z,	// g3 = gridsize_x * gridsize_y * gridsize_z
				dockpars_fgrids,
				dockpars_num_of_atypes,
				dockpars_num_of_intraE_contributors,
				dockpars_grid_spacing,
				dockpars_coeff_elec,
				dockpars_qasp,
				dockpars_coeff_desolv,
lvs's avatar
lvs committed
607
608
				dockpars_smooth,

lvs's avatar
lvs committed
609
				candidate_genotype, /*genotype,*/ /*WARNING: use "genotype" ONLY to reproduce results*/
Leonardo Solis's avatar
Leonardo Solis committed
610
611
612
613
614
615
616
617
618
619
				&candidate_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,
lvs's avatar
lvs committed
620
				#if defined (DEBUG_ENERGY_KERNEL)
Leonardo Solis's avatar
Leonardo Solis committed
621
622
623
				partial_interE,
				partial_intraE,
				#endif
lvs's avatar
lvs committed
624
625
626
#if 0
				true,
#endif
lvs's avatar
lvs committed
627
628
629
630
631
				kerconst_interintra,
				kerconst_intracontrib,
				kerconst_intra,
				kerconst_rotlist,
				kerconst_conform
Leonardo Solis's avatar
Leonardo Solis committed
632
633
				);
		// =============================================================
634

635
		#if defined (DEBUG_ENERGY_KERNEL5)
lvs's avatar
lvs committed
636
		if (/*(get_group_id(0) == 0) &&*/ (get_local_id(0) == 0)) {
637
			#if defined (PRINT_ENERGIES)
638
639
640
641
			printf("\n");
			printf("%-10s %-10.6f \n", "intra: ",  partial_intraE[0]);
			printf("%-10s %-10.6f \n", "grids: ",  partial_interE[0]);
			printf("%-10s %-10.6f \n", "Energy: ", (partial_intraE[0] + partial_interE[0]));
642
			#endif
643

644
			#if defined (PRINT_GENES_AND_GRADS)
lvs's avatar
lvs committed
645
			for(uint i = 0; i < dockpars_num_of_genes; i++) {
646
647
				if (i == 0) {
					printf("\n%s\n", "----------------------------------------------------------");
lvs's avatar
lvs committed
648
					printf("%13s %13s %5s %15s %15s\n", "gene_id", "cand-gene.value"/* "gene.value"*/, "|", "gene.grad", "(autodockdevpy units)");
649
				}
lvs's avatar
lvs committed
650
				printf("%13u %13.6f %5s %15.6f %15.6f\n", i, candidate_genotype[i] /*genotype[i]*/, "|", gradient[i], (i<3)? (gradient[i]/0.375f):(gradient[i]*180.0f/PI_FLOAT));
651
			}
652
			#endif
653

654
			#if defined (PRINT_ATOMIC_COORDS)
655
656
657
			for(uint i = 0; i < dockpars_num_of_atoms; i++) {
				if (i == 0) {
					printf("\n%s\n", "----------------------------------------------------------");
658
					printf("%s\n", "Coordinates calculated by calcenergy.cl");
lvs's avatar
lvs committed
659
					printf("%12s %12s %12s %12s\n", "atom_id", "coords.x", "coords.y", "coords.z");
660
				}
lvs's avatar
lvs committed
661
				printf("%12u %12.6f %12.6f %12.6f\n", i, calc_coords_x[i], calc_coords_y[i], calc_coords_z[i]);
lvs's avatar
lvs committed
662
			}
663
			printf("\n");
664
			#endif
lvs's avatar
lvs committed
665
		}
666
		barrier(CLK_LOCAL_MEM_FENCE);
lvs's avatar
lvs committed
667
668
		#endif

lvs's avatar
lvs committed
669
670
671
		#if defined (DEBUG_MINIMIZER)
		if (get_local_id(0) == 0) {
			printf("\n");
672
673
674
			printf("%s\n", "After calculating energies:");
			printf("%13s %5s %15s\n", "energy", "|", "cand.energy");
			printf("%13.6f %5s %15.6f\n", energy, "|", candidate_energy);
lvs's avatar
lvs committed
675
676
677
678
679
			printf("\n");
		}
		barrier(CLK_LOCAL_MEM_FENCE);
		#endif

Leonardo Solis's avatar
Leonardo Solis committed
680
681
		// Checking if E(candidate_genotype) < E(genotype)
		if (candidate_energy < energy){
682
			
Leonardo Solis's avatar
Leonardo Solis committed
683
684
685
			for(uint i = get_local_id(0); 
			 	 i < dockpars_num_of_genes; 
		 	 	 i+= NUM_OF_THREADS_PER_BLOCK) {
lvs's avatar
lvs committed
686
687
688
689
690
691


				#if defined (DEBUG_MINIMIZER)
				//printf("(%-3u) %-15.7f %-10.7f %-10.7f %-10.7f\n", i, stepsize, genotype[i], gradient[i], candidate_genotype[i]);

				if (i == 0) {
692
					printf("%s\n", "Energy IMPROVED! ... then update genotype:");
lvs's avatar
lvs committed
693
694
695
696
					printf("%13s %13s %5s %15s\n", "gene_id", "old.gene", "|", "new.gene");
				}
				printf("%13u %13.6f %5s %15.6f\n", i, genotype[i], "|", candidate_genotype[i]);

697
				#endif
lvs's avatar
lvs committed
698
699
				if (i == 0) {
					#if defined (DEBUG_MINIMIZER)
700
					printf("\n%s\n", "Energy IMPROVED! ... then increase stepsize and update energy:");
lvs's avatar
lvs committed
701
702
					#endif

Leonardo Solis's avatar
Leonardo Solis committed
703
704
					// Increase stepsize
					stepsize *= STEP_INCREASE;
705
706
707

					// Updating energy
					energy = candidate_energy;
708
				}
Leonardo Solis's avatar
Leonardo Solis committed
709

Leonardo Solis's avatar
Leonardo Solis committed
710
711
				// Updating genotype
				genotype[i] = candidate_genotype[i];
712
			}
Leonardo Solis's avatar
Leonardo Solis committed
713
714
		}
		else { 
lvs's avatar
lvs committed
715
716
			#if defined (DEBUG_MINIMIZER)
			if (get_local_id(0) == 0) {
717
				printf("%s\n", "NO energy improvement! ... then decrease stepsize:");
lvs's avatar
lvs committed
718
719
720
			}
			#endif		

Leonardo Solis's avatar
Leonardo Solis committed
721
722
723
724
			if (get_local_id(0) == 0) {
				stepsize *= STEP_DECREASE;
			}
		}
725
		barrier(CLK_LOCAL_MEM_FENCE);
726

Leonardo Solis's avatar
Leonardo Solis committed
727
		// Updating number of stepest-descent iterations (energy evaluations)
728
		if (get_local_id(0) == 0) {
Leonardo Solis's avatar
Leonardo Solis committed
729
	    		iteration_cnt = iteration_cnt + 1;
730

lvs's avatar
lvs committed
731
			#if defined (DEBUG_MINIMIZER) || defined (PRINT_MINIMIZER_ENERGY_EVOLUTION)
732
			printf("%20s %10.10f %20s %10.6f\n", "new.stepsize: ", stepsize, "new.energy: ", energy);
733
			#endif
lvs's avatar
lvs committed
734

735
			#if defined (DEBUG_ENERGY_KERNEL5)
Leonardo Solis's avatar
Leonardo Solis committed
736
			printf("%-18s [%-5s]---{%-5s}   [%-10.7f]---{%-10.7f}\n", "-ENERGY-KERNEL5-", "GRIDS", "INTRA", partial_interE[0], partial_intraE[0]);
737
738
			#endif
		}
lvs's avatar
lvs committed
739
		barrier(CLK_LOCAL_MEM_FENCE);
740

lvs's avatar
lvs committed
741
  	} while ((iteration_cnt < dockpars_max_num_of_iters) && (stepsize > 1E-8));
742
743
	// -----------------------------------------------------------------------------
	// -----------------------------------------------------------------------------
744
	// -----------------------------------------------------------------------------
745

746
  	// Updating eval counter and energy
747
	if (get_local_id(0) == 0) {
Leonardo Solis's avatar
Leonardo Solis committed
748
		dockpars_evals_of_new_entities[run_id*dockpars_pop_size+entity_id] += iteration_cnt;
Leonardo Solis's avatar
Leonardo Solis committed
749
		dockpars_energies_next[run_id*dockpars_pop_size+entity_id] = energy;
750

lvs's avatar
lvs committed
751
		#if defined (DEBUG_MINIMIZER) || defined (PRINT_MINIMIZER_ENERGY_EVOLUTION)
752
753
754
		printf("\n");
		printf("Termination criteria: ( stepsize <= %10.10f ) OR ( #sd-iters >= %-3u )\n", 1E-8, dockpars_max_num_of_iters);
		printf("-------> End of SD minimization cycle, num of energy evals: %u, final energy: %.6f\n", iteration_cnt, energy);
755
		#endif
756
	}
757

758
759
760
761
	// Mapping torsion angles
	for (uint gene_counter = get_local_id(0);
	     	  gene_counter < dockpars_num_of_genes;
	          gene_counter+= NUM_OF_THREADS_PER_BLOCK) {
762
		   if (gene_counter >= 3) {
Leonardo Solis's avatar
Leonardo Solis committed
763
			    map_angle(&(genotype[gene_counter]));
764
		   }
765
766
	}

767
	// Updating old offspring in population
768
769
	barrier(CLK_LOCAL_MEM_FENCE);

770
771
772
773
774
775
	event_t ev2 = async_work_group_copy(dockpars_conformations_next+(run_id*dockpars_pop_size+entity_id)*GENOTYPE_LENGTH_IN_GLOBMEM,
			                    genotype,
			                    dockpars_num_of_genes, 0);

	// Asynchronous copy should be finished by here
	wait_group_events(1, &ev2);
776
}