harmonicBondForce.cl

float4 delta = pos2-pos1;
float2 bondParams = PARAMS[index];
float r = SQRT(delta.x*delta.x + delta.y*delta.y + delta.z*delta.z);
float deltaIdeal = r-bondParams.x;
energy += 0.5f * bondParams.y*deltaIdeal*deltaIdeal;
float dEdR = bondParams.y * deltaIdeal;
dEdR = (r > 0.0f) ? (dEdR / r) : 0.0f;
delta.xyz *= dEdR;
float4 force1 = delta;
float4 force2 = -delta;