Fixed errors in DrudeLangevinIntegrator in mixed or double precision mode

plugins/drude/platforms/opencl/src/kernels/drudeLangevin.cl

@@ -13,10 +13,12 @@ __kernel void integrateDrudeLangevinPart1(__global mixed4* restrict velm, __glob
         int index = normalParticles[i];
         mixed4 velocity = velm[index];
         if (velocity.w != 0) {
-            mixed sqrtInvMass = SQRT(velocity.w);
+            mixed sqrtInvMass = sqrt(velocity.w);
             float4 rand = random[randomIndex+index];
             real4 f = force[index];
-            velocity.xyz = vscale*velocity.xyz + fscale*velocity.w*f.xyz + noisescale*sqrtInvMass*rand.xyz;
+            velocity.x = vscale*velocity.x + fscale*velocity.w*f.x + noisescale*sqrtInvMass*rand.x;
+            velocity.y = vscale*velocity.y + fscale*velocity.w*f.y + noisescale*sqrtInvMass*rand.y;
+            velocity.z = vscale*velocity.z + fscale*velocity.w*f.z + noisescale*sqrtInvMass*rand.z;
             velm[index] = velocity;
             posDelta[index] = (mixed4) (stepSize*velocity.x, stepSize*velocity.y, stepSize*velocity.z, 0);
         }
@@ -29,24 +31,28 @@ __kernel void integrateDrudeLangevinPart1(__global mixed4* restrict velm, __glob
         int2 particles = pairParticles[i];
         mixed4 velocity1 = velm[particles.x];
         mixed4 velocity2 = velm[particles.y];
-        mixed mass1 = RECIP(velocity1.w);
-        mixed mass2 = RECIP(velocity2.w);
-        mixed invTotalMass = RECIP(mass1+mass2);
+        mixed mass1 = 1/velocity1.w;
+        mixed mass2 = 1/velocity2.w;
+        mixed invTotalMass = 1/(mass1+mass2);
         mixed invReducedMass = (mass1+mass2)*velocity1.w*velocity2.w;
         mixed mass1fract = invTotalMass*mass1;
         mixed mass2fract = invTotalMass*mass2;
-        mixed sqrtInvTotalMass = SQRT(invTotalMass);
-        mixed sqrtInvReducedMass = SQRT(invReducedMass);
+        mixed sqrtInvTotalMass = sqrt(invTotalMass);
+        mixed sqrtInvReducedMass = sqrt(invReducedMass);
         mixed4 cmVel = velocity1*mass1fract+velocity2*mass2fract;
         mixed4 relVel = velocity2-velocity1;
-        mixed4 force1 = force[particles.x];
-        mixed4 force2 = force[particles.y];
+        mixed4 force1 = convert_mixed4(force[particles.x]);
+        mixed4 force2 = convert_mixed4(force[particles.y]);
         mixed4 cmForce = force1+force2;
         mixed4 relForce = force2*mass1fract - force1*mass2fract;
         float4 rand1 = random[randomIndex+2*i];
         float4 rand2 = random[randomIndex+2*i+1];
-        cmVel.xyz = vscale*cmVel.xyz + fscale*invTotalMass*cmForce.xyz + noisescale*sqrtInvTotalMass*rand1.xyz;
-        relVel.xyz = vscaleDrude*relVel.xyz + fscaleDrude*invReducedMass*relForce.xyz + noisescaleDrude*sqrtInvReducedMass*rand2.xyz;
+        cmVel.x = vscale*cmVel.x + fscale*invTotalMass*cmForce.x + noisescale*sqrtInvTotalMass*rand1.x;
+        cmVel.y = vscale*cmVel.y + fscale*invTotalMass*cmForce.y + noisescale*sqrtInvTotalMass*rand1.y;
+        cmVel.z = vscale*cmVel.z + fscale*invTotalMass*cmForce.z + noisescale*sqrtInvTotalMass*rand1.z;
+        relVel.x = vscaleDrude*relVel.x + fscaleDrude*invReducedMass*relForce.x + noisescaleDrude*sqrtInvReducedMass*rand2.x;
+        relVel.y = vscaleDrude*relVel.y + fscaleDrude*invReducedMass*relForce.y + noisescaleDrude*sqrtInvReducedMass*rand2.y;
+        relVel.z = vscaleDrude*relVel.z + fscaleDrude*invReducedMass*relForce.z + noisescaleDrude*sqrtInvReducedMass*rand2.z;
         velocity1.xyz = cmVel.xyz-relVel.xyz*mass2fract;
         velocity2.xyz = cmVel.xyz+relVel.xyz*mass1fract;
         velm[particles.x] = velocity1;