Krnl_InterE.cl 12.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
__kernel __attribute__ ((reqd_work_group_size(1,1,1)))
10
void Krnl_InterE(
11
12
13
             __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
14
15
16
17
			    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
18
19
20
			    float                   	     DockConst_gridsize_x_minus1,
			    float                    	     DockConst_gridsize_y_minus1,
			    float                    	     DockConst_gridsize_z_minus1,
21
	    		    unsigned int                     Host_mul_tmp2,
22
			    unsigned int                     Host_mul_tmp3
23
)
24
{
Leonardo Solis's avatar
Leonardo Solis committed
25
	char active = 0x01;
26

27
28
	__global const float* GlobFgrids2 = & GlobFgrids [Host_mul_tmp2];
	__global const float* GlobFgrids3 = & GlobFgrids [Host_mul_tmp3];
29

30
__attribute__((xcl_pipeline_loop))
Leonardo Solis's avatar
Leonardo Solis committed
31
LOOP_WHILE_INTERE_MAIN:
32
while(active) {
33

Leonardo Solis's avatar
Leonardo Solis committed
34
35
	char mode;

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

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

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

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

Leonardo Solis's avatar
Leonardo Solis committed
67
68
69
70
71
72
		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;
	}

73
74
75
	// --------------------------------------------------------------
	//printf("AFTER In INTER CHANNEL\n");

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

Leonardo Solis's avatar
Leonardo Solis committed
80
81
	float interE = 0.0f;

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

89
		float3 loc_coords_atid1 = loc_coords[atom1_id];
90

91
92
93
		float x = loc_coords_atid1.x;
		float y = loc_coords_atid1.y;
		float z = loc_coords_atid1.z;
94

95
		float q = KerConstStatic_atom_charges_const [atom1_id];
96

Leonardo Solis's avatar
Leonardo Solis committed
97
98
99
		float partialE1;
		float partialE2;
		float partialE3;
100

101
		// If the atom is outside of the grid
Leonardo Solis's avatar
Leonardo Solis committed
102
103
104
		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))	{
105

106
			// Penalty is 2^24 for each atom outside the grid
107
			/*
108
			interE += 16777216.0f; 
109
110
111
112
			*/
			partialE1 = 16777216.0f; 
			partialE2 = 0.0f;
			partialE3 = 0.0f;
113
114
115
		} 
		else 
		{
Leonardo Solis's avatar
Leonardo Solis committed
116
117
118
119
120
121
			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));
122

Leonardo Solis's avatar
Leonardo Solis committed
123
124
125
			float dx = x - x_low; 
			float dy = y - y_low; 
			float dz = z - z_low;
126
127
128

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

139
			#if defined (DEBUG_KRNL_INTERE)
140
			printf("\n\nPartial results for atom with id %i:\n", atom1_id);
141
142
143
144
145
146
147
148
149
150
151
152
153
			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

154
			// Added temporal variables
Leonardo Solis's avatar
Leonardo Solis committed
155
			uint cube_000, cube_100, cube_010, cube_110, cube_001, cube_101, cube_011, cube_111;
156

Leonardo Solis's avatar
Leonardo Solis committed
157
158
159
160
161
			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;

162
163
164
165
166
167
168
169
        	        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
170

171
			uint mul_tmp = atom1_typeid * DockConst_g3;
172

173
			// Energy contribution of the current grid type
Leonardo Solis's avatar
Leonardo Solis committed
174
			float cube [2][2][2];
175
176
177
178
179
180
181
182
	                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];
183
		
184
			#if defined (DEBUG_KRNL_INTERE)
185
186
187
188
189
190
191
192
193
194
195
			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

196
197
198
199
200
201
202
203
204
			/*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];
205

206
			#if defined (DEBUG_KRNL_INTERE)
207
			printf("interpolated value = %f\n\n", TRILININTERPOL(cube, weights));
208
209
			#endif

210
			// Energy contribution of the electrostatic grid
211
212
213
214
215
216
217
218
			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]*/;
219

220
			#if defined (DEBUG_KRNL_INTERE)
221
222
223
224
225
226
227
228
229
230
231
			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

232
233
234
235
236
237
238
239
240
			/*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]);
241
		
242
			#if defined (DEBUG_KRNL_INTERE)
Leonardo Solis's avatar
Leonardo Solis committed
243
			printf("interpolated value = %f, multiplied by q = %f\n\n", TRILININTERPOL(cube, weights), q*TRILININTERPOL(cube, weights));
244
245
			#endif

246
			// Energy contribution of the desolvation grid
247
248
249
250
251
252
253
254
			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]*/;
255

256
			#if defined (DEBUG_KRNL_INTERE)
257
258
259
260
261
262
263
264
265
266
267
			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

268
269
270
271
272
273
274
275
276
			/*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]);
277

278
			#if defined (DEBUG_KRNL_INTERE)
Leonardo Solis's avatar
Leonardo Solis committed
279
			printf("interpolated value = %f, multiplied by abs(q) = %f\n\n", TRILININTERPOL(cube, weights), fabs(q) * trilin_interpol(cube, weights));
280
			printf("Current value of intermolecular energy = %f\n\n\n", interE);
281
282
			#endif
		}
283
284

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

287
288
289
	// --------------------------------------------------------------
	// Send intermolecular energy to chanel
	// --------------------------------------------------------------
290
291
	float final_interE = interE;

292
	switch (mode) {
Leonardo Solis's avatar
Leonardo Solis committed
293
294
		case 'I':  write_pipe_block(chan_Intere2StoreIC_intere,     &final_interE); break;
		case 'G':  write_pipe_block(chan_Intere2StoreGG_intere,     &final_interE); break;
295
	}
296
	// --------------------------------------------------------------
Leonardo Solis's avatar
Leonardo Solis committed
297
298
299
300

	#if defined (DEBUG_KRNL_INTERE)
	printf("AFTER Out INTERE CHANNEL\n");
	#endif
301
 	
302
} // End of while(active)
303
304
305
306

	#if defined (DEBUG_ACTIVE_KERNEL)
	printf("	%-20s: %s\n", "Krnl_InterE", "disabled");
	#endif
307
308
309
}
// --------------------------------------------------------------------------
// --------------------------------------------------------------------------
310