I have a grid of 30x30x30 periodic cells for which I need to calculate Lennard Jones Potential Force using CUDA. A global function calls this device function and supplies it with the data structure at each point and (x, y, z) position for each point.
__device__ void calculateForce(point a[][DATAYSIZE][DATAZSIZE], int idx, int idy, int idz){
a[idx][idy][idz].force_a[0] = 0;
a[idx][idy][idz].force_a[1] = 0;
a[idx][idy][idz].force_a[2] = 0;
a[idx][idy][idz].force_b[0] = 0;
a[idx][idy][idz].force_b[1] = 0;
a[idx][idy][idz].force_b[2] = 0;
int dx, dy, dz;
float fax = 0, fay = 0, faz = 0;
float fbx = 0, fby = 0, fbz = 0;
float e1 = 24 * epsilon_a/da;
float e2 = 24 * eab/dab;
float e3 = 24 * epsilon_b/db;
for (int i = idx-CUTOFF; i < idx+CUTOFF; i++){
for(int j = idy-CUTOFF; j < idy+CUTOFF; j++){
for(int k = idz-CUTOFF; k < idz+CUTOFF; k++){
dx = idx - i;
dy = idy - j;
dz = idz - k;
float distance = sqrtf(dx*dx+dy*dy+dz*dz);
if(distance > CUTOFF){
continue;
}
int x=i,y=j,z=k;
if(x>DATAXSIZE-1){
x -= DATAXSIZE;
}else if(x<0){
x += DATAXSIZE;
}
if(y>DATAYSIZE-1){
y -= DATAYSIZE;
}else if(y<0){
y += DATAYSIZE;
}
if(z>DATAZSIZE-1){
z -= DATAZSIZE;
}else if(z<0){
z += DATAZSIZE;
}
float sqa = da/distance;
float sqab = dab/distance;
float sqb = db/distance;
float faa = e1 * (2 * powf(sqa, 13) - powf(sqa, 7));
float fab = e2 *(2 * powf(sqab, 13) - powf(sqab, 7)) ;
float fbb = e3 * (2 * powf(sqb, 13) - powf(sqb, 7));
float na = (a[x][y][z].na * faa + a[x][y][z].nb * fab)/distance;
float nb = (a[x][y][z].na * fab + a[x][y][z].nb * fbb)/distance;
fax += na * dx;
fay += na * dy;
faz += na * dz;
fbx += nb * dx;
fby += nb * dy;
fbz += nb * dz;
}
}
}
a[idx][idy][idz].force_a[0] = fax;
a[idx][idy][idz].force_a[1] = fay;
a[idx][idy][idz].force_a[2] = faz;
a[idx][idy][idz].force_b[0] = fbx;
a[idx][idy][idz].force_b[1] = fby;
a[idx][idy][idz].force_b[2] = fbz;
}
For a cube of length (2 x CUTOFF radius), I calculate distance of each cell from the centre of the cube. If the distance is less than the CUTOFF radius, then I calculate Lennard Jones Potential Force. I have already minimised repetitive interactions for this code.
But here, addition assignment operations from line 55 to 61 seem to be extreme bottleneck. Is it possible to optimise it further ?