Fix mixed precision bead propagation

platforms/cuda/src/CudaKernels.cpp

@@ -8353,7 +8353,10 @@ void CudaApplyAndersenThermostatKernel::execute(ContextImpl& context) {
 void CudaNoseHooverChainKernel::initialize() {
     cu.setAsCurrent();
 
+    bool useDouble = cu.getUseDoublePrecision() || cu.getUseMixedPrecision();
+
     map<string, string> defines;
+    defines["MIXEDEXP"] = useDouble ? "exp" : "expf";
     defines["BEGIN_YS_LOOP"] = "for(const real & ys : {1}) {";
     defines["END_YS_LOOP"] = "}";
     CUmodule module = cu.createModule(CudaKernelSources::noseHooverChain, defines, "-std=c++11");

platforms/cuda/src/kernels/noseHooverChain.cu

@@ -25,11 +25,11 @@ extern "C" __global__ void propagateNoseHooverChain(mixed2* __restrict__ chainDa
             mixed wdt = ys * timeOverMTS;
             chainData[chainLength-1].y += 0.25f * wdt * chainForces[chainLength-1];
             for (int bead = chainLength - 2; bead >= 0; --bead) {
-                mixed aa = EXP(-0.125f * wdt * chainData[bead + 1].y);
+                mixed aa = MIXEDEXP(-0.125f * wdt * chainData[bead + 1].y);
                 chainData[bead].y = aa * (chainData[bead].y * aa + 0.25f * wdt * chainForces[bead]);
             }
             // update particle velocities
-            mixed aa = EXP(-0.5f * wdt * chainData[0].y);
+            mixed aa = MIXEDEXP(-0.5f * wdt * chainData[0].y);
             scale *= aa;
             // update the thermostat positions
             for (int bead = 0; bead < chainLength; ++bead) {
@@ -39,7 +39,7 @@ extern "C" __global__ void propagateNoseHooverChain(mixed2* __restrict__ chainDa
             chainForces[0] = (scale * scale * KE2 - numDOFs * kT) / chainMasses[0];
             // update thermostat velocities
             for (int bead = 0; bead < chainLength - 1; ++bead) {
-                mixed aa = EXP(-0.125f * wdt * chainData[bead + 1].y);
+                mixed aa = MIXEDEXP(-0.125f * wdt * chainData[bead + 1].y);
                 chainData[bead].y = aa * (aa * chainData[bead].y + 0.25f * wdt * chainForces[bead]);
                 chainForces[bead + 1] = (chainMasses[bead] * chainData[bead].y * chainData[bead].y - kT) / chainMasses[bead + 1];
             }