Krnl_InterE.cl 13.5 KB
Newer Older
1
// --------------------------------------------------------------------------
2
3
4
5
6
// InterE calculates the intermolecular energy of a ligand given by 
// ligand xyz-positions, and a receptor represented as a grid. 
// The grid point values must be stored at the location which starts at GlobFgrids. 
// If an atom is remains outside the grid, 
// a very high value will be added to the current energy as a penalty. 
7
8
// Originally from: processligand.c
// --------------------------------------------------------------------------
lvs's avatar
lvs committed
9
/*
10
__kernel __attribute__ ((max_global_work_dim(0)))
lvs's avatar
lvs committed
11
12
*/
__kernel __attribute__ ((reqd_work_group_size(1,1,1)))
13
void Krnl_InterE(
14
15
16
             __global const float* restrict GlobFgrids,
 	     __constant float*     restrict KerConstStatic_atom_charges_const,
 	     __constant char*      restrict KerConstStatic_atom_types_const,
Leonardo Solis's avatar
Leonardo Solis committed
17
18
19
20
			    unsigned char                    DockConst_g1,
  			    unsigned int                     DockConst_g2,
			    unsigned int                     DockConst_g3,
   			    unsigned char                    DockConst_num_of_atoms,
Leonardo Solis's avatar
Leonardo Solis committed
21
22
23
			    float                   	     DockConst_gridsize_x_minus1,
			    float                    	     DockConst_gridsize_y_minus1,
			    float                    	     DockConst_gridsize_z_minus1,
24
	    		    unsigned int                     Host_mul_tmp2,
25
			    unsigned int                     Host_mul_tmp3
26
)
27
{
Leonardo Solis's avatar
Leonardo Solis committed
28
	char active = 0x01;
29

30
31
	__global const float* GlobFgrids2 = & GlobFgrids [Host_mul_tmp2];
	__global const float* GlobFgrids3 = & GlobFgrids [Host_mul_tmp3];
32

33
__attribute__((xcl_pipeline_loop))
Leonardo Solis's avatar
Leonardo Solis committed
34
LOOP_WHILE_INTERE_MAIN:
35
while(active) {
36

Leonardo Solis's avatar
Leonardo Solis committed
37
38
	char mode;

lvs's avatar
lvs committed
39
/*
40
41
42
43
44
45
46
47
	float3 __attribute__ ((
			      memory,
			      numbanks(2),
			      bankwidth(16),
			      singlepump,
			      numreadports(1),
			      numwriteports(1)
			    )) loc_coords[MAX_NUM_OF_ATOMS];
lvs's avatar
lvs committed
48
49
*/
	float3 loc_coords[MAX_NUM_OF_ATOMS];
50

51
52
53
54
	//printf("BEFORE In INTER CHANNEL\n");
	// --------------------------------------------------------------
	// Wait for ligand atomic coordinates in channel
	// --------------------------------------------------------------
Leonardo Solis's avatar
Leonardo Solis committed
55

56
	char actmode;
lvs's avatar
lvs committed
57
58
	read_pipe_block(chan_Conf2Intere_actmode, &actmode);
/*
59
	mem_fence(CLK_CHANNEL_MEM_FENCE);
lvs's avatar
lvs committed
60
*/
61
62
	active = actmode;
	mode   = actmode;
63

64
	__attribute__((xcl_pipeline_loop))
Leonardo Solis's avatar
Leonardo Solis committed
65
	LOOP_FOR_INTERE_READ_XYZ:
Leonardo Solis's avatar
Leonardo Solis committed
66
	for (uchar pipe_cnt=0; pipe_cnt<DockConst_num_of_atoms; pipe_cnt+=2) {
lvs's avatar
lvs committed
67
68
69
		float8 tmp;
		read_pipe_block(chan_Conf2Intere_xyz, &tmp);

Leonardo Solis's avatar
Leonardo Solis committed
70
71
72
73
74
75
		float3 tmp1 = {tmp.s0, tmp.s1, tmp.s2};
		float3 tmp2 = {tmp.s4, tmp.s5, tmp.s6};
		loc_coords[pipe_cnt] = tmp1;
		loc_coords[pipe_cnt+1] = tmp2;
	}

76
77
78
	// --------------------------------------------------------------
	//printf("AFTER In INTER CHANNEL\n");

79
	#if defined (DEBUG_ACTIVE_KERNEL)
Leonardo Solis's avatar
Leonardo Solis committed
80
	if (active == 0x00) {printf("	%-20s: %s\n", "Krnl_InterE", "must be disabled");}
81
	#endif
82

Leonardo Solis's avatar
Leonardo Solis committed
83
84
	float interE = 0.0f;

85
	// For each ligand atom
86
	__attribute__((xcl_pipeline_loop))
Leonardo Solis's avatar
Leonardo Solis committed
87
	LOOP_FOR_INTER_MAIN:
Leonardo Solis's avatar
Leonardo Solis committed
88
	for (uchar atom1_id=0; atom1_id<DockConst_num_of_atoms; atom1_id++)
89
	{
90
		char atom1_typeid = KerConstStatic_atom_types_const [atom1_id];
91

92
		float3 loc_coords_atid1 = loc_coords[atom1_id];
93

94
95
96
		float x = loc_coords_atid1.x;
		float y = loc_coords_atid1.y;
		float z = loc_coords_atid1.z;
97

98
		float q = KerConstStatic_atom_charges_const [atom1_id];
99

Leonardo Solis's avatar
Leonardo Solis committed
100
101
102
		float partialE1;
		float partialE2;
		float partialE3;
103

104
		// If the atom is outside of the grid
Leonardo Solis's avatar
Leonardo Solis committed
105
106
107
		if ((x < 0.0f) || (x >= DockConst_gridsize_x_minus1) || 
		    (y < 0.0f) || (y >= DockConst_gridsize_y_minus1) ||
		    (z < 0.0f) || (z >= DockConst_gridsize_z_minus1))	{
108

109
			// Penalty is 2^24 for each atom outside the grid
110
			/*
111
			interE += 16777216.0f; 
112
113
114
115
			*/
			partialE1 = 16777216.0f; 
			partialE2 = 0.0f;
			partialE3 = 0.0f;
116
117
118
		} 
		else 
		{
Leonardo Solis's avatar
Leonardo Solis committed
119
120
121
122
123
124
			int x_low  = convert_int(floor(x));
			int y_low  = convert_int(floor(y));
			int z_low  = convert_int(floor(z));
			int x_high = convert_int(ceil(x));	 
			int y_high = convert_int(ceil(y));
			int z_high = convert_int(ceil(z));
125

Leonardo Solis's avatar
Leonardo Solis committed
126
127
128
			float dx = x - x_low; 
			float dy = y - y_low; 
			float dz = z - z_low;
129
130
131

			// Calculates the weights for trilinear interpolation
			// based on the location of the point inside
Leonardo Solis's avatar
Leonardo Solis committed
132
			float weights [2][2][2];
133
134
135
136
137
138
139
140
141
			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;

142
			#if defined (DEBUG_KRNL_INTERE)
143
			printf("\n\nPartial results for atom with id %i:\n", atom1_id);
144
145
146
147
148
149
150
151
152
153
154
155
156
			printf("x_low = %d, x_high = %d, x_frac = %f\n", x_low, x_high, dx);
			printf("y_low = %d, y_high = %d, y_frac = %f\n", y_low, y_high, dy);
			printf("z_low = %d, z_high = %d, z_frac = %f\n\n", z_low, z_high, dz);
			printf("coeff(0,0,0) = %f\n", weights [0][0][0]);
			printf("coeff(1,0,0) = %f\n", weights [1][0][0]);
			printf("coeff(0,1,0) = %f\n", weights [0][1][0]);
			printf("coeff(1,1,0) = %f\n", weights [1][1][0]);
			printf("coeff(0,0,1) = %f\n", weights [0][0][1]);
			printf("coeff(1,0,1) = %f\n", weights [1][0][1]);
			printf("coeff(0,1,1) = %f\n", weights [0][1][1]);
			printf("coeff(1,1,1) = %f\n", weights [1][1][1]);
			#endif

157
			// Added temporal variables
Leonardo Solis's avatar
Leonardo Solis committed
158
			uint cube_000, cube_100, cube_010, cube_110, cube_001, cube_101, cube_011, cube_111;
159

Leonardo Solis's avatar
Leonardo Solis committed
160
161
162
163
164
			uint ylow_times_g1  = y_low  * DockConst_g1;	
			uint yhigh_times_g1 = y_high * DockConst_g1;
        	        uint zlow_times_g2  = z_low  * DockConst_g2;	
			uint zhigh_times_g2 = z_high * DockConst_g2;

165
166
167
168
169
170
171
172
        	        cube_000 = x_low  + ylow_times_g1  + zlow_times_g2;
        	        cube_100 = x_high + ylow_times_g1  + zlow_times_g2;
        	        cube_010 = x_low  + yhigh_times_g1 + zlow_times_g2;
        	        cube_110 = x_high + yhigh_times_g1 + zlow_times_g2;
        	        cube_001 = x_low  + ylow_times_g1  + zhigh_times_g2;
        	        cube_101 = x_high + ylow_times_g1  + zhigh_times_g2;
        	        cube_011 = x_low  + yhigh_times_g1 + zhigh_times_g2;
        	        cube_111 = x_high + yhigh_times_g1 + zhigh_times_g2;
Leonardo Solis's avatar
Leonardo Solis committed
173

174
			uint mul_tmp = atom1_typeid * DockConst_g3;
175

176
			// Energy contribution of the current grid type
Leonardo Solis's avatar
Leonardo Solis committed
177
			float cube [2][2][2];
178
179
180
181
182
183
184
185
	                cube [0][0][0] = GlobFgrids[cube_000 + mul_tmp];
        	        cube [1][0][0] = GlobFgrids[cube_100 + mul_tmp];
        	        cube [0][1][0] = GlobFgrids[cube_010 + mul_tmp];
        	        cube [1][1][0] = GlobFgrids[cube_110 + mul_tmp];
        	        cube [0][0][1] = GlobFgrids[cube_001 + mul_tmp];
        	        cube [1][0][1] = GlobFgrids[cube_101 + mul_tmp];
        	        cube [0][1][1] = GlobFgrids[cube_011 + mul_tmp];
        	        cube [1][1][1] = GlobFgrids[cube_111 + mul_tmp];
186
		
187
			#if defined (DEBUG_KRNL_INTERE)
188
189
190
191
192
193
194
195
196
197
198
			printf("Interpolation of van der Waals map:\n");
			printf("cube(0,0,0) = %f\n", cube [0][0][0]);
			printf("cube(1,0,0) = %f\n", cube [1][0][0]);
			printf("cube(0,1,0) = %f\n", cube [0][1][0]);
			printf("cube(1,1,0) = %f\n", cube [1][1][0]);
			printf("cube(0,0,1) = %f\n", cube [0][0][1]);
			printf("cube(1,0,1) = %f\n", cube [1][0][1]);
			printf("cube(0,1,1) = %f\n", cube [0][1][1]);
			printf("cube(1,1,1) = %f\n", cube [1][1][1]);
			#endif

199
200
201
202
203
204
205
206
207
			/*partialE1 = TRILININTERPOL(cube, weights);*/
			partialE1 = cube[0][0][0] * weights[0][0][0] +
				    cube[1][0][0] * weights[1][0][0] +
				    cube[0][1][0] * weights[0][1][0] +
				    cube[1][1][0] * weights[1][1][0] + 
				    cube[0][0][1] * weights[0][0][1] +
				    cube[1][0][1] * weights[1][0][1] + 
				    cube[0][1][1] * weights[0][1][1] +
				    cube[1][1][1] * weights[1][1][1];
208

209
			#if defined (DEBUG_KRNL_INTERE)
210
			printf("interpolated value = %f\n\n", TRILININTERPOL(cube, weights));
211
212
			#endif

213
			// Energy contribution of the electrostatic grid
Leonardo Solis's avatar
Leonardo Solis committed
214
215
216
			/*
			#if defined(SEPARATE_FGRID_INTERE)
			#else
217
			uint mul_tmp2 = Host_mul_tmp2;
Leonardo Solis's avatar
Leonardo Solis committed
218
219
			#endif
			*/
220
221
222
223
224
225
226
227
			cube [0][0][0] = GlobFgrids2[cube_000] /*GlobFgrids [Host_mul_tmp2 + cube_000]*/;
                        cube [1][0][0] = GlobFgrids2[cube_100] /*GlobFgrids [Host_mul_tmp2 + cube_100]*/;
                        cube [0][1][0] = GlobFgrids2[cube_010] /*GlobFgrids [Host_mul_tmp2 + cube_010]*/;
                        cube [1][1][0] = GlobFgrids2[cube_110] /*GlobFgrids [Host_mul_tmp2 + cube_110]*/;
                        cube [0][0][1] = GlobFgrids2[cube_001] /*GlobFgrids [Host_mul_tmp2 + cube_001]*/;
                        cube [1][0][1] = GlobFgrids2[cube_101] /*GlobFgrids [Host_mul_tmp2 + cube_101]*/;
                        cube [0][1][1] = GlobFgrids2[cube_011] /*GlobFgrids [Host_mul_tmp2 + cube_011]*/;
                        cube [1][1][1] = GlobFgrids2[cube_111] /*GlobFgrids [Host_mul_tmp2 + cube_111]*/;
228

229
			#if defined (DEBUG_KRNL_INTERE)
230
231
232
233
234
235
236
237
238
239
240
			printf("Interpolation of electrostatic map:\n");
			printf("cube(0,0,0) = %f\n", cube [0][0][0]);
			printf("cube(1,0,0) = %f\n", cube [1][0][0]);
			printf("cube(0,1,0) = %f\n", cube [0][1][0]);
			printf("cube(1,1,0) = %f\n", cube [1][1][0]);
			printf("cube(0,0,1) = %f\n", cube [0][0][1]);
			printf("cube(1,0,1) = %f\n", cube [1][0][1]);
			printf("cube(0,1,1) = %f\n", cube [0][1][1]);
			printf("cube(1,1,1) = %f\n", cube [1][1][1]);
			#endif

241
242
243
244
245
246
247
248
249
			/*partialE2 = q * TRILININTERPOL(cube, weights);*/
			partialE2 = q * (cube[0][0][0] * weights[0][0][0] +
				    	 cube[1][0][0] * weights[1][0][0] +
				    	 cube[0][1][0] * weights[0][1][0] +
				    	 cube[1][1][0] * weights[1][1][0] + 
				    	 cube[0][0][1] * weights[0][0][1] +
				    	 cube[1][0][1] * weights[1][0][1] + 
				    	 cube[0][1][1] * weights[0][1][1] +
				    	 cube[1][1][1] * weights[1][1][1]);
250
		
251
			#if defined (DEBUG_KRNL_INTERE)
Leonardo Solis's avatar
Leonardo Solis committed
252
			printf("interpolated value = %f, multiplied by q = %f\n\n", TRILININTERPOL(cube, weights), q*TRILININTERPOL(cube, weights));
253
254
			#endif

255
			// Energy contribution of the desolvation grid
Leonardo Solis's avatar
Leonardo Solis committed
256
257
258
			/*
			#if defined(SEPARATE_FGRID_INTERE)
			#else
259
			uint mul_tmp3 = Host_mul_tmp3;
Leonardo Solis's avatar
Leonardo Solis committed
260
261
			#endif
			*/
262

263
264
265
266
267
268
269
270
			cube [0][0][0] = GlobFgrids3[cube_000] /*GlobFgrids [Host_mul_tmp3 + cube_000]*/;
                        cube [1][0][0] = GlobFgrids3[cube_100] /*GlobFgrids [Host_mul_tmp3 + cube_100]*/;
                        cube [0][1][0] = GlobFgrids3[cube_010] /*GlobFgrids [Host_mul_tmp3 + cube_010]*/;
                        cube [1][1][0] = GlobFgrids3[cube_110] /*GlobFgrids [Host_mul_tmp3 + cube_110]*/;
                        cube [0][0][1] = GlobFgrids3[cube_001] /*GlobFgrids [Host_mul_tmp3 + cube_001]*/;
                        cube [1][0][1] = GlobFgrids3[cube_101] /*GlobFgrids [Host_mul_tmp3 + cube_101]*/;
                        cube [0][1][1] = GlobFgrids3[cube_011] /*GlobFgrids [Host_mul_tmp3 + cube_011]*/;
                        cube [1][1][1] = GlobFgrids3[cube_111] /*GlobFgrids [Host_mul_tmp3 + cube_111]*/;
271

272
			#if defined (DEBUG_KRNL_INTERE)
273
274
275
276
277
278
279
280
281
282
283
			printf("Interpolation of desolvation map:\n");
			printf("cube(0,0,0) = %f\n", cube [0][0][0]);
			printf("cube(1,0,0) = %f\n", cube [1][0][0]);
			printf("cube(0,1,0) = %f\n", cube [0][1][0]);
			printf("cube(1,1,0) = %f\n", cube [1][1][0]);
			printf("cube(0,0,1) = %f\n", cube [0][0][1]);
			printf("cube(1,0,1) = %f\n", cube [1][0][1]);
			printf("cube(0,1,1) = %f\n", cube [0][1][1]);
			printf("cube(1,1,1) = %f\n", cube [1][1][1]);
			#endif

284
285
286
287
288
289
290
291
292
			/*partialE3 = fabs(q) * TRILININTERPOL(cube, weights);*/
			partialE3 = fabs(q) * (cube[0][0][0] * weights[0][0][0] +
				    	       cube[1][0][0] * weights[1][0][0] +
				    	       cube[0][1][0] * weights[0][1][0] +
				    	       cube[1][1][0] * weights[1][1][0] + 
				    	       cube[0][0][1] * weights[0][0][1] +
				    	       cube[1][0][1] * weights[1][0][1] + 
				    	       cube[0][1][1] * weights[0][1][1] +
				    	       cube[1][1][1] * weights[1][1][1]);
293

294
			#if defined (DEBUG_KRNL_INTERE)
Leonardo Solis's avatar
Leonardo Solis committed
295
			printf("interpolated value = %f, multiplied by abs(q) = %f\n\n", TRILININTERPOL(cube, weights), fabs(q) * trilin_interpol(cube, weights));
296
			printf("Current value of intermolecular energy = %f\n\n\n", interE);
297
298
			#endif
		}
299
300

		interE += partialE1 + partialE2 + partialE3;
Leonardo Solis's avatar
Leonardo Solis committed
301
	} // End of atom1_id for-loop
302

303
304
305
	// --------------------------------------------------------------
	// Send intermolecular energy to chanel
	// --------------------------------------------------------------
306
307
	float final_interE = interE;

308
	switch (mode) {
Leonardo Solis's avatar
Leonardo Solis committed
309
310
		case 'I':  write_pipe_block(chan_Intere2StoreIC_intere,     &final_interE); break;
		case 'G':  write_pipe_block(chan_Intere2StoreGG_intere,     &final_interE); break;
lvs's avatar
lvs committed
311
312
313
314
315
316
317
318
319
		case 0x01: write_pipe_block(chan_Intere2StoreLS_LS1_intere, &final_interE); break;
		case 0x02: write_pipe_block(chan_Intere2StoreLS_LS2_intere, &final_interE); break;
		case 0x03: write_pipe_block(chan_Intere2StoreLS_LS3_intere, &final_interE); break;
		case 0x04: write_pipe_block(chan_Intere2StoreLS_LS4_intere, &final_interE); break;
		case 0x05: write_pipe_block(chan_Intere2StoreLS_LS5_intere, &final_interE); break;
		case 0x06: write_pipe_block(chan_Intere2StoreLS_LS6_intere, &final_interE); break;
		case 0x07: write_pipe_block(chan_Intere2StoreLS_LS7_intere, &final_interE); break;
		case 0x08: write_pipe_block(chan_Intere2StoreLS_LS8_intere, &final_interE); break;
		case 0x09: write_pipe_block(chan_Intere2StoreLS_LS9_intere, &final_interE); break;
320
	}
321
	// --------------------------------------------------------------
Leonardo Solis's avatar
Leonardo Solis committed
322
323
324
325

	#if defined (DEBUG_KRNL_INTERE)
	printf("AFTER Out INTERE CHANNEL\n");
	#endif
326
 	
327
} // End of while(active)
328
329
330
331

	#if defined (DEBUG_ACTIVE_KERNEL)
	printf("	%-20s: %s\n", "Krnl_InterE", "disabled");
	#endif
332
333
334
}
// --------------------------------------------------------------------------
// --------------------------------------------------------------------------
335