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
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
	#if defined (DEBUG_INITIAL_2BRT)
	// 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;

	// 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
212
213
			dockpars_smooth,

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

240
241
242
	// 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)
243
244
	if (get_local_id(0) == 0)
	{
245
246
247
248
		printf("\n");
		printf("%20s \n", "hardcoded genotype: ");
		printf("%20s %.6f\n", "initial energy: ", energy);
		printf("%20s %.6f\n\n", "initial stepsize: ", stepsize);		
249
250
251
252
	}
	barrier(CLK_LOCAL_MEM_FENCE);
	#endif

253
	// Perform gradient-descent iterations
254

255
	#if 0
Leonardo Solis's avatar
Leonardo Solis committed
256
257
258
259
260
261
262
263
264
265
	// 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;
266
	#endif
Leonardo Solis's avatar
Leonardo Solis committed
267

268
269
270
271
272
273
274
275
276
	#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
277

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

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

Leonardo Solis's avatar
Leonardo Solis committed
286
287
288
	// 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
289
	// (IEEE-754 single float has a precision of about 6 decimal digits)
290
	do {
291
		#if 0
292
		// Specific input genotypes for a ligand with no rotatable bonds (1ac8).
Leonardo Solis's avatar
Leonardo Solis committed
293
		// Translation genes must be expressed in grids in OCLADock (genotype [0|1|2]).
294
295
		// 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
296
		// Rotation genes are expresed in the Shoemake space: genotype [3|4|5]
297
		// xyz_gene_gridspace = gridcenter_gridspace + (input_gene_realspace - gridcenter_realspace)/gridsize
298
299

		// 1ac8				
Leonardo Solis's avatar
Leonardo Solis committed
300
301
302
303
304
305
		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;
306
		#endif
307

308
		#if 0
309
		// 3tmn
Leonardo Solis's avatar
Leonardo Solis committed
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
		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;
331
		#endif
332

333
334
		#if 0
		// 2j5s
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
		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;
356
357
		#endif

lvs's avatar
lvs committed
358
		#if 0
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
		// 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
382

383
384
385
386
387
388
389
390
391
392
393
394
395
		// 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

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

399
400
401
402
403
404
405
		// =============================================================
		gpu_calc_gradient(
				dockpars_rotbondlist_length,
				dockpars_num_of_atoms,
				dockpars_gridsize_x,
				dockpars_gridsize_y,
				dockpars_gridsize_z,
406
407
408
								    	// 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
409
410
411
412
413
414
415
				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
416
417
				dockpars_smooth,

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

426
427
428
429
				calc_coords_x,
				calc_coords_y,
				calc_coords_z,

lvs's avatar
lvs committed
430
431
432
433
434
435
				kerconst_interintra,
				kerconst_intracontrib,
				kerconst_intra,
				kerconst_rotlist,
				kerconst_conform
				,
436
437
438
				rotbonds_const,
				rotbonds_atoms_const,
				num_rotating_atoms_per_rotbond_const
439
440
441
442
				,
	     			angle_const,
	     			dependence_on_theta_const,
	     			dependence_on_rotangle_const
443
444
445
446
447
448
449
450
			 	// 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,
451
452
453
				gradient_intra_x,
				gradient_intra_y,
				gradient_intra_z,
Leonardo Solis's avatar
Leonardo Solis committed
454
				gradient
455
456
				);
		// =============================================================
457

458
459
460
		// This could be enabled back for double checking
		#if 0
		#if defined (DEBUG_ENERGY_KERNEL5)	
lvs's avatar
lvs committed
461
		if (/*(get_group_id(0) == 0) &&*/ (get_local_id(0) == 0)) {
462
463
		
			#if defined (PRINT_GENES_AND_GRADS)
lvs's avatar
lvs committed
464
			for(uint i = 0; i < dockpars_num_of_genes; i++) {
465
466
467
468
469
				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
470
			}
471
			#endif
472

473
			#if defined (PRINT_ATOMIC_COORDS)
474
475
476
			for(uint i = 0; i < dockpars_num_of_atoms; i++) {
				if (i == 0) {
					printf("\n%s\n", "----------------------------------------------------------");
477
					printf("%s\n", "Coordinates calculated by calcgradient.cl");
lvs's avatar
lvs committed
478
					printf("%12s %12s %12s %12s\n", "atom_id", "coords.x", "coords.y", "coords.z");
479
				}
lvs's avatar
lvs committed
480
				printf("%12u %12.6f %12.6f %12.6f\n", i, calc_coords_x[i], calc_coords_y[i], calc_coords_z[i]);
481
482
			}
			printf("\n");
483
			#endif
lvs's avatar
lvs committed
484
		}
485
486
		#endif
		#endif
487
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

		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) {

517
			// This could be enabled back for details
518
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
			#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);
		}
544
		barrier(CLK_LOCAL_MEM_FENCE);
545
546
		
		for(uint i = get_local_id(0); i < dockpars_num_of_genes; i+= NUM_OF_THREADS_PER_BLOCK) {
547
	     		// Taking step
Leonardo Solis's avatar
Leonardo Solis committed
548
			candidate_genotype[i] = genotype[i] - stepsize * gradient[i];	
549
550

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

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
			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
584
			#endif
585
	   	}
586
		
587
		// Evaluating candidate
588
		barrier(CLK_LOCAL_MEM_FENCE);
589

Leonardo Solis's avatar
Leonardo Solis committed
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
		// =============================================================
		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
606
607
				dockpars_smooth,

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

634
		#if defined (DEBUG_ENERGY_KERNEL5)
lvs's avatar
lvs committed
635
		if (/*(get_group_id(0) == 0) &&*/ (get_local_id(0) == 0)) {
636
			#if defined (PRINT_ENERGIES)
637
638
639
640
			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]));
641
			#endif
642

643
			#if defined (PRINT_GENES_AND_GRADS)
lvs's avatar
lvs committed
644
			for(uint i = 0; i < dockpars_num_of_genes; i++) {
645
646
				if (i == 0) {
					printf("\n%s\n", "----------------------------------------------------------");
lvs's avatar
lvs committed
647
					printf("%13s %13s %5s %15s %15s\n", "gene_id", "cand-gene.value"/* "gene.value"*/, "|", "gene.grad", "(autodockdevpy units)");
648
				}
lvs's avatar
lvs committed
649
				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));
650
			}
651
			#endif
652

653
			#if defined (PRINT_ATOMIC_COORDS)
654
655
656
			for(uint i = 0; i < dockpars_num_of_atoms; i++) {
				if (i == 0) {
					printf("\n%s\n", "----------------------------------------------------------");
657
					printf("%s\n", "Coordinates calculated by calcenergy.cl");
lvs's avatar
lvs committed
658
					printf("%12s %12s %12s %12s\n", "atom_id", "coords.x", "coords.y", "coords.z");
659
				}
lvs's avatar
lvs committed
660
				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
661
			}
662
			printf("\n");
663
			#endif
lvs's avatar
lvs committed
664
		}
665
		barrier(CLK_LOCAL_MEM_FENCE);
lvs's avatar
lvs committed
666
667
		#endif

lvs's avatar
lvs committed
668
669
670
		#if defined (DEBUG_MINIMIZER)
		if (get_local_id(0) == 0) {
			printf("\n");
671
672
673
			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
674
675
676
677
678
			printf("\n");
		}
		barrier(CLK_LOCAL_MEM_FENCE);
		#endif

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


				#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) {
691
					printf("%s\n", "Energy IMPROVED! ... then update genotype:");
lvs's avatar
lvs committed
692
693
694
695
					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]);

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

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

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

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

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

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

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

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

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

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

lvs's avatar
lvs committed
750
		#if defined (DEBUG_MINIMIZER) || defined (PRINT_MINIMIZER_ENERGY_EVOLUTION)
751
752
753
		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);
754
		#endif
755
	}
756

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

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

769
770
771
772
773
774
	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);
775
}