diff --git a/src/particles/ReferenceParticle.H b/src/particles/ReferenceParticle.H
index 7f1b7ae19..226025b74 100644
--- a/src/particles/ReferenceParticle.H
+++ b/src/particles/ReferenceParticle.H
@@ -66,7 +66,7 @@ namespace impactx
             using namespace amrex::literals;
 
             amrex::ParticleReal const ref_gamma = -pt;
-            amrex::ParticleReal const ref_beta = sqrt(1.0_prt - 1.0_prt/pow(ref_gamma,2));
+            amrex::ParticleReal const ref_beta = std::sqrt(1.0_prt - 1.0_prt/pow(ref_gamma,2));
             return ref_beta;
         }
 
@@ -81,7 +81,7 @@ namespace impactx
             using namespace amrex::literals;
 
             amrex::ParticleReal const ref_gamma = -pt;
-            amrex::ParticleReal const ref_betagamma = sqrt(pow(ref_gamma, 2) - 1.0_prt);
+            amrex::ParticleReal const ref_betagamma = std::sqrt(std::pow(ref_gamma, 2) - 1.0_prt);
             return ref_betagamma;
         }
 
@@ -118,7 +118,7 @@ namespace impactx
             if (pt != 0.0_prt)
             {
                 pt = -kin_energy_MeV() / massE - 1.0_prt;
-                pz = sqrt(pow(pt, 2) - 1.0_prt);
+                pz = std::sqrt(std::pow(pt, 2) - 1.0_prt);
             }
 
             return *this;
@@ -155,7 +155,7 @@ namespace impactx
             px = 0.0;
             py = 0.0;
             pt = -kin_energy / mass_MeV() - 1.0_prt;
-            pz = sqrt(pow(pt, 2) - 1.0_prt);
+            pz = std::sqrt(std::pow(pt, 2) - 1.0_prt);
 
             return *this;
         }
@@ -171,7 +171,7 @@ namespace impactx
             using namespace amrex::literals;
 
             amrex::ParticleReal const ref_gamma = -pt;
-            amrex::ParticleReal const ref_betagamma = sqrt(pow(ref_gamma, 2) - 1.0_prt);
+            amrex::ParticleReal const ref_betagamma = std::sqrt(std::pow(ref_gamma, 2) - 1.0_prt);
             //amrex::ParticleReal const ref_rigidity = mass*ref_betagamma*(ablastr::constant::SI::c)/charge; //fails due to "charge"
             amrex::ParticleReal const ref_rigidity = mass*ref_betagamma*(ablastr::constant::SI::c)/(ablastr::constant::SI::q_e);
             return ref_rigidity;
diff --git a/src/particles/diagnostics/NonlinearLensInvariants.H b/src/particles/diagnostics/NonlinearLensInvariants.H
index 87a5f640e..9ce4eda3e 100644
--- a/src/particles/diagnostics/NonlinearLensInvariants.H
+++ b/src/particles/diagnostics/NonlinearLensInvariants.H
@@ -92,7 +92,7 @@ namespace impactx::diagnostics
             Complex const re1(1.0_prt, 0.0_prt);
             Complex const im1(0.0_prt, 1.0_prt);
 
-            // compute croot = sqrt(1-zeta**2)
+            // compute croot = std::sqrt(1-zeta**2)
             Complex croot = amrex::pow(zeta, 2);
             croot = re1 - croot;
             croot = amrex::sqrt(croot);
diff --git a/src/particles/distribution/Gaussian.H b/src/particles/distribution/Gaussian.H
index d6baf9640..24397c0b0 100644
--- a/src/particles/distribution/Gaussian.H
+++ b/src/particles/distribution/Gaussian.H
@@ -100,34 +100,34 @@ namespace impactx::distribution
 
             u1 = amrex::Random(engine);
             u2 = amrex::Random(engine);
-            ln1 = sqrt(-2_prt*log(u1));
-            x = ln1*cos(2_prt*pi*u2);
-            px = ln1*sin(2_prt*pi*u2);
+            ln1 = std::sqrt(-2_prt*std::log(u1));
+            x = ln1 * std::cos(2_prt*pi*u2);
+            px = ln1 * std::sin(2_prt*pi*u2);
 
             u1 = amrex::Random(engine);
             u2 = amrex::Random(engine);
-            ln1 = sqrt(-2_prt*log(u1));
-            y = ln1*cos(2_prt*pi*u2);
-            py = ln1*sin(2_prt*pi*u2);
+            ln1 = std::sqrt(-2_prt*std::log(u1));
+            y = ln1 * std::cos(2_prt*pi*u2);
+            py = ln1 * std::sin(2_prt*pi*u2);
 
             u1 = amrex::Random(engine);
             u2 = amrex::Random(engine);
-            ln1 = sqrt(-2_prt*log(u1));
-            t = ln1*cos(2_prt*pi*u2);
-            pt = ln1*sin(2_prt*pi*u2);
+            ln1 = std::sqrt(-2_prt*std::log(u1));
+            t = ln1 * std::cos(2_prt*pi*u2);
+            pt = ln1 * std::sin(2_prt*pi*u2);
 
             // Transform to produce the desired second moments/correlations:
-            root = sqrt(1.0_prt-m_muxpx*m_muxpx);
+            root = std::sqrt(1.0_prt-m_muxpx*m_muxpx);
             a1 = m_lambdaX * x / root;
             a2 = m_lambdaPx * (-m_muxpx * x / root + px);
             x = a1;
             px = a2;
-            root = sqrt(1.0_prt-m_muypy*m_muypy);
+            root = std::sqrt(1.0_prt-m_muypy*m_muypy);
             a1 = m_lambdaY * y / root;
             a2 = m_lambdaPy * (-m_muypy * y / root + py);
             y = a1;
             py = a2;
-            root = sqrt(1.0_prt-m_mutpt*m_mutpt);
+            root = std::sqrt(1.0_prt-m_mutpt*m_mutpt);
             a1 = m_lambdaT * t / root;
             a2 = m_lambdaPt * (-m_mutpt * t / root + pt);
             t = a1;
diff --git a/src/particles/distribution/KVdist.H b/src/particles/distribution/KVdist.H
index 87e343399..3554013eb 100644
--- a/src/particles/distribution/KVdist.H
+++ b/src/particles/distribution/KVdist.H
@@ -101,27 +101,27 @@ namespace impactx::distribution
             v = amrex::Random(engine);
             phi = amrex::Random(engine);
             phi = 2_prt*pi*phi;
-            r = sqrt(v);
-            x = r*cos(phi);
-            y = r*sin(phi);
+            r = std::sqrt(v);
+            x = r * std::cos(phi);
+            y = r * std::sin(phi);
 
             // Sample and transform to define (px,py):
             beta = amrex::Random(engine);
             beta = 2_prt*pi*beta;
-            p = sqrt(1_prt-pow(r,2));
-            px = p*cos(beta);
-            py = p*sin(beta);
+            p = std::sqrt(1_prt-pow(r,2));
+            px = p * std::cos(beta);
+            py = p * std::sin(beta);
 
             // Sample and transform to define (t,pt):
             t = amrex::Random(engine);
             t = 2.0_prt*(t-0.5_prt);
             u1 = amrex::Random(engine);
             u2 = amrex::Random(engine);
-            ln1 = sqrt(-2_prt*log(u1));
-            pt = ln1*cos(2_prt*pi*u2);
+            ln1 = std::sqrt(-2_prt*std::log(u1));
+            pt = ln1 * std::cos(2_prt*pi*u2);
 
             // Scale to produce the identity covariance matrix:
-            amrex::ParticleReal const c = sqrt(3.0_prt);
+            amrex::ParticleReal const c = std::sqrt(3.0_prt);
             x = 2_prt*x;
             y = 2_prt*y;
             t = c*t;
@@ -130,17 +130,17 @@ namespace impactx::distribution
             // pt = pt;
 
             // Transform to produce the desired second moments/correlations:
-            root = sqrt(1.0_prt-m_muxpx*m_muxpx);
+            root = std::sqrt(1.0_prt-m_muxpx*m_muxpx);
             a1 = m_lambdaX * x / root;
             a2 = m_lambdaPx * (-m_muxpx * x / root + px);
             x = a1;
             px = a2;
-            root = sqrt(1.0_prt-m_muypy*m_muypy);
+            root = std::sqrt(1.0_prt-m_muypy*m_muypy);
             a1 = m_lambdaY * y / root;
             a2 = m_lambdaPy * (-m_muypy * y / root + py);
             y = a1;
             py = a2;
-            root = sqrt(1.0_prt-m_mutpt*m_mutpt);
+            root = std::sqrt(1.0_prt-m_mutpt*m_mutpt);
             a1 = m_lambdaT * t / root;
             a2 = m_lambdaPt * (-m_mutpt * t / root + pt);
             t = a1;
diff --git a/src/particles/distribution/Kurth4D.H b/src/particles/distribution/Kurth4D.H
index cd254d2f9..d5b40e2d7 100644
--- a/src/particles/distribution/Kurth4D.H
+++ b/src/particles/distribution/Kurth4D.H
@@ -102,9 +102,9 @@ namespace impactx::distribution
             v = amrex::Random(engine);
             phi = amrex::Random(engine);
             phi = 2_prt*pi*phi;
-            r = sqrt(v);
-            x = r*cos(phi);
-            y = r*sin(phi);
+            r = std::sqrt(v);
+            x = r * std::cos(phi);
+            y = r * std::sin(phi);
 
             // Random samples used to define Lz:
             u = amrex::Random(engine);
@@ -113,25 +113,25 @@ namespace impactx::distribution
             // Random samples used to define pr:
             alpha = amrex::Random(engine);
             alpha = pi*alpha;
-            pmax = 1.0_prt - pow((Lz/r),2) - pow(r,2) + pow(Lz,2);
-            pmax = sqrt(pmax);
-            pr = pmax*cos(alpha);
+            pmax = 1.0_prt - std::pow((Lz/r),2) - std::pow(r,2) + std::pow(Lz,2);
+            pmax = std::sqrt(pmax);
+            pr = pmax * std::cos(alpha);
             pphi = Lz/r;
 
             // Transformations used to obtain (px,py):
-            px = pr*cos(phi)-pphi*sin(phi);
-            py = pr*sin(phi)+pphi*cos(phi);
+            px = pr * std::cos(phi)-pphi * std::sin(phi);
+            py = pr * std::sin(phi)+pphi * std::cos(phi);
 
             // Sample and transform to define (t,pt):
             t = amrex::Random(engine);
             t = 2.0_prt*(t-0.5_prt);
             u1 = amrex::Random(engine);
             u2 = amrex::Random(engine);
-            ln1 = sqrt(-2_prt*log(u1));
-            pt = ln1*cos(2_prt*pi*u2);
+            ln1 = std::sqrt(-2_prt*std::log(u1));
+            pt = ln1 * std::cos(2_prt*pi*u2);
 
             // Scale to produce the identity covariance matrix:
-            amrex::ParticleReal const c = sqrt(3.0_prt);
+            amrex::ParticleReal const c = std::sqrt(3.0_prt);
             x = 2_prt*x;
             y = 2_prt*y;
             t = c*t;
@@ -140,17 +140,17 @@ namespace impactx::distribution
             // pt = pt;
 
             // Transform to produce the desired second moments/correlations:
-            root = sqrt(1.0_prt-m_muxpx*m_muxpx);
+            root = std::sqrt(1.0_prt-m_muxpx*m_muxpx);
             a1 = m_lambdaX * x / root;
             a2 = m_lambdaPx * (-m_muxpx * x / root + px);
             x = a1;
             px = a2;
-            root = sqrt(1.0_prt-m_muypy*m_muypy);
+            root = std::sqrt(1.0_prt-m_muypy*m_muypy);
             a1 = m_lambdaY * y / root;
             a2 = m_lambdaPy * (-m_muypy * y / root + py);
             y = a1;
             py = a2;
-            root = sqrt(1.0_prt-m_mutpt*m_mutpt);
+            root = std::sqrt(1.0_prt-m_mutpt*m_mutpt);
             a1 = m_lambdaT * t / root;
             a2 = m_lambdaPt * (-m_mutpt * t / root + pt);
             t = a1;
diff --git a/src/particles/distribution/Kurth6D.H b/src/particles/distribution/Kurth6D.H
index e3d1b8c42..2e29bd606 100644
--- a/src/particles/distribution/Kurth6D.H
+++ b/src/particles/distribution/Kurth6D.H
@@ -104,14 +104,14 @@ namespace impactx::distribution
             v = amrex::Random(engine);
             costheta = amrex::Random(engine);
             costheta = 2_prt*(costheta-0.5_prt);
-            sintheta = sqrt(1_prt-pow(costheta,2));
+            sintheta = std::sqrt(1_prt-pow(costheta,2));
             phi = amrex::Random(engine);
             phi = 2_prt*pi*phi;
 
             // Transformations for (x,y,t):
-            r = pow(v,1_prt/3_prt);
-            x = r*sintheta*cos(phi);
-            y = r*sintheta*sin(phi);
+            r = std::pow(v,1_prt/3_prt);
+            x = r*sintheta * std::cos(phi);
+            y = r*sintheta * std::sin(phi);
             t = r*costheta;
 
             // Random samples used to define L:
@@ -121,23 +121,23 @@ namespace impactx::distribution
             // Random samples used to define pr:
             alpha = amrex::Random(engine);
             alpha = pi*alpha;
-            pmax = 1_prt - pow(L/r,2) - pow(r,2) + pow(L,2);
-            pmax = sqrt(pmax);
-            pr = pmax*cos(alpha);
+            pmax = 1_prt - std::pow(L/r,2) - std::pow(r,2) + std::pow(L,2);
+            pmax = std::sqrt(pmax);
+            pr = pmax * std::cos(alpha);
 
             // Random samples used to define ptangent:
             beta = amrex::Random(engine);
             beta = 2_prt*pi*beta;
-            p1 = L/r*cos(beta);  // This is phi component
-            p2 = L/r*sin(beta);  // This is theta component
+            p1 = L/r * std::cos(beta);  // This is phi component
+            p2 = L/r * std::sin(beta);  // This is theta component
 
             // Transformation from spherical to Cartesian coord.:
-            px = pr*sintheta*cos(phi) + p2*costheta*cos(phi) - p1*sin(phi);
-            py = pr*sintheta*sin(phi) + p2*costheta*sin(phi) + p1*cos(phi);
+            px = pr*sintheta * std::cos(phi) + p2*costheta * std::cos(phi) - p1 * std::sin(phi);
+            py = pr*sintheta * std::sin(phi) + p2*costheta * std::sin(phi) + p1 * std::cos(phi);
             pt = pr*costheta - p2*sintheta;
 
             // Scale to produce the identity covariance matrix:
-            amrex::ParticleReal const c = sqrt(5.0_prt);
+            amrex::ParticleReal const c = std::sqrt(5.0_prt);
             x = c*x;
             y = c*y;
             t = c*t;
@@ -146,17 +146,17 @@ namespace impactx::distribution
             pt = c*pt;
 
             // Transform to produce the desired second moments/correlations:
-            root = sqrt(1.0_prt-m_muxpx*m_muxpx);
+            root = std::sqrt(1.0_prt-m_muxpx*m_muxpx);
             a1 = m_lambdaX * x / root;
             a2 = m_lambdaPx * (-m_muxpx * x / root + px);
             x = a1;
             px = a2;
-            root = sqrt(1.0_prt-m_muypy*m_muypy);
+            root = std::sqrt(1.0_prt-m_muypy*m_muypy);
             a1 = m_lambdaY * y / root;
             a2 = m_lambdaPy * (-m_muypy * y / root + py);
             y = a1;
             py = a2;
-            root = sqrt(1.0_prt-m_mutpt*m_mutpt);
+            root = std::sqrt(1.0_prt-m_mutpt*m_mutpt);
             a1 = m_lambdaT * t / root;
             a2 = m_lambdaPt * (-m_mutpt * t / root + pt);
             t = a1;
diff --git a/src/particles/distribution/Semigaussian.H b/src/particles/distribution/Semigaussian.H
index 63d985657..f6c3deaa1 100644
--- a/src/particles/distribution/Semigaussian.H
+++ b/src/particles/distribution/Semigaussian.H
@@ -103,14 +103,14 @@ namespace impactx::distribution
             phi = amrex::Random(engine);
             phi = 2_prt*pi*phi;
             v = amrex::Random(engine);
-            r = sqrt(v);
-            x = r*cos(phi);
-            y = r*sin(phi);
+            r = std::sqrt(v);
+            x = r * std::cos(phi);
+            y = r * std::sin(phi);
             t = amrex::Random(engine);
             t = 2_prt*(t-0.5_prt);
 
             // Scale to produce the identity covariance matrix:
-            amrex::ParticleReal const c = sqrt(3.0_prt);
+            amrex::ParticleReal const c = std::sqrt(3.0_prt);
             x = 2_prt*x;
             y = 2_prt*y;
             t = c*t;
@@ -119,27 +119,27 @@ namespace impactx::distribution
 
             u1 = amrex::Random(engine);
             u2 = amrex::Random(engine);
-            ln1 = sqrt(-2_prt*log(u1));
-            px = ln1*cos(2_prt*pi*u2);
-            py = ln1*sin(2_prt*pi*u2);
+            ln1 = std::sqrt(-2_prt*std::log(u1));
+            px = ln1 * std::cos(2_prt*pi*u2);
+            py = ln1 * std::sin(2_prt*pi*u2);
 
             u1 = amrex::Random(engine);
             u2 = amrex::Random(engine);
-            ln1 = sqrt(-2_prt*log(u1));
-            pt = ln1*cos(2_prt*pi*u2);
+            ln1 = std::sqrt(-2_prt*std::log(u1));
+            pt = ln1 * std::cos(2_prt*pi*u2);
 
             // Transform to produce the desired second moments/correlations:
-            root = sqrt(1.0_prt-m_muxpx*m_muxpx);
+            root = std::sqrt(1.0_prt-m_muxpx*m_muxpx);
             a1 = m_lambdaX * x / root;
             a2 = m_lambdaPx * (-m_muxpx * x / root + px);
             x = a1;
             px = a2;
-            root = sqrt(1.0_prt-m_muypy*m_muypy);
+            root = std::sqrt(1.0_prt-m_muypy*m_muypy);
             a1 = m_lambdaY * y / root;
             a2 = m_lambdaPy * (-m_muypy * y / root + py);
             y = a1;
             py = a2;
-            root = sqrt(1.0_prt-m_mutpt*m_mutpt);
+            root = std::sqrt(1.0_prt-m_mutpt*m_mutpt);
             a1 = m_lambdaT * t / root;
             a2 = m_lambdaPt * (-m_mutpt * t / root + pt);
             t = a1;
diff --git a/src/particles/distribution/Thermal.H b/src/particles/distribution/Thermal.H
index 63169332a..3017cad31 100644
--- a/src/particles/distribution/Thermal.H
+++ b/src/particles/distribution/Thermal.H
@@ -97,7 +97,7 @@ namespace distribution
             amrex::ParticleReal q_e = refpart.charge_qe();
 
             // Set space charge intensity
-            m_Cintensity = q_e*bunch_charge/(pow(bg,2)*Erest*ablastr::constant::SI::ep0);
+            m_Cintensity = q_e*bunch_charge/(std::pow(bg,2)*Erest*ablastr::constant::SI::ep0);
 
             // Set minimum and maximum radius
             amrex::ParticleReal r_scale = matched_scale_radius();
@@ -106,8 +106,8 @@ namespace distribution
             // amrex::PrintToFile("equilibrium_params.out") << r_scale << " " << data.Cintensity << "\n";
 
             // Scale the parameters p1 and p2
-            amrex::ParticleReal rt2pi = sqrt(2.0_prt*pi);
-            amrex::ParticleReal p_scale = pow(r_scale*rt2pi,-3);
+            amrex::ParticleReal rt2pi = std::sqrt(2.0_prt*pi);
+            amrex::ParticleReal p_scale = std::pow(r_scale*rt2pi,-3);
             m_p1 = m_p1*p_scale;
             m_p2 = m_p2*p_scale;
 
@@ -162,7 +162,7 @@ namespace distribution
             amrex::ParticleReal k = m_k;
             amrex::ParticleReal kT = (1.0_prt - m_w) * m_T1 + m_w * m_T2;
             amrex::ParticleReal a = m_Cintensity/(4.0_prt*pi*5.0_prt*sqrt(5.0_prt));
-            amrex::ParticleReal rscale = sqrt(kT + pow(a*k,2.0/3.0))/k;
+            amrex::ParticleReal rscale = std::sqrt(kT + std::pow(a*k,2.0/3.0))/k;
 
             return rscale;
         }
@@ -252,17 +252,17 @@ namespace distribution
 
             // Define intermediate quantities
             amrex::ParticleReal potential = 0.0_prt;
-            potential = pow(k*r,2.0)/2.0_prt + c1*phi1 + c2*phi2;
-            amrex::ParticleReal Pdensity1 = m_p1*exp(-potential/T1);
-            amrex::ParticleReal Pdensity2 = m_p2*exp(-potential/T2);
+            potential = std::pow(k*r,2.0)/2.0_prt + c1*phi1 + c2*phi2;
+            amrex::ParticleReal Pdensity1 = m_p1 * std::exp(-potential/T1);
+            amrex::ParticleReal Pdensity2 = m_p2 * std::exp(-potential/T2);
             // amrex::ParticleReal Pdensity_tot = (1.0_prt-w)*Pdensity1 + w*Pdensity2;
             // amrex::PrintToFile("Pdensity.out") << reval << " " << Pdensity_tot << "\n";
 
             // Apply map to update f1 and f2:
             m_phi1 = phi1;
             m_phi2 = phi2;
-            m_f1 = f1 + tau*4.0_prt*pi*pow(r,2.0)*Pdensity1;
-            m_f2 = f2 + tau*4.0_prt*pi*pow(r,2.0)*Pdensity2;
+            m_f1 = f1 + tau*4.0_prt*pi * std::pow(r,2.0)*Pdensity1;
+            m_f2 = f2 + tau*4.0_prt*pi * std::pow(r,2.0)*Pdensity2;
             reval = r;
         }
     };
@@ -371,28 +371,28 @@ namespace distribution
             // Generate six standard normal random variables using Box-Muller:
             u1 = amrex::Random(engine);
             u2 = amrex::Random(engine);
-            ln1 = sqrt(-2_prt*log(u1));
-            g1 = ln1*cos(2_prt*pi*u2);
-            g2 = ln1*sin(2_prt*pi*u2);
+            ln1 = std::sqrt(-2_prt*std::log(u1));
+            g1 = ln1 * std::cos(2_prt*pi*u2);
+            g2 = ln1 * std::sin(2_prt*pi*u2);
             u1 = amrex::Random(engine);
             u2 = amrex::Random(engine);
-            ln1 = sqrt(-2_prt*log(u1));
-            g3 = ln1*cos(2_prt*pi*u2);
-            g4 = ln1*sin(2_prt*pi*u2);
+            ln1 = std::sqrt(-2_prt*std::log(u1));
+            g3 = ln1 * std::cos(2_prt*pi*u2);
+            g4 = ln1 * std::sin(2_prt*pi*u2);
             u1 = amrex::Random(engine);
             u2 = amrex::Random(engine);
-            ln1 = sqrt(-2_prt*log(u1));
-            g5 = ln1*cos(2_prt*pi*u2);
-            g6 = ln1*sin(2_prt*pi*u2);
+            ln1 = std::sqrt(-2_prt*std::log(u1));
+            g5 = ln1 * std::cos(2_prt*pi*u2);
+            g6 = ln1 * std::sin(2_prt*pi*u2);
 
             // Scale the last three variables to produce the momenta:
             amrex::ParticleReal kT = (uhalo > m_w) ? m_T1 : m_T2;  // select core or halo value
-            px = sqrt(kT)*g4;
-            py = sqrt(kT)*g5;
-            pz = sqrt(kT)*g6;
+            px = std::sqrt(kT)*g4;
+            py = std::sqrt(kT)*g5;
+            pz = std::sqrt(kT)*g6;
 
             // Normalize the first three variables to produce uniform samples on the unit 3-sphere:
-            norm = sqrt(g1*g1+g2*g2+g3*g3);
+            norm = std::sqrt(g1*g1+g2*g2+g3*g3);
             g1 /= norm;
             g2 /= norm;
             g3 /= norm;
diff --git a/src/particles/distribution/Triangle.H b/src/particles/distribution/Triangle.H
index a1f40dc14..3625fa4fa 100644
--- a/src/particles/distribution/Triangle.H
+++ b/src/particles/distribution/Triangle.H
@@ -99,29 +99,29 @@ namespace impactx::distribution
             // Sample the t coordinate for a ramped triangular profile (unit
             // variance):
             u0 = amrex::Random(engine);
-            t = sqrt(2_prt)*(2_prt-3_prt*sqrt(u0));
+            t = std::sqrt(2_prt)*(2_prt-3_prt*sqrt(u0));
 
             // Generate five standard normal random variables using Box-Muller:
             u1 = amrex::Random(engine);
             u2 = amrex::Random(engine);
-            ln1 = sqrt(-2_prt*log(u1));
-            g1 = ln1*cos(2_prt*pi*u2);
-            g2 = ln1*sin(2_prt*pi*u2);
+            ln1 = std::sqrt(-2_prt*std::log(u1));
+            g1 = ln1 * std::cos(2_prt*pi*u2);
+            g2 = ln1 * std::sin(2_prt*pi*u2);
             u1 = amrex::Random(engine);
             u2 = amrex::Random(engine);
-            ln1 = sqrt(-2_prt*log(u1));
-            g3 = ln1*cos(2_prt*pi*u2);
-            g4 = ln1*sin(2_prt*pi*u2);
+            ln1 = std::sqrt(-2_prt*std::log(u1));
+            g3 = ln1 * std::cos(2_prt*pi*u2);
+            g4 = ln1 * std::sin(2_prt*pi*u2);
             u1 = amrex::Random(engine);
             u2 = amrex::Random(engine);
-            ln1 = sqrt(-2_prt*log(u1));
-            g5 = ln1*cos(2_prt*pi*u2);
+            ln1 = std::sqrt(-2_prt*std::log(u1));
+            g5 = ln1 * std::cos(2_prt*pi*u2);
 
             // Use one of these normal random variables for pt:
             pt = g5;
 
             // Normalize the rest to produce uniform samples on the unit sphere:
-            norm = sqrt(g1*g1+g2*g2+g3*g3+g4*g4);
+            norm = std::sqrt(g1*g1+g2*g2+g3*g3+g4*g4);
             g1 /= norm;
             g2 /= norm;
             g3 /= norm;
@@ -130,24 +130,24 @@ namespace impactx::distribution
             // Scale to produce uniform samples in a 4D ball (unit variance):
             d = 4_prt;  // unit ball dimension
             u1 = amrex::Random(engine);   // uniform sample
-            u2 = sqrt(d+2_prt)*pow(u1,1_prt/d);
+            u2 = std::sqrt(d+2_prt) * std::pow(u1,1_prt/d);
             x = g1*u2;
             y = g2*u2;
             px = g3*u2;
             py = g4*u2;
 
             // Transform to produce the desired second moments/correlations:
-            root = sqrt(1.0_prt-m_muxpx*m_muxpx);
+            root = std::sqrt(1.0_prt-m_muxpx*m_muxpx);
             a1 = m_lambdaX * x / root;
             a2 = m_lambdaPx * (-m_muxpx * x / root + px);
             x = a1;
             px = a2;
-            root = sqrt(1.0_prt-m_muypy*m_muypy);
+            root = std::sqrt(1.0_prt-m_muypy*m_muypy);
             a1 = m_lambdaY * y / root;
             a2 = m_lambdaPy * (-m_muypy * y / root + py);
             y = a1;
             py = a2;
-            root = sqrt(1.0_prt-m_mutpt*m_mutpt);
+            root = std::sqrt(1.0_prt-m_mutpt*m_mutpt);
             a1 = m_lambdaT * t / root;
             a2 = m_lambdaPt * (-m_mutpt * t / root + pt);
             t = a1;
diff --git a/src/particles/elements/Aperture.H b/src/particles/elements/Aperture.H
index 2f1e6e9eb..199712fb9 100644
--- a/src/particles/elements/Aperture.H
+++ b/src/particles/elements/Aperture.H
@@ -102,13 +102,13 @@ namespace impactx
             switch (m_shape)
             {
                 case Shape::rectangular :  // default
-                  if (pow(u,2)>1 || pow(v,2) > 1_prt) {
+                  if (std::pow(u,2)>1 || std::pow(v,2) > 1_prt) {
                       amrex::ParticleIDWrapper{idcpu}.make_invalid();
                   }
                   break;
 
                case Shape::elliptical :
-                  if (pow(u,2)+pow(v,2) > 1_prt) {
+                  if (std::pow(u,2) + std::pow(v,2) > 1_prt) {
                       amrex::ParticleIDWrapper{idcpu}.make_invalid();
                   }
                   break;
diff --git a/src/particles/elements/Buncher.H b/src/particles/elements/Buncher.H
index 1dc83d4ed..731682c0f 100644
--- a/src/particles/elements/Buncher.H
+++ b/src/particles/elements/Buncher.H
@@ -86,7 +86,7 @@ namespace impactx
 
             // access reference particle values to find (beta*gamma)^2
             amrex::ParticleReal const pt_ref = refpart.pt;
-            amrex::ParticleReal const betgam2 = pow(pt_ref, 2) - 1.0_prt;
+            amrex::ParticleReal const betgam2 = std::pow(pt_ref, 2) - 1.0_prt;
 
             // intialize output values of momenta
             amrex::ParticleReal pxout = px;
diff --git a/src/particles/elements/CFbend.H b/src/particles/elements/CFbend.H
index 2db863ce0..b9007986c 100644
--- a/src/particles/elements/CFbend.H
+++ b/src/particles/elements/CFbend.H
@@ -109,18 +109,18 @@ namespace impactx
 
             // access reference particle values to find beta*gamma^2
             amrex::ParticleReal const pt_ref = refpart.pt;
-            amrex::ParticleReal const betgam2 = pow(pt_ref, 2) - 1.0_prt;
-            amrex::ParticleReal const bet = sqrt(betgam2/(1.0_prt + betgam2));
+            amrex::ParticleReal const betgam2 = std::pow(pt_ref, 2) - 1.0_prt;
+            amrex::ParticleReal const bet = std::sqrt(betgam2/(1.0_prt + betgam2));
 
             // update horizontal and longitudinal phase space variables
-            amrex::ParticleReal const gx = m_k + pow(m_rc,-2);
-            amrex::ParticleReal const omegax = sqrt(std::abs(gx));
+            amrex::ParticleReal const gx = m_k + std::pow(m_rc,-2);
+            amrex::ParticleReal const omegax = std::sqrt(std::abs(gx));
 
             if(gx > 0.0) {
                 // calculate expensive terms once
                 auto const [sinx, cosx] = amrex::Math::sincos(omegax * slice_ds);
                 amrex::ParticleReal const r56 = slice_ds/betgam2
-                    + (sinx - omegax*slice_ds)/(gx*omegax*pow(bet,2)*pow(m_rc,2));
+                    + (sinx - omegax*slice_ds)/(gx*omegax * std::pow(bet,2) * std::pow(m_rc,2));
 
                 // advance position and momentum (focusing)
                 x = cosx*xout + sinx/omegax*px - (1.0_prt - cosx)/(gx*bet*m_rc)*pt;
@@ -131,10 +131,10 @@ namespace impactx
                 ptout = pt;
             } else {
                 // calculate expensive terms once
-                amrex::ParticleReal const sinhx = sinh(omegax * slice_ds);
-                amrex::ParticleReal const coshx = cosh(omegax * slice_ds);
+                amrex::ParticleReal const sinhx = std::sinh(omegax * slice_ds);
+                amrex::ParticleReal const coshx = std::cosh(omegax * slice_ds);
                 amrex::ParticleReal const r56 = slice_ds/betgam2
-                    + (sinhx - omegax*slice_ds)/(gx*omegax*pow(bet,2)*pow(m_rc,2));
+                    + (sinhx - omegax*slice_ds)/(gx*omegax * std::pow(bet,2) * std::pow(m_rc,2));
 
                 // advance position and momentum (defocusing)
                 x = coshx*xout + sinhx/omegax*px - (1.0_prt - coshx)/(gx*bet*m_rc)*pt;
@@ -147,7 +147,7 @@ namespace impactx
 
             // update vertical phase space variables
             amrex::ParticleReal const gy = -m_k;
-            amrex::ParticleReal const omegay = sqrt(std::abs(gy));
+            amrex::ParticleReal const omegay = std::sqrt(std::abs(gy));
 
             if(gy > 0.0) {
                 // calculate expensive terms once
@@ -159,8 +159,8 @@ namespace impactx
 
             } else {
                 // calculate expensive terms once
-                amrex::ParticleReal const sinhy = sinh(omegay * slice_ds);
-                amrex::ParticleReal const coshy = cosh(omegay * slice_ds);
+                amrex::ParticleReal const sinhy = std::sinh(omegay * slice_ds);
+                amrex::ParticleReal const coshy = std::cosh(omegay * slice_ds);
 
                 // advance position and momentum (defocusing)
                 y = coshy*yout + sinhy/omegay*py;
@@ -201,7 +201,7 @@ namespace impactx
 
             // assign intermediate parameter
             amrex::ParticleReal const theta = slice_ds/m_rc;
-            amrex::ParticleReal const B = sqrt(pow(pt,2)-1.0_prt)/m_rc;
+            amrex::ParticleReal const B = std::sqrt(std::pow(pt,2)-1.0_prt)/m_rc;
 
             // calculate expensive terms once
             auto const [sin_theta, cos_theta] = amrex::Math::sincos(theta);
diff --git a/src/particles/elements/ChrDrift.H b/src/particles/elements/ChrDrift.H
index 35edd625a..3080b41cf 100644
--- a/src/particles/elements/ChrDrift.H
+++ b/src/particles/elements/ChrDrift.H
@@ -106,7 +106,7 @@ namespace impactx
 
             // compute particle momentum deviation delta + 1
             amrex::ParticleReal delta1;
-            delta1 = sqrt(1_prt - 2_prt*pt/bet + pow(pt,2));
+            delta1 = std::sqrt(1_prt - 2_prt*pt/bet + std::pow(pt,2));
 
             // advance transverse position and momentum (drift)
             x = xout + slice_ds * px / delta1;
@@ -115,12 +115,12 @@ namespace impactx
             // pyout = py;
 
             // the corresponding symplectic update to t
-            amrex::ParticleReal term = 2_prt*pow(pt,2)+pow(px,2)+pow(py,2);
-            term = 2_prt - 4_prt*bet*pt + pow(bet,2)*term;
-            term = -2_prt + pow(gam,2)*term;
+            amrex::ParticleReal term = 2_prt * std::pow(pt,2) + std::pow(px,2) + std::pow(py,2);
+            term = 2_prt - 4_prt*bet*pt + std::pow(bet,2)*term;
+            term = -2_prt + std::pow(gam,2)*term;
             term = (-1_prt+bet*pt)*term;
-            term = term/(2_prt*pow(bet,3)*pow(gam,2));
-            t = tout - slice_ds * (1_prt / bet + term / pow(delta1, 3));
+            term = term/(2_prt * std::pow(bet,3) * std::pow(gam,2));
+            t = tout - slice_ds * (1_prt / bet + term /std::pow(delta1, 3));
             // ptout = pt;
 
             // assign updated momenta
@@ -156,7 +156,7 @@ namespace impactx
             amrex::ParticleReal const slice_ds = m_ds / nslice();
 
             // assign intermediate parameter
-            amrex::ParticleReal const step = slice_ds / sqrt(pow(pt,2)-1.0_prt);
+            amrex::ParticleReal const step = slice_ds /std::sqrt(std::pow(pt,2)-1.0_prt);
 
             // advance position and momentum (drift)
             refpart.x = x + step*px;
diff --git a/src/particles/elements/ChrPlasmaLens.H b/src/particles/elements/ChrPlasmaLens.H
index 0d95d7583..ca389b2e4 100644
--- a/src/particles/elements/ChrPlasmaLens.H
+++ b/src/particles/elements/ChrPlasmaLens.H
@@ -110,12 +110,12 @@ namespace impactx
 
             // compute particle momentum deviation delta + 1
             amrex::ParticleReal delta1;
-            delta1 = sqrt(1_prt - 2_prt*pt/bet + pow(pt,2));
+            delta1 = std::sqrt(1_prt - 2_prt*pt/bet + std::pow(pt,2));
             amrex::ParticleReal const delta = delta1 - 1_prt;
 
             // compute phase advance per unit length in s (in rad/m)
             // chromatic dependence on delta is included
-            amrex::ParticleReal const omega = sqrt(std::abs(g)/delta1);
+            amrex::ParticleReal const omega = std::sqrt(std::abs(g)/delta1);
 
             // initialize output values
             amrex::ParticleReal xout = x;
@@ -148,13 +148,13 @@ namespace impactx
             amrex::ParticleReal const t0 = t - term*slice_ds/delta1;
 
             amrex::ParticleReal const w = omega*delta1;
-            amrex::ParticleReal const term1 = -(pow(p2,2)-pow(q2,2)*pow(w,2))*sin(2_prt*slice_ds*omega);
-            amrex::ParticleReal const term2 = -(pow(p1,2)-pow(q1,2)*pow(w,2))*sin(2_prt*slice_ds*omega);
-            amrex::ParticleReal const term3 = -2_prt*q2*p2*w*cos(2_prt*slice_ds*omega);
-            amrex::ParticleReal const term4 = -2_prt*q1*p1*w*cos(2_prt*slice_ds*omega);
+            amrex::ParticleReal const term1 = -(std::pow(p2,2)-pow(q2,2) * std::pow(w,2)) * std::sin(2_prt*slice_ds*omega);
+            amrex::ParticleReal const term2 = -(std::pow(p1,2)-pow(q1,2) * std::pow(w,2)) * std::sin(2_prt*slice_ds*omega);
+            amrex::ParticleReal const term3 = -2_prt*q2*p2*w * std::cos(2_prt*slice_ds*omega);
+            amrex::ParticleReal const term4 = -2_prt*q1*p1*w * std::cos(2_prt*slice_ds*omega);
             amrex::ParticleReal const term5 = 2_prt*omega*(q1*p1*delta1 + q2*p2*delta1
-                                        -(pow(p1,2)+pow(p2,2))*slice_ds - (pow(q1,2)+pow(q2,2))*pow(w,2)*slice_ds);
-            t = t0 + (-1_prt+bet*pt)/(8_prt*bet*pow(delta1,3)*omega)
+                                        -(std::pow(p1,2) + std::pow(p2,2))*slice_ds - (std::pow(q1,2) + std::pow(q2,2)) * std::pow(w,2)*slice_ds);
+            t = t0 + (-1_prt+bet*pt)/(8_prt*bet * std::pow(delta1,3)*omega)
                      *(term1+term2+term3+term4+term5);
 
             // ptout = pt;
@@ -194,7 +194,7 @@ namespace impactx
             amrex::ParticleReal const slice_ds = m_ds / nslice();
 
             // assign intermediate parameter
-            amrex::ParticleReal const step = slice_ds / sqrt(pow(pt,2)-1.0_prt);
+            amrex::ParticleReal const step = slice_ds /std::sqrt(std::pow(pt,2)-1.0_prt);
 
             // advance position and momentum (straight element)
             refpart.x = x + step*px;
diff --git a/src/particles/elements/ChrQuad.H b/src/particles/elements/ChrQuad.H
index 2c86cc809..04a0787f4 100644
--- a/src/particles/elements/ChrQuad.H
+++ b/src/particles/elements/ChrQuad.H
@@ -118,12 +118,12 @@ namespace impactx
 
             // compute particle momentum deviation delta + 1
             amrex::ParticleReal delta1;
-            delta1 = sqrt(1_prt - 2_prt*pt/bet + pow(pt,2));
+            delta1 = std::sqrt(1_prt - 2_prt*pt/bet + std::pow(pt,2));
             amrex::ParticleReal const delta = delta1 - 1_prt;
 
             // compute phase advance per unit length in s (in rad/m)
             // chromatic dependence on delta is included
-            amrex::ParticleReal const omega = sqrt(std::abs(g)/delta1);
+            amrex::ParticleReal const omega = std::sqrt(std::abs(g)/delta1);
 
             // intialize output values of momenta
             amrex::ParticleReal pxout = px;
@@ -138,23 +138,23 @@ namespace impactx
 
             if(g > 0.0) {
                // advance transverse position and momentum (focusing quad)
-               x = cos(omega*slice_ds) * xout +
-                                   sin(omega*slice_ds)/(omega*delta1)*px;
-               pxout = -omega * delta1 * sin(omega*slice_ds) * xout + cos(omega * slice_ds) * px;
+               x = std::cos(omega*slice_ds) * xout +
+                                  std::sin(omega*slice_ds)/(omega*delta1)*px;
+               pxout = -omega * delta1 *std::sin(omega*slice_ds) * xout + std::cos(omega * slice_ds) * px;
 
-               y = cosh(omega*slice_ds) * yout +
-                                   sinh(omega*slice_ds)/(omega*delta1)*py;
-               pyout = omega * delta1 * sinh(omega*slice_ds) * yout + cosh(omega * slice_ds) * py;
+               y = std::cosh(omega*slice_ds) * yout +
+                                   std::sinh(omega*slice_ds)/(omega*delta1)*py;
+               pyout = omega * delta1 * std::sinh(omega*slice_ds) * yout + std::cosh(omega * slice_ds) * py;
 
             } else {
                // advance transverse position and momentum (defocusing quad)
-               x = cosh(omega*slice_ds) * xout +
-                                   sinh(omega*slice_ds)/(omega*delta1)*px;
-               pxout = omega * delta1 * sinh(omega*slice_ds) * xout + cosh(omega * slice_ds) * px;
+               x = std::cosh(omega*slice_ds) * xout +
+                                   std::sinh(omega*slice_ds)/(omega*delta1)*px;
+               pxout = omega * delta1 * std::sinh(omega*slice_ds) * xout + std::cosh(omega * slice_ds) * px;
 
-               y = cos(omega*slice_ds) * yout +
-                                   sin(omega*slice_ds)/(omega*delta1)*py;
-               pyout = -omega * delta1 * sin(omega*slice_ds) * yout + cos(omega * slice_ds) * py;
+               y = std::cos(omega*slice_ds) * yout +
+                                  std::sin(omega*slice_ds)/(omega*delta1)*py;
+               pyout = -omega * delta1 *std::sin(omega*slice_ds) * yout + std::cos(omega * slice_ds) * py;
 
                q1 = yout;
                q2 = xout;
@@ -169,13 +169,13 @@ namespace impactx
             amrex::ParticleReal const t0 = tout - term * slice_ds / delta1;
 
             amrex::ParticleReal const w = omega*delta1;
-            amrex::ParticleReal const term1 = -(pow(p2,2)+pow(q2,2)*pow(w,2))*sinh(2_prt*slice_ds*omega);
-            amrex::ParticleReal const term2 = -(pow(p1,2)-pow(q1,2)*pow(w,2))*sin(2_prt*slice_ds*omega);
-            amrex::ParticleReal const term3 = -2_prt*q2*p2*w*cosh(2_prt*slice_ds*omega);
-            amrex::ParticleReal const term4 = -2_prt*q1*p1*w*cos(2_prt*slice_ds*omega);
+            amrex::ParticleReal const term1 = -(std::pow(p2,2) + std::pow(q2,2) * std::pow(w,2))*std::sinh(2_prt*slice_ds*omega);
+            amrex::ParticleReal const term2 = -(std::pow(p1,2)-pow(q1,2) * std::pow(w,2)) * std::sin(2_prt*slice_ds*omega);
+            amrex::ParticleReal const term3 = -2_prt*q2*p2*w*std::cosh(2_prt*slice_ds*omega);
+            amrex::ParticleReal const term4 = -2_prt*q1*p1*w * std::cos(2_prt*slice_ds*omega);
             amrex::ParticleReal const term5 = 2_prt*omega*(q1*p1*delta1 + q2*p2*delta1
-                                        -(pow(p1,2)+pow(p2,2))*slice_ds - (pow(q1,2)-pow(q2,2))*pow(w,2)*slice_ds);
-            t = t0 + (-1_prt+bet*pt)/(8_prt*bet*pow(delta1,3)*omega)
+                                        -(std::pow(p1,2) + std::pow(p2,2))*slice_ds - (std::pow(q1,2)-pow(q2,2)) * std::pow(w,2)*slice_ds);
+            t = t0 + (-1_prt+bet*pt)/(8_prt*bet * std::pow(delta1,3)*omega)
                      *(term1+term2+term3+term4+term5);
 
             // ptout = pt;
@@ -213,7 +213,7 @@ namespace impactx
             amrex::ParticleReal const slice_ds = m_ds / nslice();
 
             // assign intermediate parameter
-            amrex::ParticleReal const step = slice_ds / sqrt(pow(pt,2)-1.0_prt);
+            amrex::ParticleReal const step = slice_ds /std::sqrt(std::pow(pt,2)-1.0_prt);
 
             // advance position and momentum (straight element)
             refpart.x = x + step*px;
diff --git a/src/particles/elements/ChrUniformAcc.H b/src/particles/elements/ChrUniformAcc.H
index 465b480c4..540221ddf 100644
--- a/src/particles/elements/ChrUniformAcc.H
+++ b/src/particles/elements/ChrUniformAcc.H
@@ -100,8 +100,8 @@ namespace impactx
             // access reference particle values (final, initial):
             amrex::ParticleReal const ptf_ref = refpart.pt;
             amrex::ParticleReal const pti_ref = ptf_ref + m_ez*slice_ds;
-            amrex::ParticleReal const bgf = sqrt(pow(ptf_ref, 2) - 1.0_prt);
-            amrex::ParticleReal const bgi = sqrt(pow(pti_ref, 2) - 1.0_prt);
+            amrex::ParticleReal const bgf = std::sqrt(std::pow(ptf_ref, 2) - 1.0_prt);
+            amrex::ParticleReal const bgi = std::sqrt(std::pow(pti_ref, 2) - 1.0_prt);
 
             // initial conversion from static to dynamic units:
             px = px*bgi;
@@ -111,17 +111,17 @@ namespace impactx
             // compute intermediate quantities related to acceleration
             amrex::ParticleReal const pti_tot = pti_ref + pt;
             amrex::ParticleReal const ptf_tot = ptf_ref + pt;
-            amrex::ParticleReal const pzi_tot = sqrt(pow(pti_tot,2)-1_prt);
-            amrex::ParticleReal const pzf_tot = sqrt(pow(ptf_tot,2)-1_prt);
-            amrex::ParticleReal const pzi_ref = sqrt(pow(pti_ref,2)-1_prt);
-            amrex::ParticleReal const pzf_ref = sqrt(pow(ptf_ref,2)-1_prt);
+            amrex::ParticleReal const pzi_tot = std::sqrt(std::pow(pti_tot,2)-1_prt);
+            amrex::ParticleReal const pzf_tot = std::sqrt(std::pow(ptf_tot,2)-1_prt);
+            amrex::ParticleReal const pzi_ref = std::sqrt(std::pow(pti_ref,2)-1_prt);
+            amrex::ParticleReal const pzf_ref = std::sqrt(std::pow(ptf_ref,2)-1_prt);
 
             amrex::ParticleReal const numer = -ptf_tot + pzf_tot;
             amrex::ParticleReal const denom = -pti_tot + pzi_tot;
 
             // compute focusing constant (1/m) and rotation angle (in rad)
             amrex::ParticleReal const alpha = m_bz/2.0_prt;
-            amrex::ParticleReal const theta = alpha/m_ez*log(numer/denom);
+            amrex::ParticleReal const theta = alpha/m_ez*std::log(numer/denom);
 
             // intialize output values
             amrex::ParticleReal xout = x;
@@ -132,15 +132,15 @@ namespace impactx
             amrex::ParticleReal ptout = pt;
 
             // advance positions and momenta using map for focusing
-            xout = cos(theta)*x + sin(theta)/alpha*px;
-            pxout = -alpha*sin(theta)*x + cos(theta)*px;
+            xout = std::cos(theta)*x + std::sin(theta)/alpha*px;
+            pxout = -alpha * std::sin(theta)*x + std::cos(theta)*px;
 
-            yout = cos(theta)*y + sin(theta)/alpha*py;
-            pyout = -alpha*sin(theta)*y + cos(theta)*py;
+            yout = std::cos(theta)*y + std::sin(theta)/alpha*py;
+            pyout = -alpha * std::sin(theta)*y + std::cos(theta)*py;
 
             // the correct symplectic update for t
             tout = t + (pzf_tot - pzf_ref - pzi_tot + pzi_ref)/m_ez;
-            tout = tout + (1_prt/pzi_tot - 1_prt/pzf_tot)*(pow(py-alpha*x,2)+pow(px+alpha*y,2))/(2_prt*m_ez);
+            tout = tout + (1_prt/pzi_tot - 1_prt/pzf_tot)*(std::pow(py-alpha*x,2) + std::pow(px+alpha*y,2))/(2_prt*m_ez);
             ptout = pt;
 
             // assign intermediate momenta
@@ -149,11 +149,11 @@ namespace impactx
             pt = ptout;
 
             // advance positions and momenta using map for rotation
-            x = cos(theta)*xout + sin(theta)*yout;
-            pxout = cos(theta)*px + sin(theta)*py;
+            x = std::cos(theta)*xout + std::sin(theta)*yout;
+            pxout = std::cos(theta)*px + std::sin(theta)*py;
 
-            y = -sin(theta)*xout + cos(theta)*yout;
-            pyout = -sin(theta)*px + cos(theta)*py;
+            y = -sin(theta)*xout + std::cos(theta)*yout;
+            pyout = -sin(theta)*px + std::cos(theta)*py;
 
             t = tout;
             ptout = pt;
@@ -196,14 +196,14 @@ namespace impactx
             amrex::ParticleReal const slice_ds = m_ds / nslice();
 
             // compute intial value of beta*gamma
-            amrex::ParticleReal const bgi = sqrt(pow(pt, 2) - 1.0_prt);
+            amrex::ParticleReal const bgi = std::sqrt(std::pow(pt, 2) - 1.0_prt);
 
             // advance pt (uniform acceleration)
             refpart.pt = pt - m_ez*slice_ds;
 
             // compute final value of beta*gamma
             amrex::ParticleReal const ptf = refpart.pt;
-            amrex::ParticleReal const bgf = sqrt(pow(ptf, 2) - 1.0_prt);
+            amrex::ParticleReal const bgf = std::sqrt(std::pow(ptf, 2) - 1.0_prt);
 
             // update t
             refpart.t = t + (bgf - bgi)/m_ez;
diff --git a/src/particles/elements/ConstF.H b/src/particles/elements/ConstF.H
index ae1009e08..c8af6217e 100644
--- a/src/particles/elements/ConstF.H
+++ b/src/particles/elements/ConstF.H
@@ -92,7 +92,7 @@ namespace impactx
 
             // access reference particle values to find beta*gamma^2
             amrex::ParticleReal const pt_ref = refpart.pt;
-            amrex::ParticleReal const betgam2 = pow(pt_ref, 2) - 1.0_prt;
+            amrex::ParticleReal const betgam2 = std::pow(pt_ref, 2) - 1.0_prt;
 
             // intialize output values
             amrex::ParticleReal xout = x;
@@ -106,14 +106,14 @@ namespace impactx
             amrex::ParticleReal const slice_ds = m_ds / nslice();
 
             // advance position and momentum
-            xout = cos(m_kx*slice_ds)*x + sin(m_kx*slice_ds)/m_kx*px;
-            pxout = -m_kx*sin(m_kx*slice_ds)*x + cos(m_kx*slice_ds)*px;
+            xout = std::cos(m_kx*slice_ds)*x + std::sin(m_kx*slice_ds)/m_kx*px;
+            pxout = -m_kx * std::sin(m_kx*slice_ds)*x + std::cos(m_kx*slice_ds)*px;
 
-            yout = cos(m_ky*slice_ds)*y + sin(m_ky*slice_ds)/m_ky*py;
-            pyout = -m_ky*sin(m_ky*slice_ds)*y + cos(m_ky*slice_ds)*py;
+            yout = std::cos(m_ky*slice_ds)*y + std::sin(m_ky*slice_ds)/m_ky*py;
+            pyout = -m_ky * std::sin(m_ky*slice_ds)*y + std::cos(m_ky*slice_ds)*py;
 
-            tout = cos(m_kt*slice_ds)*t + sin(m_kt*slice_ds)/(betgam2*m_kt)*pt;
-            ptout = -(m_kt*betgam2)*sin(m_kt*slice_ds)*t + cos(m_kt*slice_ds)*pt;
+            tout = std::cos(m_kt*slice_ds)*t + std::sin(m_kt*slice_ds)/(betgam2*m_kt)*pt;
+            ptout = -(m_kt*betgam2) * std::sin(m_kt*slice_ds)*t + std::cos(m_kt*slice_ds)*pt;
 
             // assign updated values
             x = xout;
@@ -151,7 +151,7 @@ namespace impactx
             amrex::ParticleReal const slice_ds = m_ds / nslice();
 
             // assign intermediate parameter
-            amrex::ParticleReal const step = slice_ds / sqrt(pow(pt, 2)-1.0_prt);
+            amrex::ParticleReal const step = slice_ds /std::sqrt(std::pow(pt, 2)-1.0_prt);
 
             // advance position and momentum (straight element)
             refpart.x = x + step*px;
diff --git a/src/particles/elements/DipEdge.H b/src/particles/elements/DipEdge.H
index 15dba97bf..ca0bc3b9e 100644
--- a/src/particles/elements/DipEdge.H
+++ b/src/particles/elements/DipEdge.H
@@ -96,13 +96,13 @@ namespace impactx
             shift_in(x, y, px, py);
 
             // edge focusing matrix elements (zero gap)
-            amrex::ParticleReal const R21 = tan(m_psi)/m_rc;
+            amrex::ParticleReal const R21 = std::tan(m_psi)/m_rc;
             amrex::ParticleReal R43 = -R21;
             amrex::ParticleReal vf = 0;
 
             // first-order effect of nonzero gap
-            vf = (1.0_prt + pow(sin(m_psi),2))/(pow(cos(m_psi),3));
-            vf *= m_g * m_K2/(pow(m_rc,2));
+            vf = (1.0_prt + std::pow(sin(m_psi),2))/(std::pow(cos(m_psi),3));
+            vf *= m_g * m_K2/(std::pow(m_rc,2));
             R43 += vf;
 
             // apply edge focusing
diff --git a/src/particles/elements/Drift.H b/src/particles/elements/Drift.H
index c04c6d230..d20f466c1 100644
--- a/src/particles/elements/Drift.H
+++ b/src/particles/elements/Drift.H
@@ -97,7 +97,7 @@ namespace impactx
 
             // access reference particle values to find beta*gamma^2
             amrex::ParticleReal const pt_ref = refpart.pt;
-            amrex::ParticleReal const betgam2 = pow(pt_ref, 2) - 1.0_prt;
+            amrex::ParticleReal const betgam2 = std::pow(pt_ref, 2) - 1.0_prt;
 
             // advance position and momentum (drift)
             xout = x + slice_ds * px;
@@ -143,7 +143,7 @@ namespace impactx
             amrex::ParticleReal const slice_ds = m_ds / nslice();
 
             // assign intermediate parameter
-            amrex::ParticleReal const step = slice_ds / sqrt(pow(pt,2)-1.0_prt);
+            amrex::ParticleReal const step = slice_ds /std::sqrt(std::pow(pt,2)-1.0_prt);
 
             // advance position and momentum (drift)
             refpart.x = x + step*px;
diff --git a/src/particles/elements/ExactDrift.H b/src/particles/elements/ExactDrift.H
index cf10c12b8..8e8da94d9 100644
--- a/src/particles/elements/ExactDrift.H
+++ b/src/particles/elements/ExactDrift.H
@@ -103,8 +103,8 @@ namespace impactx
             amrex::ParticleReal const betgam = refpart.beta_gamma();
 
             // compute the radical in the denominator (= pz):
-            amrex::ParticleReal const pzden = sqrt(pow(pt-1_prt/bet,2) -
-                                1_prt/pow(betgam,2) - pow(px,2) - pow(py,2));
+            amrex::ParticleReal const pzden = std::sqrt(std::pow(pt-1_prt/bet,2) -
+                                1_prt/pow(betgam,2) - std::pow(px,2) - std::pow(py,2));
 
             // advance position and momentum (exact drift)
             x = xout + slice_ds * px / pzden;
@@ -147,7 +147,7 @@ namespace impactx
             amrex::ParticleReal const slice_ds = m_ds / nslice();
 
             // assign intermediate parameter
-            amrex::ParticleReal const step = slice_ds / sqrt(pow(pt,2)-1.0_prt);
+            amrex::ParticleReal const step = slice_ds /std::sqrt(std::pow(pt,2)-1.0_prt);
 
             // advance position and momentum (drift)
             refpart.x = x + step*px;
diff --git a/src/particles/elements/ExactSbend.H b/src/particles/elements/ExactSbend.H
index dbef6bc58..9967ed3f2 100644
--- a/src/particles/elements/ExactSbend.H
+++ b/src/particles/elements/ExactSbend.H
@@ -131,8 +131,8 @@ namespace impactx
             amrex::ParticleReal ptout = pt;
 
             // assign intermediate quantities
-            amrex::ParticleReal const pperp = sqrt(pow(pt,2)-2.0_prt/bet*pt-pow(py,2)+1.0_prt);
-            amrex::ParticleReal const pzi = sqrt(pow(pperp,2)-pow(px,2));
+            amrex::ParticleReal const pperp = std::sqrt(std::pow(pt,2)-2.0_prt/bet*pt-pow(py,2)+1.0_prt);
+            amrex::ParticleReal const pzi = std::sqrt(std::pow(pperp,2)-pow(px,2));
             amrex::ParticleReal const rho = rc + xout;
             auto const [sin_phi, cos_phi] = amrex::Math::sincos(slice_phi);
 
@@ -142,8 +142,8 @@ namespace impactx
             ptout = pt;
 
             // angle of momentum rotation
-            amrex::ParticleReal const pzf = sqrt(pow(pperp,2)-pow(pxout,2));
-            amrex::ParticleReal const theta = slice_phi + asin(px/pperp) - asin(pxout/pperp);
+            amrex::ParticleReal const pzf = std::sqrt(std::pow(pperp,2)-pow(pxout,2));
+            amrex::ParticleReal const theta = slice_phi + std::asin(px/pperp) - std::asin(pxout/pperp);
 
             // update position coordinates
             x = -rc + rho*cos_phi + rc*(pzf + px*sin_phi - pzi*cos_phi);
diff --git a/src/particles/elements/Multipole.H b/src/particles/elements/Multipole.H
index 04791d482..c0e9fa693 100644
--- a/src/particles/elements/Multipole.H
+++ b/src/particles/elements/Multipole.H
@@ -98,7 +98,7 @@ namespace impactx
 
             // access reference particle values to find (beta*gamma)^2
             //amrex::ParticleReal const pt_ref = refpart.pt;
-            //amrex::ParticleReal const betgam2 = pow(pt_ref, 2) - 1.0_prt;
+            //amrex::ParticleReal const betgam2 = std::pow(pt_ref, 2) - 1.0_prt;
 
             // intialize output values
             amrex::ParticleReal xout = x;
diff --git a/src/particles/elements/NonlinearLens.H b/src/particles/elements/NonlinearLens.H
index 93940b067..c3624cbb6 100644
--- a/src/particles/elements/NonlinearLens.H
+++ b/src/particles/elements/NonlinearLens.H
@@ -95,7 +95,7 @@ namespace impactx
 
             // access reference particle values to find (beta*gamma)^2
             //amrex::ParticleReal const pt_ref = refpart.pt;
-            //amrex::ParticleReal const betgam2 = pow(pt_ref, 2) - 1.0_prt;
+            //amrex::ParticleReal const betgam2 = std::pow(pt_ref, 2) - 1.0_prt;
 
             // intialize output values
             amrex::ParticleReal xout = x;
@@ -111,18 +111,18 @@ namespace impactx
             Complex const re1(1.0_prt, 0.0_prt);
             Complex const im1(0.0_prt, 1.0_prt);
 
-            // compute croot = sqrt(1-zeta**2)
+            // compute croot = std::sqrt(1-zeta**2)
             Complex croot = amrex::pow(zeta, 2);
             croot = re1 - croot;
             croot = amrex::sqrt(croot);
 
             // compute carcsin = arcsin(zeta)
-            Complex carcsin = im1*zeta + croot;
-            carcsin = -im1*amrex::log(carcsin);
+            Complex carcsin = im1 * zeta + croot;
+            carcsin = -im1 * amrex::log(carcsin);
 
             // compute complex function F'(zeta)
-            Complex dF = zeta/amrex::pow(croot, 2);
-            dF = dF + carcsin/amrex::pow(croot,3);
+            Complex dF = zeta / amrex::pow(croot, 2);
+            dF = dF + carcsin / amrex::pow(croot, 3);
 
             // compute momentum kick
             amrex::ParticleReal const kick = -m_knll/m_cnll;
diff --git a/src/particles/elements/PRot.H b/src/particles/elements/PRot.H
index ed21e74be..d08e7c009 100644
--- a/src/particles/elements/PRot.H
+++ b/src/particles/elements/PRot.H
@@ -97,8 +97,8 @@ namespace impactx
             auto const [sin_theta, cos_theta] = amrex::Math::sincos(theta);
             auto const [sin_phi_in, cos_phi_in] = amrex::Math::sincos(m_phi_in);
 
-            amrex::ParticleReal const pz = sqrt(1.0_prt - 2.0_prt*pt/beta
-               + pow(pt,2) - pow(py,2) - pow(px + sin_phi_in,2));
+            amrex::ParticleReal const pz = std::sqrt(1.0_prt - 2.0_prt*pt/beta
+               + std::pow(pt,2) - std::pow(py,2) - std::pow(px + sin_phi_in,2));
             amrex::ParticleReal const pzf = pz*cos_theta - (px + sin_phi_in)*sin_theta;
 
             // advance position and momentum
diff --git a/src/particles/elements/RFCavity.H b/src/particles/elements/RFCavity.H
index 6db8c9a12..2a83a5d66 100644
--- a/src/particles/elements/RFCavity.H
+++ b/src/particles/elements/RFCavity.H
@@ -282,7 +282,7 @@ namespace RFCavityData
             amrex::ParticleReal const slice_ds = m_ds / nslice();
 
             // compute intial value of beta*gamma
-            amrex::ParticleReal const bgi = sqrt(pow(pt, 2) - 1.0_prt);
+            amrex::ParticleReal const bgi = std::sqrt(std::pow(pt, 2) - 1.0_prt);
 
             // call integrator to advance (t,pt)
             amrex::ParticleReal const zin = s - sedge;
@@ -298,7 +298,7 @@ namespace RFCavityData
             refpart.z = z + slice_ds*pz/bgi;
 
             // compute final value of beta*gamma
-            amrex::ParticleReal const bgf = sqrt(pow(ptf, 2) - 1.0_prt);
+            amrex::ParticleReal const bgf = std::sqrt(std::pow(ptf, 2) - 1.0_prt);
 
             // advance momentum (px,py,pz)
             refpart.px = px*bgf/bgi;
@@ -367,14 +367,14 @@ namespace RFCavityData
                efieldint = z*efield;
                for (int j=1; j < m_ncoef; ++j)
                {
-                 efield = efield + cos_data[j]*cos(j*2*pi*z/zlen) +
-                     sin_data[j]*sin(j*2*pi*z/zlen);
-                 efieldp = efieldp-j*2*pi*cos_data[j]*sin(j*2*pi*z/zlen)/zlen +
-                      j*2*pi*sin_data[j]*cos(j*2*pi*z/zlen)/zlen;
-                 efieldpp = efieldpp- pow(j*2*pi*cos_data[j]/zlen,2) *cos(j*2*pi*z/zlen) -
-                      pow(j*2*pi*sin_data[j]/zlen,2) *sin(j*2*pi*z/zlen);
-                 efieldint = efieldint + zlen*cos_data[j]*sin(j*2*pi*z/zlen)/(j*2*pi) -
-                      zlen*sin_data[j]*cos(j*2*pi*z/zlen)/(j*2*pi);
+                 efield = efield + cos_data[j] * std::cos(j*2*pi*z/zlen) +
+                     sin_data[j] * std::sin(j*2*pi*z/zlen);
+                 efieldp = efieldp-j*2*pi*cos_data[j] * std::sin(j*2*pi*z/zlen)/zlen +
+                      j*2*pi*sin_data[j] * std::cos(j*2*pi*z/zlen)/zlen;
+                 efieldpp = efieldpp- std::pow(j*2*pi*cos_data[j]/zlen,2)  * std::cos(j*2*pi*z/zlen) -
+                     std::pow(j*2*pi*sin_data[j]/zlen,2)  * std::sin(j*2*pi*z/zlen);
+                 efieldint = efieldint + zlen*cos_data[j] * std::sin(j*2*pi*z/zlen)/(j*2*pi) -
+                      zlen*sin_data[j] * std::cos(j*2*pi*z/zlen)/(j*2*pi);
                }
             }
             else  // endpoint of the RF, outsize zlen
@@ -382,7 +382,7 @@ namespace RFCavityData
                efieldint = std::copysign(z, z)*zmid*0.5_prt*cos_data[0];;
                for (int j=1; j < m_ncoef; ++j)
                {
-                 efieldint = efieldint - zlen*sin_data[j]*cos(j*pi)/(j*2*pi);
+                 efieldint = efieldint - zlen*sin_data[j] * std::cos(j*pi)/(j*2*pi);
                }
             }
             return std::make_tuple(efield, efieldp, efieldint);
@@ -408,7 +408,7 @@ namespace RFCavityData
             amrex::ParticleReal const pt = refpart.pt;
 
             if (pt < -1.0_prt) {
-                refpart.t = t + tau/sqrt(1.0_prt - pow(pt, -2));
+                refpart.t = t + tau/std::sqrt(1.0_prt - std::pow(pt, -2));
                 refpart.pt = pt;
             }
             else {
@@ -459,17 +459,17 @@ namespace RFCavityData
             // push the linear map equations
             amrex::Array2D<amrex::ParticleReal, 1, 6, 1, 6> const R = refpart.map;
             amrex::ParticleReal const s = tau/refpart.beta_gamma();
-            amrex::ParticleReal const L = E0*ezp*sin(k*t+phi)/(2.0_prt*k);
+            amrex::ParticleReal const L = E0*ezp * std::sin(k*t+phi)/(2.0_prt*k);
 
             refpart.map(1,1) = (1.0_prt-s*L)*R(1,1) + s*R(2,1);
             refpart.map(1,2) = (1.0_prt-s*L)*R(1,2) + s*R(2,2);
-            refpart.map(2,1) = -s*pow(L,2)*R(1,1) + (1.0_prt+s*L)*R(2,1);
-            refpart.map(2,2) = -s*pow(L,2)*R(1,2) + (1.0_prt+s*L)*R(2,2);
+            refpart.map(2,1) = -s * std::pow(L,2)*R(1,1) + (1.0_prt+s*L)*R(2,1);
+            refpart.map(2,2) = -s * std::pow(L,2)*R(1,2) + (1.0_prt+s*L)*R(2,2);
 
             refpart.map(3,3) = (1.0_prt-s*L)*R(3,3) + s*R(4,3);
             refpart.map(3,4) = (1.0_prt-s*L)*R(3,4) + s*R(4,4);
-            refpart.map(4,3) = -s*pow(L,2)*R(3,3) + (1.0_prt+s*L)*R(4,3);
-            refpart.map(4,4) = -s*pow(L,2)*R(3,4) + (1.0_prt+s*L)*R(4,4);
+            refpart.map(4,3) = -s * std::pow(L,2)*R(3,3) + (1.0_prt+s*L)*R(4,3);
+            refpart.map(4,4) = -s * std::pow(L,2)*R(3,4) + (1.0_prt+s*L)*R(4,4);
         }
 
         /** This pushes the reference particle and the linear map matrix
@@ -506,12 +506,12 @@ namespace RFCavityData
             amrex::ignore_unused(ezf);
 
             refpart.t = t;
-            refpart.pt = pt - E0*(ezintf-ezint)*cos(k*t+phi);
+            refpart.pt = pt - E0*(ezintf-ezint) * std::cos(k*t+phi);
 
             // push the linear map equations
             amrex::Array2D<amrex::ParticleReal, 1, 6, 1, 6> const R = refpart.map;
-            amrex::ParticleReal const M = E0*(ezintf-ezint)*k*sin(k*t+phi);
-            amrex::ParticleReal const L = E0*(ezpf-ezp)*sin(k*t+phi)/(2.0_prt*k)+M/2.0_prt;
+            amrex::ParticleReal const M = E0*(ezintf-ezint)*k * std::sin(k*t+phi);
+            amrex::ParticleReal const L = E0*(ezpf-ezp) * std::sin(k*t+phi)/(2.0_prt*k)+M/2.0_prt;
 
             refpart.map(2,1) = L*R(1,1) + R(2,1);
             refpart.map(2,2) = L*R(1,2) + R(2,2);
diff --git a/src/particles/elements/Sbend.H b/src/particles/elements/Sbend.H
index 605be632b..e4e5c75be 100644
--- a/src/particles/elements/Sbend.H
+++ b/src/particles/elements/Sbend.H
@@ -103,8 +103,8 @@ namespace impactx
 
             // access reference particle values to find beta*gamma^2
             amrex::ParticleReal const pt_ref = refpart.pt;
-            amrex::ParticleReal const betgam2 = pow(pt_ref, 2) - 1.0_prt;
-            amrex::ParticleReal const bet = sqrt(betgam2/(1.0_prt + betgam2));
+            amrex::ParticleReal const betgam2 = std::pow(pt_ref, 2) - 1.0_prt;
+            amrex::ParticleReal const bet = std::sqrt(betgam2/(1.0_prt + betgam2));
 
             // calculate expensive terms once
             amrex::ParticleReal const theta = slice_ds/m_rc;
@@ -162,7 +162,7 @@ namespace impactx
 
             // assign intermediate parameter
             amrex::ParticleReal const theta = slice_ds/m_rc;
-            amrex::ParticleReal const B = sqrt(pow(pt,2)-1.0_prt)/m_rc;
+            amrex::ParticleReal const B = std::sqrt(std::pow(pt,2)-1.0_prt)/m_rc;
 
             // calculate expensive terms once
             auto const [sin_theta, cos_theta] = amrex::Math::sincos(theta);
diff --git a/src/particles/elements/ShortRF.H b/src/particles/elements/ShortRF.H
index b7fa316f4..543a66e27 100644
--- a/src/particles/elements/ShortRF.H
+++ b/src/particles/elements/ShortRF.H
@@ -99,9 +99,9 @@ namespace impactx
 
             // access reference particle values (final, initial):
             amrex::ParticleReal const ptf_ref = refpart.pt;
-            amrex::ParticleReal const pti_ref = ptf_ref + m_V*cos(phi);
-            amrex::ParticleReal const bgf = sqrt(pow(ptf_ref, 2) - 1.0_prt);
-            amrex::ParticleReal const bgi = sqrt(pow(pti_ref, 2) - 1.0_prt);
+            amrex::ParticleReal const pti_ref = ptf_ref + m_V * std::cos(phi);
+            amrex::ParticleReal const bgf = std::sqrt(std::pow(ptf_ref, 2) - 1.0_prt);
+            amrex::ParticleReal const bgi = std::sqrt(std::pow(pti_ref, 2) - 1.0_prt);
 
             // initial conversion from static to dynamic units:
             px = px*bgi;
@@ -124,7 +124,7 @@ namespace impactx
             pyout = py;
 
             // tout = t;
-            ptout = pt - m_V*cos(k*t + phi) + m_V*cos(phi);
+            ptout = pt - m_V * std::cos(k*t + phi) + m_V * std::cos(phi);
 
             // assign updated values
             x = xout;
@@ -170,14 +170,14 @@ namespace impactx
             amrex::ParticleReal const phi = m_phase*(pi/180.0_prt);
 
             // compute intial value of beta*gamma
-            amrex::ParticleReal const bgi = sqrt(pow(pt, 2) - 1.0_prt);
+            amrex::ParticleReal const bgi = std::sqrt(std::pow(pt, 2) - 1.0_prt);
 
             // advance pt
-            refpart.pt = pt - m_V*cos(phi);
+            refpart.pt = pt - m_V * std::cos(phi);
 
             // compute final value of beta*gamma
             amrex::ParticleReal const ptf = refpart.pt;
-            amrex::ParticleReal const bgf = sqrt(pow(ptf, 2) - 1.0_prt);
+            amrex::ParticleReal const bgf = std::sqrt(std::pow(ptf, 2) - 1.0_prt);
 
             // advance position (x,y,z,t)
             refpart.x = x;
diff --git a/src/particles/elements/SoftQuad.H b/src/particles/elements/SoftQuad.H
index 89a6c4755..cfd50eb44 100644
--- a/src/particles/elements/SoftQuad.H
+++ b/src/particles/elements/SoftQuad.H
@@ -285,7 +285,7 @@ namespace SoftQuadrupoleData
             amrex::ParticleReal const slice_ds = m_ds / nslice();
 
             // compute intial value of beta*gamma
-            amrex::ParticleReal const bgi = sqrt(pow(pt, 2) - 1.0_prt);
+            amrex::ParticleReal const bgi = std::sqrt(std::pow(pt, 2) - 1.0_prt);
 
             // call integrator to advance (t,pt)
             amrex::ParticleReal const zin = s - sedge;
@@ -312,7 +312,7 @@ namespace SoftQuadrupoleData
             refpart.z = z + slice_ds*pz/bgi;
 
             // compute final value of beta*gamma
-            amrex::ParticleReal const bgf = sqrt(pow(ptf, 2) - 1.0_prt);
+            amrex::ParticleReal const bgf = std::sqrt(std::pow(ptf, 2) - 1.0_prt);
 
             // advance momentum (px,py,pz)
             refpart.px = px*bgf/bgi;
@@ -362,12 +362,12 @@ namespace SoftQuadrupoleData
                bfieldint = z*bfield;
                for (int j=1; j < m_ncoef; ++j)
                {
-                 bfield = bfield + cos_data[j] * cos(j * 2 * pi * z / zlen) +
-                         sin_data[j] * sin(j * 2 * pi * z / zlen);
-                 bfieldp = bfieldp - j * 2 * pi * cos_data[j] * sin(j * 2 * pi * z / zlen) / zlen +
-                           j * 2 * pi * sin_data[j] * cos(j * 2 * pi * z / zlen) / zlen;
-                 bfieldint = bfieldint + zlen * cos_data[j] * sin(j * 2 * pi * z / zlen) / (j * 2 * pi) -
-                             zlen * sin_data[j] * cos(j * 2 * pi * z / zlen) / (j * 2 * pi);
+                 bfield = bfield + cos_data[j] *std::cos(j * 2 * pi * z / zlen) +
+                         sin_data[j] *std::sin(j * 2 * pi * z / zlen);
+                 bfieldp = bfieldp - j * 2 * pi * cos_data[j] *std::sin(j * 2 * pi * z / zlen) / zlen +
+                           j * 2 * pi * sin_data[j] *std::cos(j * 2 * pi * z / zlen) / zlen;
+                 bfieldint = bfieldint + zlen * cos_data[j] *std::sin(j * 2 * pi * z / zlen) / (j * 2 * pi) -
+                             zlen * sin_data[j] *std::cos(j * 2 * pi * z / zlen) / (j * 2 * pi);
                }
             }
             return std::make_tuple(bfield, bfieldp, bfieldint);
@@ -394,7 +394,7 @@ namespace SoftQuadrupoleData
             amrex::ParticleReal const z = zeval;
 
             if (pt < -1.0_prt) {
-                refpart.t = t + tau/sqrt(1.0_prt - pow(pt, -2));
+                refpart.t = t + tau/std::sqrt(1.0_prt - std::pow(pt, -2));
                 refpart.pt = pt;
             }
             else {
diff --git a/src/particles/elements/SoftSol.H b/src/particles/elements/SoftSol.H
index e3f79efc2..0151fc9c6 100644
--- a/src/particles/elements/SoftSol.H
+++ b/src/particles/elements/SoftSol.H
@@ -296,7 +296,7 @@ namespace SoftSolenoidData
             amrex::ParticleReal const slice_ds = m_ds / nslice();
 
             // compute intial value of beta*gamma
-            amrex::ParticleReal const bgi = sqrt(pow(pt, 2) - 1.0_prt);
+            amrex::ParticleReal const bgi = std::sqrt(std::pow(pt, 2) - 1.0_prt);
 
             // call integrator to advance (t,pt)
             amrex::ParticleReal const zin = s - sedge;
@@ -321,7 +321,7 @@ namespace SoftSolenoidData
             refpart.z = z + slice_ds*pz/bgi;
 
             // compute final value of beta*gamma
-            amrex::ParticleReal const bgf = sqrt(pow(ptf, 2) - 1.0_prt);
+            amrex::ParticleReal const bgf = std::sqrt(std::pow(ptf, 2) - 1.0_prt);
 
             // advance momentum (px,py,pz)
             refpart.px = px*bgf/bgi;
@@ -371,12 +371,12 @@ namespace SoftSolenoidData
                bfieldint = z*bfield;
                for (int j=1; j < m_ncoef; ++j)
                {
-                 bfield = bfield + cos_data[j]*cos(j*2*pi*z/zlen) +
-                     sin_data[j]*sin(j*2*pi*z/zlen);
-                 bfieldp = bfieldp-j*2*pi*cos_data[j]*sin(j*2*pi*z/zlen)/zlen +
-                      j*2*pi*sin_data[j]*cos(j*2*pi*z/zlen)/zlen;
-                 bfieldint = bfieldint + zlen*cos_data[j]*sin(j*2*pi*z/zlen)/(j*2*pi) -
-                      zlen*sin_data[j]*cos(j*2*pi*z/zlen)/(j*2*pi);
+                 bfield = bfield + cos_data[j] * std::cos(j*2*pi*z/zlen) +
+                     sin_data[j] * std::sin(j*2*pi*z/zlen);
+                 bfieldp = bfieldp-j*2*pi*cos_data[j] * std::sin(j*2*pi*z/zlen)/zlen +
+                      j*2*pi*sin_data[j] * std::cos(j*2*pi*z/zlen)/zlen;
+                 bfieldint = bfieldint + zlen*cos_data[j] * std::sin(j*2*pi*z/zlen)/(j*2*pi) -
+                      zlen*sin_data[j] * std::cos(j*2*pi*z/zlen)/(j*2*pi);
                }
             }
             return std::make_tuple(bfield, bfieldp, bfieldint);
@@ -403,7 +403,7 @@ namespace SoftSolenoidData
             amrex::ParticleReal const z = zeval;
 
             if (pt < -1.0_prt) {
-                refpart.t = t + tau/sqrt(1.0_prt - pow(pt, -2));
+                refpart.t = t + tau/std::sqrt(1.0_prt - std::pow(pt, -2));
                 refpart.pt = pt;
             }
             else {
@@ -466,7 +466,7 @@ namespace SoftSolenoidData
             // push the linear map equations
             amrex::Array2D<amrex::ParticleReal, 1, 6, 1, 6> const R = refpart.map;
             amrex::ParticleReal const alpha = B0*bz/2.0_prt;
-            amrex::ParticleReal const alpha2 = pow(alpha,2);
+            amrex::ParticleReal const alpha2 = std::pow(alpha,2);
 
             refpart.map(2,1) = R(2,1) - tau*alpha2*R(1,1);
             refpart.map(2,2) = R(2,2) - tau*alpha2*R(1,2);
@@ -515,8 +515,8 @@ namespace SoftSolenoidData
             // push the linear map equations
             amrex::Array2D<amrex::ParticleReal, 1, 6, 1, 6> const R = refpart.map;
             amrex::ParticleReal const theta = tau*B0*bz/2.0_prt;
-            amrex::ParticleReal const cs = cos(theta);
-            amrex::ParticleReal const sn = sin(theta);
+            amrex::ParticleReal const cs = std::cos(theta);
+            amrex::ParticleReal const sn = std::sin(theta);
 
             refpart.map(1,1) = R(1,1)*cs + R(3,1)*sn;
             refpart.map(1,2) = R(1,2)*cs + R(3,2)*sn;
diff --git a/src/particles/elements/Sol.H b/src/particles/elements/Sol.H
index 3753aafff..88d8e9be4 100644
--- a/src/particles/elements/Sol.H
+++ b/src/particles/elements/Sol.H
@@ -94,7 +94,7 @@ namespace impactx
 
             // access reference particle values to find beta*gamma^2
             amrex::ParticleReal const pt_ref = refpart.pt;
-            amrex::ParticleReal const betgam2 = pow(pt_ref, 2) - 1.0_prt;
+            amrex::ParticleReal const betgam2 = std::pow(pt_ref, 2) - 1.0_prt;
 
             // compute phase advance per unit length (in rad/m) and
             // rotation angle (in rad)
@@ -168,7 +168,7 @@ namespace impactx
             amrex::ParticleReal const slice_ds = m_ds / nslice();
 
             // assign intermediate parameter
-            amrex::ParticleReal const step = slice_ds / sqrt(pow(pt,2)-1.0_prt);
+            amrex::ParticleReal const step = slice_ds /std::sqrt(std::pow(pt,2)-1.0_prt);
 
             // advance position and momentum (straight element)
             refpart.x = x + step*px;
diff --git a/src/particles/elements/ThinDipole.H b/src/particles/elements/ThinDipole.H
index ecfd78e0d..5bd1e6633 100644
--- a/src/particles/elements/ThinDipole.H
+++ b/src/particles/elements/ThinDipole.H
@@ -104,7 +104,7 @@ namespace impactx
             amrex::ParticleReal ptout = pt;
 
             // compute the function expressing dp/p in terms of pt (labeled f in Ripken etc.)
-            amrex::ParticleReal f = -1.0_prt + sqrt(1.0_prt - 2.0_prt*pt/beta_ref + pow(pt,2));
+            amrex::ParticleReal f = -1.0_prt + std::sqrt(1.0_prt - 2.0_prt*pt/beta_ref + std::pow(pt,2));
             amrex::ParticleReal fprime = (1.0_prt - beta_ref*pt)/(beta_ref*(1.0_prt + f));
 
             // compute the effective (equivalent) arc length and curvature
@@ -113,7 +113,7 @@ namespace impactx
 
             // advance position and momentum
             x = xout;
-            pxout = px - pow(kx,2) * ds * xout + kx * ds * f; //eq (3.2b)
+            pxout = px - std::pow(kx,2) * ds * xout + kx * ds * f; //eq (3.2b)
 
             y = yout;
             pyout = py;
diff --git a/src/particles/integrators/Integrators.H b/src/particles/integrators/Integrators.H
index 636f04384..1bd770baf 100644
--- a/src/particles/integrators/Integrators.H
+++ b/src/particles/integrators/Integrators.H
@@ -137,7 +137,7 @@ namespace impactx::integrators
 
         // initialize numerical integration parameters
         amrex::ParticleReal const dz = (zout-zin)/nsteps;
-        amrex::ParticleReal const alpha = 1.0_prt - pow(2.0_prt,1.0/3.0);
+        amrex::ParticleReal const alpha = 1.0_prt - std::pow(2.0_prt,1.0/3.0);
         amrex::ParticleReal const tau2 = dz/(1.0_prt + alpha);
         amrex::ParticleReal const tau1 = tau2/2.0_prt;
         amrex::ParticleReal const tau3 = alpha*tau1;
diff --git a/src/particles/transformation/CoordinateTransformation.cpp b/src/particles/transformation/CoordinateTransformation.cpp
index cbf105533..299c95972 100644
--- a/src/particles/transformation/CoordinateTransformation.cpp
+++ b/src/particles/transformation/CoordinateTransformation.cpp
@@ -62,7 +62,7 @@ namespace impactx::transformation
                     amrex::ParticleReal *const AMREX_RESTRICT part_pz = soa_real[RealSoA::pz].dataPtr();
 
                     // Design value of pz/mc = beta*gamma
-                    amrex::ParticleReal const pzd = sqrt(pow(pd, 2) - 1.0);
+                    amrex::ParticleReal const pzd = std::sqrt(std::pow(pd, 2) - 1.0);
 
                     ToFixedS const to_s(pzd);
                     amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(long i) {
diff --git a/src/particles/transformation/ToFixedS.H b/src/particles/transformation/ToFixedS.H
index 10cb1b297..9dde2c3c0 100644
--- a/src/particles/transformation/ToFixedS.H
+++ b/src/particles/transformation/ToFixedS.H
@@ -59,7 +59,7 @@ namespace impactx::transformation
             using namespace amrex::literals;
 
             // compute value of reference ptd = -gamma
-            amrex::ParticleReal const argd = 1.0_prt + pow(m_pzd, 2);
+            amrex::ParticleReal const argd = 1.0_prt + std::pow(m_pzd, 2);
             AMREX_ASSERT_WITH_MESSAGE(argd > 0.0_prt, "invalid ptd arg (<=0)");
             amrex::ParticleReal const ptdf = argd > 0.0_prt ? -sqrt(argd) : -1.0_prt;
 
@@ -70,7 +70,7 @@ namespace impactx::transformation
             pz = pz * m_pzd;
 
             // compute value of particle pt = -gamma
-            amrex::ParticleReal const arg = 1.0_prt + pow(m_pzd + pz, 2) + pow(px, 2) + pow(py, 2);
+            amrex::ParticleReal const arg = 1.0_prt + std::pow(m_pzd + pz, 2) + std::pow(px, 2) + std::pow(py, 2);
             AMREX_ASSERT_WITH_MESSAGE(arg > 0.0_prt, "invalid pt arg (<=0)");
             amrex::ParticleReal const ptf = arg > 0.0_prt ? -sqrt(arg) : -1.0_prt;
 
diff --git a/src/particles/transformation/ToFixedT.H b/src/particles/transformation/ToFixedT.H
index 577436ea3..af90cc132 100644
--- a/src/particles/transformation/ToFixedT.H
+++ b/src/particles/transformation/ToFixedT.H
@@ -62,9 +62,9 @@ namespace impactx::transformation
             constexpr amrex::ParticleReal tol = 1.0e-8_prt;
 
             // compute value of reference pzd = beta*gamma
-            amrex::ParticleReal const argd = -1.0_prt + pow(m_ptd, 2);
+            amrex::ParticleReal const argd = -1.0_prt + std::pow(m_ptd, 2);
             // AMREX_ASSERT_WITH_MESSAGE(argd > 0.0_prt, "invalid pzd arg (<=0)");
-            amrex::ParticleReal const pzdf = argd > 0.0_prt ? sqrt(argd) : tol;
+            amrex::ParticleReal const pzdf = argd > 0.0_prt ? std::sqrt(argd) : tol;
 
             // transform momenta to dynamic units (eg, so that momenta are
             // normalized by mc):
@@ -73,9 +73,9 @@ namespace impactx::transformation
             pt = pt * pzdf;
 
             // compute value of particle pz = beta*gamma
-            amrex::ParticleReal const arg = -1.0_prt + pow(m_ptd+pt, 2) - pow(px, 2) - pow(py, 2);
+            amrex::ParticleReal const arg = -1.0_prt + std::pow(m_ptd+pt, 2) - std::pow(px, 2) - std::pow(py, 2);
             // AMREX_ASSERT_WITH_MESSAGE(arg > 0.0_prt, "invalid pz arg (<=0)");
-            amrex::ParticleReal const pzf = arg > 0.0_prt ? sqrt(arg) : tol;
+            amrex::ParticleReal const pzf = arg > 0.0_prt ? std::sqrt(arg) : tol;
 
             // transform position and momentum (from fixed s to fixed t)
             x = x + px*t/(m_ptd+pt);
diff --git a/src/particles/wakefields/WakeConvolution.cpp b/src/particles/wakefields/WakeConvolution.cpp
index 680ed6d08..c18cf4fc7 100644
--- a/src/particles/wakefields/WakeConvolution.cpp
+++ b/src/particles/wakefields/WakeConvolution.cpp
@@ -42,7 +42,7 @@ namespace impactx::particles::wakefields
         using namespace amrex::literals;
 
         amrex::Real const s0 = (0.169_rt * std::pow(a, 1.79_rt) * std::pow(g, 0.38_rt)) / std::pow(L, 1.17_rt);
-        amrex::Real const term = std::sqrt(std::abs(s) / s0) * std::exp(-std::sqrt(std::abs(s) / s0));
+        amrex::Real const term = std::sqrt(std::abs(s) / s0) * std::exp(- std::sqrt(std::abs(s) / s0));
         return (4_rt * impactx::particles::wakefields::Z0 * ablastr::constant::SI::c * s0 * unit_step(s)) / (amrex::Real(M_PI) * std::pow(a, 4)) * term;
     }
 
@@ -56,7 +56,7 @@ namespace impactx::particles::wakefields
         using namespace amrex::literals;
 
         amrex::Real const s00 = g * std::pow((a / (alpha(g / L) * L)), 2) / 8.0_rt;
-        return (impactx::particles::wakefields::Z0 * ablastr::constant::SI::c * unit_step(s) * std::exp(-std::sqrt(std::abs(s) / s00))) / (amrex::Real(M_PI) * std::pow(a, 2));
+        return (impactx::particles::wakefields::Z0 * ablastr::constant::SI::c * unit_step(s) * std::exp(- std::sqrt(std::abs(s) / s00))) / (amrex::Real(M_PI) * std::pow(a, 2));
     }
 
     amrex::Gpu::DeviceVector<amrex::Real>