From 7a7325f6f71fa51ec06d4a10e040290f3d54e278 Mon Sep 17 00:00:00 2001
From: Michael Zingale <michael.zingale@stonybrook.edu>
Date: Sun, 10 Sep 2023 14:21:07 -0400
Subject: [PATCH 1/5] template out the plasma neutrino derivatives

---
 networks/rhs.H     |  10 +-
 neutrinos/sneut5.H | 564 +++++++++++++++++++++++++++------------------
 2 files changed, 350 insertions(+), 224 deletions(-)

diff --git a/networks/rhs.H b/networks/rhs.H
index 55a0060be2..e507d786aa 100644
--- a/networks/rhs.H
+++ b/networks/rhs.H
@@ -1348,9 +1348,10 @@ void rhs (burn_t& burn_state, Array1D<Real, 1, nrhs>& ydot)
 
     // Evaluate the neutrino cooling.
 #ifdef NEUTRINOS
+    constexpr int do_derivatives{0};
     Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
-    sneut5(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar,
-           sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
+    sneut5<do_derivatives>(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar,
+                           sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
 #else
     Real sneut = 0.0;
 #endif
@@ -1495,8 +1496,9 @@ void jac (burn_t& burn_state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac)
     // Evaluate the neutrino cooling.
 #ifdef NEUTRINOS
     Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
-    sneut5(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar,
-           sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar,
+                           sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
 #else
     Real sneut = 0.0, dsneutdt = 0.0, dsneutdd = 0.0, dsnuda = 0.0, dsnudz = 0.0;
     amrex::ignore_unused(sneut, dsneutdd);
diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H
index 6b066a2895..4a1b5aaf19 100644
--- a/neutrinos/sneut5.H
+++ b/neutrinos/sneut5.H
@@ -202,6 +202,8 @@ Real zfermim12(const Real x)
     return zfermim12r;
 }
 
+
+template <int do_derivatives>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void sneut5(const Real temp, const Real den,
             const Real abar, const Real zbar,
@@ -366,15 +368,18 @@ void sneut5(const Real temp, const Real den,
     rmda   = -rm*abari;
     rmdz   = den*abari;
     rmi    = 1.0e0_rt/rm;
-;
+
     a0     = rm * 1.0e-9_rt;
     a1     = std::pow(a0, oneth);
     zeta   = a1 * xlm1;
-    zetadt = -a1 * xlm2 * xldt;
-    a2     = oneth * a1*rmi * xlm1;
-    zetada = a2 * rmda;
-    zetadz = a2 * rmdz;
-;
+
+    if constexpr (do_derivatives) {
+        zetadt = -a1 * xlm2 * xldt;
+        a2     = oneth * a1*rmi * xlm1;
+        zetada = a2 * rmda;
+        zetadz = a2 * rmdz;
+    }
+
     zeta2 = zeta * zeta;
     zeta3 = zeta2 * zeta;
 
@@ -383,7 +388,9 @@ void sneut5(const Real temp, const Real den,
 
     // equation 2.8
     gl   = 1.0e0_rt - 13.04e0_rt*xl2 +133.5e0_rt*xl4 +1534.0e0_rt*xl6 +918.6e0_rt*xl8;
-    gldt = xldt*(-26.08e0_rt*xl +534.0e0_rt*xl3 +9204.0e0_rt*xl5 +7348.8e0_rt*xl7);
+    if constexpr (do_derivatives) {
+        gldt = xldt*(-26.08e0_rt*xl +534.0e0_rt*xl3 +9204.0e0_rt*xl5 +7348.8e0_rt*xl7);
+    }
 
     // equation 2.7
 
@@ -399,10 +406,12 @@ void sneut5(const Real temp, const Real den,
     }
 
     xnum   = a1 * b1;
-    c      = a2*b1 + a1*b2;
-    xnumdt = c*zetadt;
-    xnumda = c*zetada;
-    xnumdz = c*zetadz;
+    if constexpr (do_derivatives) {
+        c      = a2*b1 + a1*b2;
+        xnumdt = c*zetadt;
+        xnumda = c*zetada;
+        xnumdz = c*zetadz;
+    }
 
     if (t9 < 10.0_rt) {
        a1   = 9.383e-1_rt*xlm1 - 4.141e-1_rt*xlm2 + 5.829e-2_rt*xlm3;
@@ -412,18 +421,21 @@ void sneut5(const Real temp, const Real den,
        a2   = -1.2383e0_rt*xlm2 + 2.0e0_rt*8.141e-1_rt*xlm3;
     }
 
-    b1   = 3.0e0_rt*zeta2;
-
     xden   = zeta3 + a1;
-    xdendt = b1*zetadt + a2*xldt;
-    xdenda = b1*zetada;
-    xdendz = b1*zetadz;
+    if constexpr (do_derivatives) {
+        b1   = 3.0e0_rt*zeta2;
+        xdendt = b1*zetadt + a2*xldt;
+        xdenda = b1*zetada;
+        xdendz = b1*zetadz;
+    }
 
     a1      = 1.0e0_rt/xden;
     fpair   = xnum*a1;
-    fpairdt = (xnumdt - fpair*xdendt)*a1;
-    fpairda = (xnumda - fpair*xdenda)*a1;
-    fpairdz = (xnumdz - fpair*xdendz)*a1;
+    if constexpr (do_derivatives) {
+        fpairdt = (xnumdt - fpair*xdendt)*a1;
+        fpairda = (xnumda - fpair*xdenda)*a1;
+        fpairdz = (xnumdz - fpair*xdendz)*a1;
+    }
 
     // equation 2.6
     a1     = 10.7480e0_rt*xl2 + 0.3967e0_rt*xlp5 + 1.005e0_rt;
@@ -438,40 +450,49 @@ void sneut5(const Real temp, const Real den,
     b1     = 1.0e0_rt + rm*c;
 
     xden   = std::pow(b1, -0.3e0_rt);
-
-    d      = -0.3e0_rt*xden/b1;
-    xdendt = -d*rm*c*c*a2;
-    xdenda = d*rmda*c;
-    xdendz = d*rmdz*c;
+    if constexpr (do_derivatives) {
+        d      = -0.3e0_rt*xden/b1;
+        xdendt = -d*rm*c*c*a2;
+        xdenda = d*rmda*c;
+        xdendz = d*rmdz*c;
+    }
 
     qpair   = xnum*xden;
-    qpairdt = xnumdt*xden + xnum*xdendt;
-    qpairda = xnum*xdenda;
-    qpairdz = xnum*xdendz;
+    if constexpr (do_derivatives) {
+        qpairdt = xnumdt*xden + xnum*xdendt;
+        qpairda = xnum*xdenda;
+        qpairdz = xnum*xdendz;
+    }
 
     // equation 2.5
     a1    = std::exp(-2.0e0_rt*xlm1);
-    a2    = a1*2.0e0_rt*xlm2*xldt;
 
     spair   = a1*fpair;
-    spairdt = a2*fpair + a1*fpairdt;
-    spairda = a1*fpairda;
-    spairdz = a1*fpairdz;
+    if constexpr (do_derivatives) {
+        a2    = a1*2.0e0_rt*xlm2*xldt;
+        spairdt = a2*fpair + a1*fpairdt;
+        spairda = a1*fpairda;
+        spairdz = a1*fpairdz;
+    }
 
     a1      = spair;
     spair   = gl*a1;
-    spairdt = gl*spairdt + gldt*a1;
-    spairda = gl*spairda;
-    spairdz = gl*spairdz;
+    if constexpr (do_derivatives) {
+        spairdt = gl*spairdt + gldt*a1;
+        spairda = gl*spairda;
+        spairdz = gl*spairdz;
+    }
 
     a1      = tfac4*(1.0e0_rt + tfac3 * qpair);
     a2      = tfac4*tfac3;
 
     a3      = spair;
     spair   = a1*a3;
-    spairdt = a1*spairdt + a2*qpairdt*a3;
-    spairda = a1*spairda + a2*qpairda*a3;
-    spairdz = a1*spairdz + a2*qpairdz*a3;
+    if constexpr (do_derivatives) {
+        spairdt = a1*spairdt + a2*qpairdt*a3;
+        spairda = a1*spairda + a2*qpairda*a3;
+        spairdz = a1*spairdz + a2*qpairdz*a3;
+    }
 
     // plasma neutrino section
     // for collective reactions like gamma_plasmon => nu_e + nubar_e
@@ -488,10 +509,12 @@ void sneut5(const Real temp, const Real den,
 
     gl2   = 1.1095e11_rt * rm * c00;
 
-    gl2dt = -2.0e0_rt*gl2*tempi;
-    d     = rm*c00*b2*0.5e0_rt*b2*a3*1.019e-6_rt;
-    gl2da = 1.1095e11_rt * (rmda*c00  - d*rmda);
-    gl2dz = 1.1095e11_rt * (rmdz*c00  - d*rmdz);
+    if constexpr (do_derivatives) {
+        gl2dt = -2.0e0_rt*gl2*tempi;
+        d     = rm*c00*b2*0.5e0_rt*b2*a3*1.019e-6_rt;
+        gl2da = 1.1095e11_rt * (rmda*c00  - d*rmda);
+        gl2dz = 1.1095e11_rt * (rmdz*c00  - d*rmdz);
+    }
 
     gl    = std::sqrt(gl2);
     gl12  = std::sqrt(gl);
@@ -502,10 +525,12 @@ void sneut5(const Real temp, const Real den,
     // equation 4.7
     ft   = 2.4e0_rt + 0.6e0_rt*gl12 + 0.51e0_rt*gl + 1.25e0_rt*gl32;
     gum  = 1.0e0_rt/gl2;
-    a1   =(0.25e0_rt*0.6e0_rt*gl12 +0.5e0_rt*0.51e0_rt*gl +0.75e0_rt*1.25e0_rt*gl32)*gum;
-    ftdt = a1*gl2dt;
-    ftda = a1*gl2da;
-    ftdz = a1*gl2dz;
+    if constexpr (do_derivatives) {
+        a1   =(0.25e0_rt*0.6e0_rt*gl12 +0.5e0_rt*0.51e0_rt*gl +0.75e0_rt*1.25e0_rt*gl32)*gum;
+        ftdt = a1*gl2dt;
+        ftda = a1*gl2da;
+        ftdz = a1*gl2dz;
+    }
 
     // equation 4.8
     a1   = 8.6e0_rt*gl2 + 1.35e0_rt*gl72;
@@ -517,25 +542,29 @@ void sneut5(const Real temp, const Real den,
     c    = 1.0e0_rt/b1;
     fl   = a1*c;
 
-    d    = (a2 - fl*b2)*c;
-    fldt = d*gl2dt;
-    flda = d*gl2da;
-    fldz = d*gl2dz;
+    if constexpr (do_derivatives) {
+        d    = (a2 - fl*b2)*c;
+        fldt = d*gl2dt;
+        flda = d*gl2da;
+        fldz = d*gl2dz;
+    }
 
     // equation 4.9 and 4.10
     cc   = std::log10(2.0e0_rt*rm);
     xlnt = std::log10(temp);
 
     xnum   = sixth * (17.5e0_rt + cc - 3.0e0_rt*xlnt);
-    xnumdt = -iln10*0.5e0_rt*tempi;
-    a2     = iln10*sixth*rmi;
-    xnumda = a2*rmda;
-    xnumdz = a2*rmdz;
-
     xden   = sixth * (-24.5e0_rt + cc + 3.0e0_rt*xlnt);
-    xdendt = iln10*0.5e0_rt*tempi;
-    xdenda = a2*rmda;
-    xdendz = a2*rmdz;
+    if constexpr (do_derivatives) {
+        xnumdt = -iln10*0.5e0_rt*tempi;
+        a2     = iln10*sixth*rmi;
+        xnumda = a2*rmda;
+        xnumdz = a2*rmdz;
+
+        xdendt = iln10*0.5e0_rt*tempi;
+        xdenda = a2*rmda;
+        xdendz = a2*rmdz;
+    }
 
     // equation 4.11
     if (std::abs(xnum) > 0.7e0_rt || xden < 0.0e0_rt) {
@@ -554,14 +583,16 @@ void sneut5(const Real temp, const Real den,
        b2  = -b1*2.0e0_rt*(4.5e0_rt*xnum + 0.9e0_rt)*4.5e0_rt;
 
        c   = amrex::min(0.0e0_rt, xden - 1.6e0_rt + 1.25e0_rt*xnum);
-       if (c == 0.0_rt) {
-          dumdt = 0.0e0_rt;
-          dumda = 0.0e0_rt;
-          dumdz = 0.0e0_rt;
-       } else {
-          dumdt = xdendt + 1.25e0_rt*xnumdt;
-          dumda = xdenda + 1.25e0_rt*xnumda;
-          dumdz = xdendz + 1.25e0_rt*xnumdz;
+       if constexpr (do_derivatives) {
+           if (c == 0.0_rt) {
+               dumdt = 0.0e0_rt;
+               dumda = 0.0e0_rt;
+               dumdz = 0.0e0_rt;
+           } else {
+               dumdt = xdendt + 1.25e0_rt*xnumdt;
+               dumda = xdenda + 1.25e0_rt*xnumda;
+               dumdz = xdendz + 1.25e0_rt*xnumdz;
+           }
        }
 
        d   = 0.57e0_rt - 0.25e0_rt*xnum;
@@ -569,49 +600,58 @@ void sneut5(const Real temp, const Real den,
        c00 = std::exp(-a3*a3);
 
        f1  = -c00*2.0e0_rt*a3/d;
-       c01 = f1*(dumdt + a3*0.25e0_rt*xnumdt);
-       c03 = f1*(dumda + a3*0.25e0_rt*xnumda);
-       c04 = f1*(dumdz + a3*0.25e0_rt*xnumdz);
-
        fxy   = 1.05e0_rt + (a1 - b1)*c00;
-       fxydt = (a2*xnumdt -  b2*xnumdt)*c00 + (a1-b1)*c01;
-       fxyda = (a2*xnumda -  b2*xnumda)*c00 + (a1-b1)*c03;
-       fxydz = (a2*xnumdz -  b2*xnumdz)*c00 + (a1-b1)*c04;
-
+       if constexpr (do_derivatives) {
+           c01 = f1*(dumdt + a3*0.25e0_rt*xnumdt);
+           c03 = f1*(dumda + a3*0.25e0_rt*xnumda);
+           c04 = f1*(dumdz + a3*0.25e0_rt*xnumdz);
+
+           fxydt = (a2*xnumdt -  b2*xnumdt)*c00 + (a1-b1)*c01;
+           fxyda = (a2*xnumda -  b2*xnumda)*c00 + (a1-b1)*c03;
+           fxydz = (a2*xnumdz -  b2*xnumdz)*c00 + (a1-b1)*c04;
+       }
     }
 
     // equation 4.1 and 4.5
     splas   = (ft + fl) * fxy;
-    splasdt = (ftdt + fldt)*fxy + (ft+fl)*fxydt;
-    splasda = (ftda + flda)*fxy + (ft+fl)*fxyda;
-    splasdz = (ftdz + fldz)*fxy + (ft+fl)*fxydz;
+    if constexpr (do_derivatives) {
+        splasdt = (ftdt + fldt)*fxy + (ft+fl)*fxydt;
+        splasda = (ftda + flda)*fxy + (ft+fl)*fxyda;
+        splasdz = (ftdz + fldz)*fxy + (ft+fl)*fxydz;
+    }
 
     a2      = std::exp(-gl);
     a3      = -0.5e0_rt*a2*gl*gum;
 
     a1      = splas;
     splas   = a2*a1;
-    splasdt = a2*splasdt + a3*gl2dt*a1;
-    splasda = a2*splasda + a3*gl2da*a1;
-    splasdz = a2*splasdz + a3*gl2dz*a1;
+    if constexpr (do_derivatives) {
+        splasdt = a2*splasdt + a3*gl2dt*a1;
+        splasda = a2*splasda + a3*gl2da*a1;
+        splasdz = a2*splasdz + a3*gl2dz*a1;
+    }
 
     a2      = gl6;
     a3      = 3.0e0_rt*gl6*gum;
 
     a1      = splas;
     splas   = a2*a1;
-    splasdt = a2*splasdt + a3*gl2dt*a1;
-    splasda = a2*splasda + a3*gl2da*a1;
-    splasdz = a2*splasdz + a3*gl2dz*a1;
+    if constexpr (do_derivatives) {
+        splasdt = a2*splasdt + a3*gl2dt*a1;
+        splasda = a2*splasda + a3*gl2da*a1;
+        splasdz = a2*splasdz + a3*gl2dz*a1;
+    }
 
     a2      = 0.93153e0_rt * 3.0e21_rt * xl9;
     a3      = 0.93153e0_rt * 3.0e21_rt * 9.0e0_rt*xl8*xldt;
 
     a1      = splas;
     splas   = a2*a1;
-    splasdt = a2*splasdt + a3*a1;
-    splasda = a2*splasda;
-    splasdz = a2*splasdz;
+    if constexpr (do_derivatives) {
+        splasdt = a2*splasdt + a3*a1;
+        splasda = a2*splasda;
+        splasdz = a2*splasdz;
+    }
 
     // photoneutrino process section
     // for reactions like e- + gamma => e- + nu_e + nubar_e
@@ -797,79 +837,104 @@ void sneut5(const Real temp, const Real den,
 
     // equation 3.4
     dum   = a0 + a1*zeta + a2*zeta2;
-    dumdt = f0 + f1*zeta + a1*zetadt + f2*zeta2 + 2.0e0_rt*a2*zeta*zetadt;
-    dumda = a1*zetada + 2.0e0_rt*a2*zeta*zetada;
-    dumdz = a1*zetadz + 2.0e0_rt*a2*zeta*zetadz;
+    if constexpr (do_derivatives) {
+        dumdt = f0 + f1*zeta + a1*zetadt + f2*zeta2 + 2.0e0_rt*a2*zeta*zetadt;
+        dumda = a1*zetada + 2.0e0_rt*a2*zeta*zetada;
+        dumdz = a1*zetadz + 2.0e0_rt*a2*zeta*zetadz;
+    }
 
     z      = std::exp(-cc*zeta);
 
     xnum   = dum*z;
-    xnumdt = dumdt*z - dum*z*cc*zetadt;
-    xnumda = dumda*z - dum*z*cc*zetada;
-    xnumdz = dumdz*z - dum*z*cc*zetadz;
+    if constexpr (do_derivatives) {
+        xnumdt = dumdt*z - dum*z*cc*zetadt;
+        xnumda = dumda*z - dum*z*cc*zetada;
+        xnumdz = dumdz*z - dum*z*cc*zetadz;
+    }
 
     xden   = zeta3 + 6.290e-3_rt*xlm1 + 7.483e-3_rt*xlm2 + 3.061e-4_rt*xlm3;
 
     dum    = 3.0e0_rt*zeta2;
-    xdendt = dum*zetadt - xldt*(6.290e-3_rt*xlm2
-         + 2.0e0_rt*7.483e-3_rt*xlm3 + 3.0e0_rt*3.061e-4_rt*xlm4);
-    xdenda = dum*zetada;
-    xdendz = dum*zetadz;
-
+    if constexpr (do_derivatives) {
+        xdendt = dum*zetadt - xldt*(6.290e-3_rt*xlm2
+                                    + 2.0e0_rt*7.483e-3_rt*xlm3 + 3.0e0_rt*3.061e-4_rt*xlm4);
+        xdenda = dum*zetada;
+        xdendz = dum*zetadz;
+    }
     dum      = 1.0e0_rt/xden;
     fphot   = xnum*dum;
-    fphotdt = (xnumdt - fphot*xdendt)*dum;
-    fphotda = (xnumda - fphot*xdenda)*dum;
-    fphotdz = (xnumdz - fphot*xdendz)*dum;
+    if constexpr (do_derivatives) {
+        fphotdt = (xnumdt - fphot*xdendt)*dum;
+        fphotda = (xnumda - fphot*xdenda)*dum;
+        fphotdz = (xnumdz - fphot*xdendz)*dum;
+    }
 
     // equation 3.3
     a0     = 1.0e0_rt + 2.045e0_rt * xl;
     xnum   = 0.666e0_rt*std::pow(a0, -2.066e0_rt);
-    xnumdt = -2.066e0_rt*xnum/a0 * 2.045e0_rt*xldt;
+    if constexpr (do_derivatives) {
+        xnumdt = -2.066e0_rt*xnum/a0 * 2.045e0_rt*xldt;
+    }
 
     dum    = 1.875e8_rt*xl + 1.653e8_rt*xl2 + 8.499e8_rt*xl3 - 1.604e8_rt*xl4;
-    dumdt  = xldt*(1.875e8_rt + 2.0e0_rt*1.653e8_rt*xl + 3.0e0_rt*8.499e8_rt*xl2
-             - 4.0e0_rt*1.604e8_rt*xl3);
+    if constexpr (do_derivatives) {
+        dumdt  = xldt*(1.875e8_rt + 2.0e0_rt*1.653e8_rt*xl + 3.0e0_rt*8.499e8_rt*xl2
+                       - 4.0e0_rt*1.604e8_rt*xl3);
+    }
 
     z      = 1.0e0_rt/dum;
     xden   = 1.0e0_rt + rm*z;
-    xdendt =  -rm*z*z*dumdt;
-    xdenda =  rmda*z;
-    xdendz =  rmdz*z;
+    if constexpr (do_derivatives) {
+        xdendt =  -rm*z*z*dumdt;
+        xdenda =  rmda*z;
+        xdendz =  rmdz*z;
+    }
 
     z      = 1.0e0_rt/xden;
     qphot = xnum*z;
-    qphotdt = (xnumdt - qphot*xdendt)*z;
+    if constexpr (do_derivatives) {
+        qphotdt = (xnumdt - qphot*xdendt)*z;
+    }
     dum      = -qphot*z;
-    qphotda = dum*xdenda;
-    qphotdz = dum*xdendz;
+    if constexpr (do_derivatives) {
+        qphotda = dum*xdenda;
+        qphotdz = dum*xdendz;
+    }
 
     // equation 3.2
     sphot   = xl5 * fphot;
-    sphotdt = 5.0e0_rt*xl4*xldt*fphot + xl5*fphotdt;
-    sphotda = xl5*fphotda;
-    sphotdz = xl5*fphotdz;
+    if constexpr (do_derivatives) {
+        sphotdt = 5.0e0_rt*xl4*xldt*fphot + xl5*fphotdt;
+        sphotda = xl5*fphotda;
+        sphotdz = xl5*fphotdz;
+    }
 
     a1      = sphot;
     sphot   = rm*a1;
-    sphotdt = rm*sphotdt;
-    sphotda = rm*sphotda + rmda*a1;
-    sphotdz = rm*sphotdz + rmdz*a1;
+    if constexpr (do_derivatives) {
+        sphotdt = rm*sphotdt;
+        sphotda = rm*sphotda + rmda*a1;
+        sphotdz = rm*sphotdz + rmdz*a1;
+    }
 
     a1      = tfac4*(1.0e0_rt - tfac3 * qphot);
     a2      = -tfac4*tfac3;
 
     a3      = sphot;
     sphot   = a1*a3;
-    sphotdt = a1*sphotdt + a2*qphotdt*a3;
-    sphotda = a1*sphotda + a2*qphotda*a3;
-    sphotdz = a1*sphotdz + a2*qphotdz*a3;
+    if constexpr (do_derivatives) {
+        sphotdt = a1*sphotdt + a2*qphotdt*a3;
+        sphotda = a1*sphotda + a2*qphotda*a3;
+        sphotdz = a1*sphotdz + a2*qphotdz*a3;
+    }
 
     if (sphot <= 0.0_rt) {
        sphot   = 0.0e0_rt;
-       sphotdt = 0.0e0_rt;
-       sphotda = 0.0e0_rt;
-       sphotdz = 0.0e0_rt;
+       if constexpr (do_derivatives) {
+           sphotdt = 0.0e0_rt;
+           sphotda = 0.0e0_rt;
+           sphotdz = 0.0e0_rt;
+       }
     }
 
     // bremsstrahlung neutrino section
@@ -900,13 +965,17 @@ void sneut5(const Real temp, const Real den,
 
        // equation 5.3
        dum   = 7.05e6_rt * t832 + 5.12e4_rt * t83;
-       dumdt = (1.5e0_rt*7.05e6_rt*t812 + 3.0e0_rt*5.12e4_rt*t82)*1.0e-8_rt;
+       if constexpr (do_derivatives) {
+           dumdt = (1.5e0_rt*7.05e6_rt*t812 + 3.0e0_rt*5.12e4_rt*t82)*1.0e-8_rt;
+       }
 
        z     = 1.0e0_rt/dum;
        eta   = rm*z;
-       etadt = -rm*z*z*dumdt;
-       etada = rmda*z;
-       etadz = rmdz*z;
+       if constexpr (do_derivatives) {
+           etadt = -rm*z*z*dumdt;
+           etada = rmda*z;
+           etadz = rmdz*z;
+       }
 
        etam1 = 1.0e0_rt/eta;
        etam2 = etam1 * etam1;
@@ -918,27 +987,35 @@ void sneut5(const Real temp, const Real den,
        xnum  = 1.0e0_rt/a0;
 
        dum   = 1.0e0_rt + 1.47e0_rt*etam1 + 3.29e-2_rt*etam2;
-       z     = -1.47e0_rt*etam2 - 2.0e0_rt*3.29e-2_rt*etam3;
-       dumdt = z*etadt;
-       dumda = z*etada;
-       dumdz = z*etadz;
+       if constexpr (do_derivatives) {
+           z     = -1.47e0_rt*etam2 - 2.0e0_rt*3.29e-2_rt*etam3;
+           dumdt = z*etadt;
+           dumda = z*etada;
+           dumdz = z*etadz;
+       }
 
        c00   = 1.26e0_rt * (1.0e0_rt+etam1);
-       z     = -1.26e0_rt*etam2;
-       c01   = z*etadt;
-       c03   = z*etada;
-       c04   = z*etadz;
+       if constexpr (do_derivatives) {
+           z     = -1.26e0_rt*etam2;
+           c01   = z*etadt;
+           c03   = z*etada;
+           c04   = z*etadz;
+       }
 
        z      = 1.0e0_rt/dum;
        xden   = c00*z;
-       xdendt = (c01 - xden*dumdt)*z;
-       xdenda = (c03 - xden*dumda)*z;
-       xdendz = (c04 - xden*dumdz)*z;
+       if constexpr (do_derivatives) {
+           xdendt = (c01 - xden*dumdt)*z;
+           xdenda = (c03 - xden*dumda)*z;
+           xdendz = (c04 - xden*dumdz)*z;
+       }
 
        fbrem   = xnum + xden;
-       fbremdt = -xnum*xnum*f0 + xdendt;
-       fbremda = xdenda;
-       fbremdz = xdendz;
+       if constexpr (do_derivatives) {
+           fbremdt = -xnum*xnum*f0 + xdendt;
+           fbremda = xdenda;
+           fbremdz = xdendz;
+       }
 
        // equation 5.9
        a0    = 230.0e0_rt + 6.7e5_rt*t8m2 + 7.66e9_rt*t8m5;
@@ -946,16 +1023,20 @@ void sneut5(const Real temp, const Real den,
 
        z     = 1.0e0_rt + rm*1.0e-9_rt;
        dum   = a0*z;
-       dumdt = f0*z;
-       z     = a0*1.0e-9_rt;
-       dumda = z*rmda;
-       dumdz = z*rmdz;
+       if constexpr (do_derivatives) {
+           dumdt = f0*z;
+           z     = a0*1.0e-9_rt;
+           dumda = z*rmda;
+           dumdz = z*rmdz;
+       }
 
        xnum   = 1.0e0_rt/dum;
-       z      = -xnum*xnum;
-       xnumdt = z*dumdt;
-       xnumda = z*dumda;
-       xnumdz = z*dumdz;
+       if constexpr (do_derivatives) {
+           z      = -xnum*xnum;
+           xnumdt = z*dumdt;
+           xnumda = z*dumda;
+           xnumdz = z*dumdz;
+       }
 
        c00   = 7.75e5_rt*t832 + 247.0e0_rt * std::pow(t8, 3.85e0_rt);
        dd00  = (1.5e0_rt*7.75e5_rt*t812 + 3.85e0_rt*247.0e0_rt*std::pow(t8, 2.85e0_rt))*1.0e-8_rt;
@@ -968,33 +1049,43 @@ void sneut5(const Real temp, const Real den,
 
        z     = std::pow(den, 0.656e0_rt);
        dum   = c00*rmi  + c01  + c02*z;
-       dumdt = dd00*rmi + dd01 + dd02*z;
-       z     = -c00*rmi*rmi;
-       dumda = z*rmda;
-       dumdz = z*rmdz;
+       if constexpr (do_derivatives) {
+           dumdt = dd00*rmi + dd01 + dd02*z;
+           z     = -c00*rmi*rmi;
+           dumda = z*rmda;
+           dumdz = z*rmdz;
+       }
 
        xden  = 1.0e0_rt/dum;
-       z      = -xden*xden;
-       xdendt = z*dumdt;
-       xdenda = z*dumda;
-       xdendz = z*dumdz;
+       if constexpr (do_derivatives) {
+           z      = -xden*xden;
+           xdendt = z*dumdt;
+           xdenda = z*dumda;
+           xdendz = z*dumdz;
+       }
 
        gbrem   = xnum + xden;
-       gbremdt = xnumdt + xdendt;
-       gbremda = xnumda + xdenda;
-       gbremdz = xnumdz + xdendz;
+       if constexpr (do_derivatives) {
+           gbremdt = xnumdt + xdendt;
+           gbremda = xnumda + xdenda;
+           gbremdz = xnumdz + xdendz;
+       }
 
        // equation 5.1
        dum    = 0.5738e0_rt*zbar*ye*t86*den;
-       dumdt  = 0.5738e0_rt*zbar*ye*6.0e0_rt*t85*den*1.0e-8_rt;
-       dumda  = -dum*abari;
-       dumdz  = 0.5738e0_rt*2.0e0_rt*ye*t86*den;
+       if constexpr (do_derivatives) {
+           dumdt  = 0.5738e0_rt*zbar*ye*6.0e0_rt*t85*den*1.0e-8_rt;
+           dumda  = -dum*abari;
+           dumdz  = 0.5738e0_rt*2.0e0_rt*ye*t86*den;
+       }
 
        z       = tfac4*fbrem - tfac5*gbrem;
        sbrem   = dum * z;
-       sbremdt = dumdt*z + dum*(tfac4*fbremdt - tfac5*gbremdt);
-       sbremda = dumda*z + dum*(tfac4*fbremda - tfac5*gbremda);
-       sbremdz = dumdz*z + dum*(tfac4*fbremdz - tfac5*gbremdz);
+       if constexpr (do_derivatives) {
+           sbremdt = dumdt*z + dum*(tfac4*fbremdt - tfac5*gbremdt);
+           sbremda = dumda*z + dum*(tfac4*fbremda - tfac5*gbremda);
+           sbremdz = dumdz*z + dum*(tfac4*fbremdz - tfac5*gbremdz);
+       }
 
        // liquid metal with c12 parameters (not too different for other elements)
        // equation 5.18 and 5.16
@@ -1078,9 +1169,11 @@ void sneut5(const Real temp, const Real den,
        //      - 0.00069e0_rt*sin5*5.0e0_rt);
 
        dum   = 2.275e-1_rt * zbar * zbar*t8m1 * std::pow(den6*abari, oneth);
-       dumdt = -dum*tempi;
-       dumda = -oneth*dum*abari;
-       dumdz = 2.0e0_rt*dum*zbari;
+       if constexpr (do_derivatives) {
+           dumdt = -dum*tempi;
+           dumda = -oneth*dum*abari;
+           dumdz = 2.0e0_rt*dum*zbari;
+       }
 
        gm1   = 1.0e0_rt/dum;
        gm2   = gm1*gm1;
@@ -1098,46 +1191,56 @@ void sneut5(const Real temp, const Real den,
 
        // equation 5.19 and 5.20
        fliq   = v*fb + (1.0e0_rt - v)*ft;
-       fliqdt = a0*dumdt*(fb - ft);
-       fliqda = a0*dumda*(fb - ft);
-       fliqdz = a0*dumdz*(fb - ft);
+       if constexpr (do_derivatives) {
+           fliqdt = a0*dumdt*(fb - ft);
+           fliqda = a0*dumda*(fb - ft);
+           fliqdz = a0*dumdz*(fb - ft);
+       }
 
        gliq   = w*gb + (1.0e0_rt - w)*gt;
-       gliqdt = a1*dumdt*(gb - gt);
-       gliqda = a1*dumda*(gb - gt);
-       gliqdz = a1*dumdz*(gb - gt);
+       if constexpr (do_derivatives) {
+           gliqdt = a1*dumdt*(gb - gt);
+           gliqda = a1*dumda*(gb - gt);
+           gliqdz = a1*dumdz*(gb - gt);
+       }
 
        // equation 5.17
        dum    = 0.5738e0_rt*zbar*ye*t86*den;
-       dumdt  = 0.5738e0_rt*zbar*ye*6.0e0_rt*t85*den*1.0e-8_rt;
-       dumda  = -dum*abari;
-       dumdz  = 0.5738e0_rt*2.0e0_rt*ye*t86*den;
+       if constexpr (do_derivatives) {
+           dumdt  = 0.5738e0_rt*zbar*ye*6.0e0_rt*t85*den*1.0e-8_rt;
+           dumda  = -dum*abari;
+           dumdz  = 0.5738e0_rt*2.0e0_rt*ye*t86*den;
+       }
 
        z       = tfac4*fliq - tfac5*gliq;
        sbrem   = dum * z;
-       sbremdt = dumdt*z + dum*(tfac4*fliqdt - tfac5*gliqdt);
-       sbremda = dumda*z + dum*(tfac4*fliqda - tfac5*gliqda);
-       sbremdz = dumdz*z + dum*(tfac4*fliqdz - tfac5*gliqdz);
-
+       if constexpr (do_derivatives) {
+           sbremdt = dumdt*z + dum*(tfac4*fliqdt - tfac5*gliqdt);
+           sbremda = dumda*z + dum*(tfac4*fliqda - tfac5*gliqda);
+           sbremdz = dumdz*z + dum*(tfac4*fliqdz - tfac5*gliqdz);
+       }
     }
 
     // recombination neutrino section
     // for reactions like e- (continuum) => e- (bound) + nu_e + nubar_e
     // equation 6.11 solved for nu
     xnum   = 1.10520e8_rt * den * ye /(temp*std::sqrt(temp));
-    xnumdt = -1.50e0_rt*xnum*tempi;
-    xnumda = -xnum*abari;
-    xnumdz = xnum*zbari;
+    if constexpr (do_derivatives) {
+        xnumdt = -1.50e0_rt*xnum*tempi;
+        xnumda = -xnum*abari;
+        xnumdz = xnum*zbari;
+    }
 
     // the chemical potential
     nu   = ifermi12(xnum);
 
     // a0 is d(nu)/d(xnum)
     a0 = 1.0e0_rt/(0.5e0_rt*zfermim12(nu));
-    nudt = a0*xnumdt;
-    nuda = a0*xnumda;
-    nudz = a0*xnumdz;
-
+    if constexpr (do_derivatives) {
+        nudt = a0*xnumdt;
+        nuda = a0*xnumda;
+        nudz = a0*xnumdz;
+    }
     nu2  = nu * nu;
     nu3  = nu2 * nu;
 
@@ -1168,16 +1271,20 @@ void sneut5(const Real temp, const Real den,
     if (nu >= -20.0_rt &&  nu <= 10.0_rt) {
 
        zeta   = 1.579e5_rt*zbar*zbar*tempi;
-       zetadt = -zeta*tempi;
-       zetada = 0.0e0_rt;
-       zetadz = 2.0e0_rt*zeta*zbari;
+       if constexpr (do_derivatives) {
+           zetadt = -zeta*tempi;
+           zetada = 0.0e0_rt;
+           zetadz = 2.0e0_rt*zeta*zbari;
+       }
 
        c00    = 1.0e0_rt/(1.0e0_rt + f1*nu + f2*nu2 + f3*nu3);
        c01    = f1 + f2*2.0e0_rt*nu + f3*3.0e0_rt*nu2;
        dum    = zeta*c00;
-       dumdt  = zetadt*c00 + zeta*c01*nudt;
-       dumda  = zeta*c01*nuda;
-       dumdz  = zetadz*c00 + zeta*c01*nudz;
+       if constexpr (do_derivatives) {
+           dumdt  = zetadt*c00 + zeta*c01*nudt;
+           dumda  = zeta*c01*nuda;
+           dumdz  = zetadz*c00 + zeta*c01*nudz;
+       }
 
        z      = 1.0e0_rt/dum;
        dd00   = std::pow(dum, -2.25_rt);
@@ -1188,17 +1295,21 @@ void sneut5(const Real temp, const Real den,
        z      = std::exp(c*nu);
        dd00   = b*z*(1.0e0_rt + d*dum);
        gum    = 1.0e0_rt + dd00;
-       gumdt  = dd00*c*nudt + b*z*d*dumdt;
-       gumda  = dd00*c*nuda + b*z*d*dumda;
-       gumdz  = dd00*c*nudz + b*z*d*dumdz;
+       if constexpr (do_derivatives) {
+           gumdt  = dd00*c*nudt + b*z*d*dumdt;
+           gumda  = dd00*c*nuda + b*z*d*dumda;
+           gumdz  = dd00*c*nudz + b*z*d*dumdz;
+       }
 
        z   = std::exp(nu);
        a1  = 1.0e0_rt/gum;
 
        bigj   = c00 * z * a1;
-       bigjdt = c01*dumdt*z*a1 + c00*z*nudt*a1 - c00*z*a1*a1 * gumdt;
-       bigjda = c01*dumda*z*a1 + c00*z*nuda*a1 - c00*z*a1*a1 * gumda;
-       bigjdz = c01*dumdz*z*a1 + c00*z*nudz*a1 - c00*z*a1*a1 * gumdz;
+       if constexpr (do_derivatives) {
+           bigjdt = c01*dumdt*z*a1 + c00*z*nudt*a1 - c00*z*a1*a1 * gumdt;
+           bigjda = c01*dumda*z*a1 + c00*z*nuda*a1 - c00*z*a1*a1 * gumda;
+           bigjdz = c01*dumdz*z*a1 + c00*z*nudz*a1 - c00*z*a1*a1 * gumdz;
+       }
 
        // equation 6.5
        z     = std::exp(zeta + nu);
@@ -1207,45 +1318,58 @@ void sneut5(const Real temp, const Real den,
        a2    = 1.0e0_rt/bigj;
 
        sreco   = tfac6 * 2.649e-18_rt * ye * std::pow(zbar, 13.0_rt) * den * bigj*a1;
-       srecodt = sreco*(bigjdt*a2 - z*(zetadt + nudt)*a1);
-       srecoda = sreco*(-1.0e0_rt*abari + bigjda*a2 - z*(zetada+nuda)*a1);
-       srecodz = sreco*(14.0e0_rt*zbari + bigjdz*a2 - z*(zetadz+nudz)*a1);
-
+       if constexpr (do_derivatives) {
+           srecodt = sreco*(bigjdt*a2 - z*(zetadt + nudt)*a1);
+           srecoda = sreco*(-1.0e0_rt*abari + bigjda*a2 - z*(zetada+nuda)*a1);
+           srecodz = sreco*(14.0e0_rt*zbari + bigjdz*a2 - z*(zetadz+nudz)*a1);
+       }
     }
 
     // convert from erg/cm^3/s to erg/g/s
     // comment these out to duplicate the itoh et al plots
 
     spair   = spair*deni;
-    spairdt = spairdt*deni;
-    spairda = spairda*deni;
-    spairdz = spairdz*deni;
+    if constexpr (do_derivatives) {
+        spairdt = spairdt*deni;
+        spairda = spairda*deni;
+        spairdz = spairdz*deni;
+    }
 
     splas   = splas*deni;
-    splasdt = splasdt*deni;
-    splasda = splasda*deni;
-    splasdz = splasdz*deni;
+    if constexpr (do_derivatives) {
+        splasdt = splasdt*deni;
+        splasda = splasda*deni;
+        splasdz = splasdz*deni;
+    }
 
     sphot   = sphot*deni;
-    sphotdt = sphotdt*deni;
-    sphotda = sphotda*deni;
-    sphotdz = sphotdz*deni;
+    if constexpr (do_derivatives) {
+        sphotdt = sphotdt*deni;
+        sphotda = sphotda*deni;
+        sphotdz = sphotdz*deni;
+    }
 
     sbrem   = sbrem*deni;
-    sbremdt = sbremdt*deni;
-    sbremda = sbremda*deni;
-    sbremdz = sbremdz*deni;
+    if constexpr (do_derivatives) {
+        sbremdt = sbremdt*deni;
+        sbremda = sbremda*deni;
+        sbremdz = sbremdz*deni;
+    }
 
     sreco   = sreco*deni;
-    srecodt = srecodt*deni;
-    srecoda = srecoda*deni;
-    srecodz = srecodz*deni;
+    if constexpr (do_derivatives) {
+        srecodt = srecodt*deni;
+        srecoda = srecoda*deni;
+        srecodz = srecodz*deni;
+    }
 
     // the total neutrino loss rate
     snu    =  splas + spair + sphot + sbrem + sreco;
-    dsnudt =  splasdt + spairdt + sphotdt + sbremdt + srecodt;
-    dsnuda =  splasda + spairda + sphotda + sbremda + srecoda;
-    dsnudz =  splasdz + spairdz + sphotdz + sbremdz + srecodz;
+    if constexpr (do_derivatives) {
+        dsnudt =  splasdt + spairdt + sphotdt + sbremdt + srecodt;
+        dsnuda =  splasda + spairda + sphotda + sbremda + srecoda;
+        dsnudz =  splasdz + spairdz + sphotdz + sbremdz + srecodz;
+    }
 
 }
 

From d41e92078c87200fd4fcdfbf7a049b1ee97ca988 Mon Sep 17 00:00:00 2001
From: Michael Zingale <michael.zingale@stonybrook.edu>
Date: Sun, 10 Sep 2023 14:52:40 -0400
Subject: [PATCH 2/5] more conversion

---
 networks/CNO_extras/actual_rhs.H              |  7 +--
 networks/ECSN/actual_rhs.H                    |  7 +--
 networks/ase/actual_rhs.H                     |  7 +--
 .../C-burn-simple/actual_rhs.H                |  7 +--
 .../ignition_reaclib/URCA-simple/actual_rhs.H |  7 +--
 networks/nova/actual_rhs.H                    |  7 +--
 networks/nova2/actual_rhs.H                   |  7 +--
 networks/partition_test/actual_rhs.H          |  7 +--
 networks/sn160/actual_rhs.H                   |  7 +--
 networks/subch_approx/actual_rhs.H            |  7 +--
 networks/subch_full/actual_rhs.H              |  7 +--
 networks/subch_simple/actual_rhs.H            |  7 +--
 neutrinos/sneut5.H                            | 45 ++++++++++---------
 13 files changed, 73 insertions(+), 56 deletions(-)

diff --git a/networks/CNO_extras/actual_rhs.H b/networks/CNO_extras/actual_rhs.H
index 79e1367582..03501b4031 100644
--- a/networks/CNO_extras/actual_rhs.H
+++ b/networks/CNO_extras/actual_rhs.H
@@ -794,8 +794,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -1197,7 +1197,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);
diff --git a/networks/ECSN/actual_rhs.H b/networks/ECSN/actual_rhs.H
index b7aa913121..0fb4bc3422 100644
--- a/networks/ECSN/actual_rhs.H
+++ b/networks/ECSN/actual_rhs.H
@@ -412,8 +412,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -629,7 +629,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * dsnuda + (zion[j-1] - state.zbar) * state.abar * dsnudz);
diff --git a/networks/ase/actual_rhs.H b/networks/ase/actual_rhs.H
index 26445601af..242edac898 100644
--- a/networks/ase/actual_rhs.H
+++ b/networks/ase/actual_rhs.H
@@ -1093,8 +1093,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -1646,7 +1646,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);
diff --git a/networks/ignition_reaclib/C-burn-simple/actual_rhs.H b/networks/ignition_reaclib/C-burn-simple/actual_rhs.H
index 6a84156b9a..5422a92f35 100644
--- a/networks/ignition_reaclib/C-burn-simple/actual_rhs.H
+++ b/networks/ignition_reaclib/C-burn-simple/actual_rhs.H
@@ -212,8 +212,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -309,7 +309,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);
diff --git a/networks/ignition_reaclib/URCA-simple/actual_rhs.H b/networks/ignition_reaclib/URCA-simple/actual_rhs.H
index cc76ef7078..be1b24da39 100644
--- a/networks/ignition_reaclib/URCA-simple/actual_rhs.H
+++ b/networks/ignition_reaclib/URCA-simple/actual_rhs.H
@@ -234,8 +234,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -343,7 +343,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);
diff --git a/networks/nova/actual_rhs.H b/networks/nova/actual_rhs.H
index a9e79d1fcc..260a47e241 100644
--- a/networks/nova/actual_rhs.H
+++ b/networks/nova/actual_rhs.H
@@ -487,8 +487,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -746,7 +746,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);
diff --git a/networks/nova2/actual_rhs.H b/networks/nova2/actual_rhs.H
index 49ea50923f..9cae820c4b 100644
--- a/networks/nova2/actual_rhs.H
+++ b/networks/nova2/actual_rhs.H
@@ -727,8 +727,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -1058,7 +1058,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);
diff --git a/networks/partition_test/actual_rhs.H b/networks/partition_test/actual_rhs.H
index 37d0cac6c0..cbff936ed8 100644
--- a/networks/partition_test/actual_rhs.H
+++ b/networks/partition_test/actual_rhs.H
@@ -202,8 +202,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -335,7 +335,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);
diff --git a/networks/sn160/actual_rhs.H b/networks/sn160/actual_rhs.H
index 8f81f200bb..647a980a70 100644
--- a/networks/sn160/actual_rhs.H
+++ b/networks/sn160/actual_rhs.H
@@ -10396,8 +10396,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -18028,7 +18028,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);
diff --git a/networks/subch_approx/actual_rhs.H b/networks/subch_approx/actual_rhs.H
index a9b90e7e67..c4e168ec2c 100644
--- a/networks/subch_approx/actual_rhs.H
+++ b/networks/subch_approx/actual_rhs.H
@@ -1375,8 +1375,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -1928,7 +1928,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);
diff --git a/networks/subch_full/actual_rhs.H b/networks/subch_full/actual_rhs.H
index d082bd48de..c27ce0b3d7 100644
--- a/networks/subch_full/actual_rhs.H
+++ b/networks/subch_full/actual_rhs.H
@@ -1603,8 +1603,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -2390,7 +2390,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);
diff --git a/networks/subch_simple/actual_rhs.H b/networks/subch_simple/actual_rhs.H
index 68e84fa45f..eae7f7bbbb 100644
--- a/networks/subch_simple/actual_rhs.H
+++ b/networks/subch_simple/actual_rhs.H
@@ -1227,8 +1227,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -1723,7 +1723,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);
diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H
index 4a1b5aaf19..81449010fb 100644
--- a/neutrinos/sneut5.H
+++ b/neutrinos/sneut5.H
@@ -808,32 +808,35 @@ void sneut5(const Real temp, const Real den,
          + c01*cos1 + dd01*sin1 + c02*cos2 + dd02*sin2
          + c03*cos3 + dd03*sin3 + c04*cos4 + dd04*sin4
          + c05*cos5 + dd05*sin5 + 0.5e0_rt*c06*last;
-
-    f0 =  taudt*fac1*(-c01*sin1 + dd01*cos1 - c02*sin2*2.0e0_rt
-         + dd02*cos2*2.0e0_rt - c03*sin3*3.0e0_rt + dd03*cos3*3.0e0_rt
-         - c04*sin4*4.0e0_rt + dd04*cos4*4.0e0_rt
-         - c05*sin5*5.0e0_rt + dd05*cos5*5.0e0_rt)
-         - 0.5e0_rt*c06*xast*fac2*taudt;
-
     a1 = 0.5e0_rt*c10
          + c11*cos1 + dd11*sin1 + c12*cos2 + dd12*sin2
          + c13*cos3 + dd13*sin3 + c14*cos4 + dd14*sin4
          + c15*cos5 + dd15*sin5 + 0.5e0_rt*c16*last;
 
-    f1 = taudt*fac1*(-c11*sin1 + dd11*cos1 - c12*sin2*2.0e0_rt
-         + dd12*cos2*2.0e0_rt - c13*sin3*3.0e0_rt + dd13*cos3*3.0e0_rt
-         - c14*sin4*4.0e0_rt + dd14*cos4*4.0e0_rt - c15*sin5*5.0e0_rt
-         + dd15*cos5*5.0e0_rt) - 0.5e0_rt*c16*xast*fac2*taudt;
 
     a2 = 0.5e0_rt*c20
          + c21*cos1 + dd21*sin1 + c22*cos2 + dd22*sin2
          + c23*cos3 + dd23*sin3 + c24*cos4 + dd24*sin4
          + c25*cos5 + dd25*sin5 + 0.5e0_rt*c26*last;
 
-    f2 = taudt*fac1*(-c21*sin1 + dd21*cos1 - c22*sin2*2.0e0_rt
-         + dd22*cos2*2.0e0_rt - c23*sin3*3.0e0_rt + dd23*cos3*3.0e0_rt
-         - c24*sin4*4.0e0_rt + dd24*cos4*4.0e0_rt - c25*sin5*5.0e0_rt
-         + dd25*cos5*5.0e0_rt) - 0.5e0_rt*c26*xast*fac2*taudt;
+    if constexpr (do_derivatives) {
+        f0 =  taudt*fac1*(-c01*sin1 + dd01*cos1 - c02*sin2*2.0e0_rt
+                          + dd02*cos2*2.0e0_rt - c03*sin3*3.0e0_rt + dd03*cos3*3.0e0_rt
+                          - c04*sin4*4.0e0_rt + dd04*cos4*4.0e0_rt
+                          - c05*sin5*5.0e0_rt + dd05*cos5*5.0e0_rt)
+            - 0.5e0_rt*c06*xast*fac2*taudt;
+
+        f1 = taudt*fac1*(-c11*sin1 + dd11*cos1 - c12*sin2*2.0e0_rt
+                         + dd12*cos2*2.0e0_rt - c13*sin3*3.0e0_rt + dd13*cos3*3.0e0_rt
+                         - c14*sin4*4.0e0_rt + dd14*cos4*4.0e0_rt - c15*sin5*5.0e0_rt
+                         + dd15*cos5*5.0e0_rt) - 0.5e0_rt*c16*xast*fac2*taudt;
+
+        f2 = taudt*fac1*(-c21*sin1 + dd21*cos1 - c22*sin2*2.0e0_rt
+                         + dd22*cos2*2.0e0_rt - c23*sin3*3.0e0_rt + dd23*cos3*3.0e0_rt
+                         - c24*sin4*4.0e0_rt + dd24*cos4*4.0e0_rt - c25*sin5*5.0e0_rt
+                         + dd25*cos5*5.0e0_rt) - 0.5e0_rt*c26*xast*fac2*taudt;
+    }
+
 
     // equation 3.4
     dum   = a0 + a1*zeta + a2*zeta2;
@@ -1039,13 +1042,15 @@ void sneut5(const Real temp, const Real den,
        }
 
        c00   = 7.75e5_rt*t832 + 247.0e0_rt * std::pow(t8, 3.85e0_rt);
-       dd00  = (1.5e0_rt*7.75e5_rt*t812 + 3.85e0_rt*247.0e0_rt*std::pow(t8, 2.85e0_rt))*1.0e-8_rt;
-
        c01   = 4.07e0_rt + 0.0240e0_rt * std::pow(t8, 1.4e0_rt);
-       dd01  = 1.4e0_rt*0.0240e0_rt * std::pow(t8, 0.4e0_rt)*1.0e-8_rt;
-
        c02   = 4.59e-5_rt * std::pow(t8, -0.110e0_rt);
-       dd02  = -0.11e0_rt*4.59e-5_rt * std::pow(t8, -1.11e0_rt)*1.0e-8_rt;
+
+       if constexpr (do_derivatives) {
+           dd00  = (1.5e0_rt*7.75e5_rt*t812 + 3.85e0_rt*247.0e0_rt*std::pow(t8, 2.85e0_rt))*1.0e-8_rt;
+
+           dd01  = 1.4e0_rt*0.0240e0_rt * std::pow(t8, 0.4e0_rt)*1.0e-8_rt;
+           dd02  = -0.11e0_rt*4.59e-5_rt * std::pow(t8, -1.11e0_rt)*1.0e-8_rt;
+       }
 
        z     = std::pow(den, 0.656e0_rt);
        dum   = c00*rmi  + c01  + c02*z;

From 061f3d40e185e650db2fde0df09382311820b913 Mon Sep 17 00:00:00 2001
From: Michael Zingale <michael.zingale@stonybrook.edu>
Date: Sun, 10 Sep 2023 15:35:48 -0400
Subject: [PATCH 3/5] fix some compilation issues

---
 neutrinos/sneut5.H     | 6 +++---
 nse_solver/nse_check.H | 3 ++-
 2 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H
index 81449010fb..2a58ef4417 100644
--- a/neutrinos/sneut5.H
+++ b/neutrinos/sneut5.H
@@ -233,8 +233,8 @@ void sneut5(const Real temp, const Real den,
          a0,a1,a2,a3,b1,b2,c00,c01,c02,c03,c04,c05,c06,
          c10,c11,c12,c13,c14,c15,c16,c20,c21,c22,c23,c24,
          c25,c26,dd00,dd01,dd02,dd03,dd04,dd05,dd11,dd12,
-         dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,b,c,d,f0,
-         f1,deni,tempi,abari,zbari,f2,f3,z,ye,
+         dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,b,c,d{0.0},f0,
+         f1{0.0},deni,tempi,abari,zbari,f2,f3,z,ye,
          dum,dumdt,gum,gumdt,dumda,dumdz,gumda,gumdz;
 
     // pair production
@@ -599,9 +599,9 @@ void sneut5(const Real temp, const Real den,
        a3  = c/d;
        c00 = std::exp(-a3*a3);
 
-       f1  = -c00*2.0e0_rt*a3/d;
        fxy   = 1.05e0_rt + (a1 - b1)*c00;
        if constexpr (do_derivatives) {
+           f1  = -c00*2.0e0_rt*a3/d;
            c01 = f1*(dumdt + a3*0.25e0_rt*xnumdt);
            c03 = f1*(dumda + a3*0.25e0_rt*xnumda);
            c04 = f1*(dumdz + a3*0.25e0_rt*xnumdz);
diff --git a/nse_solver/nse_check.H b/nse_solver/nse_check.H
index 8ce878f196..1957dfe9e3 100644
--- a/nse_solver/nse_check.H
+++ b/nse_solver/nse_check.H
@@ -1070,7 +1070,8 @@ bool in_nse(burn_t& current_state, bool skip_molar_check=false) {
   composition(state);
 
   amrex::Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-  sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+  constexpr int do_derivatives{0};
+  sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 #else
   amrex::Real sneut = 0.0_rt;
 #endif

From 5d1a991db50ab6d83d17b493d6177dd9a17c2d21 Mon Sep 17 00:00:00 2001
From: Michael Zingale <michael.zingale@stonybrook.edu>
Date: Sun, 10 Sep 2023 20:56:12 -0400
Subject: [PATCH 4/5] more templating

---
 neutrinos/sneut5.H | 52 +++++++++++++++++++++++++++-------------------
 1 file changed, 31 insertions(+), 21 deletions(-)

diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H
index 2a58ef4417..56c7edd1f9 100644
--- a/neutrinos/sneut5.H
+++ b/neutrinos/sneut5.H
@@ -374,8 +374,8 @@ void sneut5(const Real temp, const Real den,
     zeta   = a1 * xlm1;
 
     if constexpr (do_derivatives) {
-        zetadt = -a1 * xlm2 * xldt;
         a2     = oneth * a1*rmi * xlm1;
+        zetadt = -a1 * xlm2 * xldt;
         zetada = a2 * rmda;
         zetadz = a2 * rmdz;
     }
@@ -395,18 +395,22 @@ void sneut5(const Real temp, const Real den,
     // equation 2.7
 
     a1     = 6.002e19_rt + 2.084e20_rt*zeta + 1.872e21_rt*zeta2;
-    a2     = 2.084e20_rt + 2.0e0_rt*1.872e21_rt*zeta;
 
     if (t9 < 10.0_rt) {
        b1     = exp(-5.5924e0_rt*zeta);
-       b2     = -b1*5.5924e0_rt;
+       if constexpr (do_derivatives) {
+           b2     = -b1*5.5924e0_rt;
+       }
     } else {
        b1     = exp(-4.9924e0_rt*zeta);
-       b2     = -b1*4.9924e0_rt;
+       if constexpr (do_derivatives) {
+           b2     = -b1*4.9924e0_rt;
+       }
     }
 
     xnum   = a1 * b1;
     if constexpr (do_derivatives) {
+        a2     = 2.084e20_rt + 2.0e0_rt*1.872e21_rt*zeta;
         c      = a2*b1 + a1*b2;
         xnumdt = c*zetadt;
         xnumda = c*zetada;
@@ -415,10 +419,14 @@ void sneut5(const Real temp, const Real den,
 
     if (t9 < 10.0_rt) {
        a1   = 9.383e-1_rt*xlm1 - 4.141e-1_rt*xlm2 + 5.829e-2_rt*xlm3;
-       a2   = -9.383e-1_rt*xlm2 + 2.0e0_rt*4.141e-1_rt*xlm3 - 3.0e0_rt*5.829e-2_rt*xlm4;
+       if constexpr (do_derivatives) {
+           a2   = -9.383e-1_rt*xlm2 + 2.0e0_rt*4.141e-1_rt*xlm3 - 3.0e0_rt*5.829e-2_rt*xlm4;
+       }
     } else {
        a1   = 1.2383e0_rt*xlm1 - 8.141e-1_rt*xlm2;
-       a2   = -1.2383e0_rt*xlm2 + 2.0e0_rt*8.141e-1_rt*xlm3;
+       if constexpr (do_derivatives) {
+           a2   = -1.2383e0_rt*xlm2 + 2.0e0_rt*8.141e-1_rt*xlm3;
+       }
     }
 
     xden   = zeta3 + a1;
@@ -439,18 +447,19 @@ void sneut5(const Real temp, const Real den,
 
     // equation 2.6
     a1     = 10.7480e0_rt*xl2 + 0.3967e0_rt*xlp5 + 1.005e0_rt;
-    a2     = xldt*(2.0e0_rt*10.7480e0_rt*xl + 0.5e0_rt*0.3967e0_rt*xlmp5);
     xnum   = 1.0e0_rt/a1;
-    xnumdt = -xnum*xnum*a2;
+    if constexpr (do_derivatives) {
+        a2     = xldt*(2.0e0_rt*10.7480e0_rt*xl + 0.5e0_rt*0.3967e0_rt*xlmp5);
+        xnumdt = -xnum*xnum*a2;
+    }
 
     a1     = 7.692e7_rt*xl3 + 9.715e6_rt*xlp5;
-    a2     = xldt*(3.0e0_rt*7.692e7_rt*xl2 + 0.5e0_rt*9.715e6_rt*xlmp5);
-
     c      = 1.0e0_rt/a1;
     b1     = 1.0e0_rt + rm*c;
 
     xden   = std::pow(b1, -0.3e0_rt);
     if constexpr (do_derivatives) {
+        a2     = xldt*(3.0e0_rt*7.692e7_rt*xl2 + 0.5e0_rt*9.715e6_rt*xlmp5);
         d      = -0.3e0_rt*xden/b1;
         xdendt = -d*rm*c*c*a2;
         xdenda = d*rmda*c;
@@ -484,11 +493,11 @@ void sneut5(const Real temp, const Real den,
     }
 
     a1      = tfac4*(1.0e0_rt + tfac3 * qpair);
-    a2      = tfac4*tfac3;
 
     a3      = spair;
     spair   = a1*a3;
     if constexpr (do_derivatives) {
+        a2      = tfac4*tfac3;
         spairdt = a1*spairdt + a2*qpairdt*a3;
         spairda = a1*spairda + a2*qpairda*a3;
         spairdz = a1*spairdz + a2*qpairdz*a3;
@@ -500,16 +509,16 @@ void sneut5(const Real temp, const Real den,
 
     a1   = 1.019e-6_rt*rm;
     a2   = std::pow(a1, twoth);
-    a3   = twoth*a2/a1;
 
     b1   = std::sqrt(1.0e0_rt + a2);
-    b2   = 1.0e0_rt/b1;
 
     c00  = 1.0e0_rt/(temp*temp*b1);
 
     gl2   = 1.1095e11_rt * rm * c00;
 
     if constexpr (do_derivatives) {
+        a3   = twoth*a2/a1;
+        b2   = 1.0e0_rt/b1;
         gl2dt = -2.0e0_rt*gl2*tempi;
         d     = rm*c00*b2*0.5e0_rt*b2*a3*1.019e-6_rt;
         gl2da = 1.1095e11_rt * (rmda*c00  - d*rmda);
@@ -534,15 +543,15 @@ void sneut5(const Real temp, const Real den,
 
     // equation 4.8
     a1   = 8.6e0_rt*gl2 + 1.35e0_rt*gl72;
-    a2   = 8.6e0_rt + 1.75e0_rt*1.35e0_rt*gl72*gum;
 
     b1   = 225.0e0_rt - 17.0e0_rt*gl + gl2;
-    b2   = -0.5e0_rt*17.0e0_rt*gl*gum + 1.0e0_rt;
 
     c    = 1.0e0_rt/b1;
     fl   = a1*c;
 
     if constexpr (do_derivatives) {
+        a2   = 8.6e0_rt + 1.75e0_rt*1.35e0_rt*gl72*gum;
+        b2   = -0.5e0_rt*17.0e0_rt*gl*gum + 1.0e0_rt;
         d    = (a2 - fl*b2)*c;
         fldt = d*gl2dt;
         flda = d*gl2da;
@@ -577,13 +586,14 @@ void sneut5(const Real temp, const Real den,
     } else {
 
        a1  = 0.39e0_rt - 1.25e0_rt*xnum - 0.35e0_rt*std::sin(4.5e0_rt*xnum);
-       a2  = -1.25e0_rt - 4.5e0_rt*0.35e0_rt*std::cos(4.5e0_rt*xnum);
 
        b1  = 0.3e0_rt * std::exp(-std::pow(4.5e0_rt*xnum + 0.9e0_rt, 2.0_rt));
-       b2  = -b1*2.0e0_rt*(4.5e0_rt*xnum + 0.9e0_rt)*4.5e0_rt;
 
        c   = amrex::min(0.0e0_rt, xden - 1.6e0_rt + 1.25e0_rt*xnum);
+
        if constexpr (do_derivatives) {
+           a2  = -1.25e0_rt - 4.5e0_rt*0.35e0_rt*std::cos(4.5e0_rt*xnum);
+           b2  = -b1*2.0e0_rt*(4.5e0_rt*xnum + 0.9e0_rt)*4.5e0_rt;
            if (c == 0.0_rt) {
                dumdt = 0.0e0_rt;
                dumda = 0.0e0_rt;
@@ -621,33 +631,33 @@ void sneut5(const Real temp, const Real den,
     }
 
     a2      = std::exp(-gl);
-    a3      = -0.5e0_rt*a2*gl*gum;
 
     a1      = splas;
     splas   = a2*a1;
     if constexpr (do_derivatives) {
+        a3      = -0.5e0_rt*a2*gl*gum;
         splasdt = a2*splasdt + a3*gl2dt*a1;
         splasda = a2*splasda + a3*gl2da*a1;
         splasdz = a2*splasdz + a3*gl2dz*a1;
     }
 
     a2      = gl6;
-    a3      = 3.0e0_rt*gl6*gum;
 
     a1      = splas;
     splas   = a2*a1;
     if constexpr (do_derivatives) {
+        a3      = 3.0e0_rt*gl6*gum;
         splasdt = a2*splasdt + a3*gl2dt*a1;
         splasda = a2*splasda + a3*gl2da*a1;
         splasdz = a2*splasdz + a3*gl2dz*a1;
     }
 
     a2      = 0.93153e0_rt * 3.0e21_rt * xl9;
-    a3      = 0.93153e0_rt * 3.0e21_rt * 9.0e0_rt*xl8*xldt;
 
     a1      = splas;
     splas   = a2*a1;
     if constexpr (do_derivatives) {
+        a3      = 0.93153e0_rt * 3.0e21_rt * 9.0e0_rt*xl8*xldt;
         splasdt = a2*splasdt + a3*a1;
         splasda = a2*splasda;
         splasdz = a2*splasdz;
@@ -921,11 +931,11 @@ void sneut5(const Real temp, const Real den,
     }
 
     a1      = tfac4*(1.0e0_rt - tfac3 * qphot);
-    a2      = -tfac4*tfac3;
 
     a3      = sphot;
     sphot   = a1*a3;
     if constexpr (do_derivatives) {
+        a2      = -tfac4*tfac3;
         sphotdt = a1*sphotdt + a2*qphotdt*a3;
         sphotda = a1*sphotda + a2*qphotda*a3;
         sphotdz = a1*sphotdz + a2*qphotdz*a3;

From f4f4b75ebf6fd47c7ed6c2e1a4af77085f637113 Mon Sep 17 00:00:00 2001
From: Michael Zingale <michael.zingale@stonybrook.edu>
Date: Thu, 14 Sep 2023 12:56:41 -0400
Subject: [PATCH 5/5] update test

---
 unit_test/test_neutrino_cooling/neutrino_util.cpp | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/unit_test/test_neutrino_cooling/neutrino_util.cpp b/unit_test/test_neutrino_cooling/neutrino_util.cpp
index 95a4799e96..efa6c70bb4 100644
--- a/unit_test/test_neutrino_cooling/neutrino_util.cpp
+++ b/unit_test/test_neutrino_cooling/neutrino_util.cpp
@@ -80,8 +80,10 @@ void neut_test_C(const Box& bx,
     Real dsnuda;
     Real dsnudz;
 
-    sneut5(temp_zone, dens_zone, abar, zbar,
-           snu, dsnudt, dsnudd, dsnuda, dsnudz);
+    constexpr int do_derivatives{1};
+
+    sneut5<do_derivatives>(temp_zone, dens_zone, abar, zbar,
+                           snu, dsnudt, dsnudd, dsnuda, dsnudz);
 
     sp(i, j, k, vars.isneut) = snu;
     sp(i, j, k, vars.isneutdt) = dsnudt;