         c = c * a;
         s12 = s12 * a;
         Scalar a11 = c*sb1*s1;
-        Scalar a22 = -sb2 * (2.0f*c0*s12 - c*(s1+s2));
+        Scalar a22 = -sb2 * (Scalar(2.0)*c0*s12 - c*(s1+s2));
         Scalar a33 = c*sb3*s2;
         Scalar a12 = -r12c1*(c1mag*c*s1 + c2mag*s12);
         Scalar a13 = -rb1*rb3*s12;

        Scalar b3mag = sqrt(b3mag2);

        Scalar ctmp = dab.x*dcb.x + dab.y*dcb.y + dab.z*dcb.z;
        Scalar r12c1 = 1.0f / (b1mag*b2mag);
        Scalar c1mag = ctmp * r12c1;

        ctmp = dcbm.x*ddc.x + dcbm.y*ddc.y + dcbm.z*ddc.z;
        Scalar r12c2 = 1.0f / (b2mag*b3mag);
        Scalar c2mag = ctmp * r12c2;

        // cos and sin of 2 angles and final c

        Scalar sin2 = 1.0f - c1mag*c1mag;
        if (sin2 < 0.0f) sin2 = 0.0f;
        Scalar sc1 = sqrtf(sin2);
        if (sc1 < SMALL) sc1 = SMALL;
        sc1 = 1.0f/sc1;

        sin2 = 1.0f - c2mag*c2mag;
        if (sin2 < 0.0f) sin2 = 0.0f;
        Scalar sc2 = sqrtf(sin2);
        if (sc2 < SMALL) sc2 = SMALL;
        sc2 = 1.0f/sc2;

        Scalar s1 = sc1 * sc1;
        Scalar s2 = sc2 * sc2;
        Scalar s12 = sc1 * sc2;
        Scalar c = (c0 + c1mag*c2mag) * s12;

        if (c > 1.0f) c = 1.0f;
        if (c < -1.0f) c = -1.0f;

        //phi
        Scalar phi = acosf(c);

        // and accumlate the energy/virial

        // compute 1/4 of the energy, 1/4 for each atom in the dihedral
        Scalar dihedral_eng = V*Scalar(1.0/4.0);

        // compute 1/4 of the virial, 1/4 for each atom in the dihedral
        // symmetrized version of virial tensor
