diff --git a/Jlqcd.hpp b/Jlqcd.hpp
index 9fac006..7b53c37 100644
--- a/Jlqcd.hpp
+++ b/Jlqcd.hpp
@@ -44,30 +44,34 @@ struct Jlqcd
     unsigned int L, T, Ls;
   };
 
-  inline static constexpr EnsembleParameters fUd3SaPar{0.003, 0.015, 1., 2., 64, 128, 8};
-  inline static constexpr EnsembleParameters mUd3SaPar{0.0042, 0.025, 1., 2., 48, 96, 8};
-  inline static constexpr EnsembleParameters mUd3SbPar{0.0042, 0.018, 1., 2., 48, 96, 8};
-  inline static constexpr EnsembleParameters cUd2SaPar{0.0035, 0.040, 1., 2., 32, 64, 12};
-  inline static constexpr EnsembleParameters cUd2SaLPar{0.0035, 0.040, 1., 2., 48, 96, 12};
-  inline static constexpr EnsembleParameters cUd3SaPar{0.007, 0.040, 1., 2., 32, 64, 12};
-  inline static constexpr EnsembleParameters cUd3SbPar{0.007, 0.030, 1., 2., 32, 64, 12};
-  inline static constexpr EnsembleParameters testPar{0.010, 0.030, 1., 2., 8, 16, 8};
+  inline static constexpr EnsembleParameters fUd3SaPar {0.003,  0.015, 1., 2., 64, 128, 8};
+  inline static constexpr EnsembleParameters mUd3SaPar {0.0042, 0.025, 1., 2., 48, 96,  8};
+  inline static constexpr EnsembleParameters mUd3SbPar {0.0042, 0.018, 1., 2., 48, 96,  8};
+  inline static constexpr EnsembleParameters cUd2SaPar {0.0035, 0.040, 1., 2., 32, 64,  12};
+  inline static constexpr EnsembleParameters cUd2SaLPar{0.0035, 0.040, 1., 2., 48, 96,  12};
+  inline static constexpr EnsembleParameters cUd3SaPar {0.007,  0.040, 1., 2., 32, 64,  12};
+  inline static constexpr EnsembleParameters cUd3SbPar {0.007,  0.030, 1., 2., 32, 64,  12};
+  inline static constexpr EnsembleParameters testPar   {0.010,  0.030, 1., 2.,  8, 16,   8};
 
-  static inline const EnsembleParameters &ensemble(const std::string &name);
+  static inline const EnsembleParameters &ensemble(const std::string& name);
 
-  static inline void addMixedPrecisionSolver(Application &app, const EnsembleParameters &par,
+  static inline void addMixedPrecisionSolver(Application &app,
+                                             const EnsembleParameters &par,
                                              const std::string &solverName,
-                                             const std::string &gaugeName, const double mass,
+                                             const std::string &gaugeName,
+                                             const double mass,
                                              const std::string &twist,
                                              const std::string &boundary = "1 1 1 -1",
                                              const double residual = 1.0e-8);
 
   static inline void addNoFailSolver(Application &app, const EnsembleParameters &par,
-                                     const std::string &solverName, const std::string &gaugeName,
-                                     const double mass, const std::string &twist,
+                                     const std::string &solverName,
+                                     const std::string &gaugeName, const double mass,
+                                     const std::string &twist,
                                      const std::string &boundary = "1 1 1 -1",
                                      const double residual = 1.0e-20);
 
+  // Light solvers (no deflation)
   static inline void addLightSolver(Application &app, const EnsembleParameters &par,
                                     const std::string &solverName, const std::string &gaugeName,
                                     const double mass, const std::string &twist,
@@ -80,52 +84,114 @@ struct Jlqcd
                                     const std::string &boundary = "1 1 1 -1",
                                     const double residual = 1.0e-8);
 
+  static inline void addfUd3SaLightSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addmUd3SaLightSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addmUd3SbLightSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd2SaLightSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd2SaLLightSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd3SaLightSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd3SbLightSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+
+  // Strange solvers
   static inline void addStrangeSolver(Application &app, const EnsembleParameters &par,
-                                      const std::string &solverName, const std::string &gaugeName,
-                                      const double mass, const std::string &twist,
+                                      const std::string &solverName,
+                                      const std::string &gaugeName, const double mass,
+                                      const std::string &twist,
                                       const std::string &boundary = "1 1 1 -1",
                                       const double residual = 1.0e-8);
 
   static inline void addStrangeSolver(Application &app, const EnsembleParameters &par,
-                                      const std::string &solverName, const std::string &gaugeName,
+                                      const std::string &solverName,
+                                      const std::string &gaugeName,
                                       const std::string &twist = "0. 0. 0. 0.",
                                       const std::string &boundary = "1 1 1 -1",
                                       const double residual = 1.0e-8);
 
+  static inline void addfUd3SaStrangeSolver(Application &app, const std::string &solverName,
+                                            const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                            const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addmUd3SaStrangeSolver(Application &app, const std::string &solverName,
+                                            const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                            const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addmUd3SbStrangeSolver(Application &app, const std::string &solverName,
+                                            const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                            const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd2SaStrangeSolver(Application &app, const std::string &solverName,
+                                            const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                            const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd2SaLStrangeSolver(Application &app, const std::string &solverName,
+                                            const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                            const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd3SaStrangeSolver(Application &app, const std::string &solverName,
+                                            const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                            const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+  static inline void addcUd3SbStrangeSolver(Application &app, const std::string &solverName,
+                                            const std::string &gaugeName, const std::string &twist = "0. 0. 0. 0.",
+                                            const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-8);
+
+  // Heavy solvers
   static inline void addHeavySolver(Application &app, const EnsembleParameters &par,
                                     const std::string &solverName, const std::string &gaugeName,
                                     const double mass, const std::string &twist,
                                     const std::string &boundary = "1 1 1 -1",
                                     const double residual = 1.0e-20);
+
+  static inline void addfUd3SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                          const double mass, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-20);
+  static inline void addmUd3SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                          const double mass, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-20);
+  static inline void addmUd3SbHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                          const double mass, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-20);
+  static inline void addcUd2SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                          const double mass, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-20);
+  static inline void addcUd2SaLHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                          const double mass, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-20);
+  static inline void addcUd3SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                          const double mass, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-20);
+  static inline void addcUd3SbHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                          const double mass, const std::string &twist = "0. 0. 0. 0.",
+                                          const std::string &boundary = "1 1 1 -1", const double residual = 1.0e-20);
+
 };
 
-inline const Jlqcd::EnsembleParameters &Jlqcd::ensemble(const std::string &name)
+inline const Jlqcd::EnsembleParameters &Jlqcd::ensemble(const std::string& name)
 {
-  if (name == "fUd3Sa")
-    return fUd3SaPar;
-  if (name == "mUd3Sa")
-    return mUd3SaPar;
-  if (name == "mUd3Sb")
-    return mUd3SbPar;
-  if (name == "cUd2Sa")
-    return cUd2SaPar;
-  if (name == "cUd2SaL")
-    return cUd2SaLPar;
-  if (name == "cUd3Sa")
-    return cUd3SaPar;
-  if (name == "cUd3Sb")
-    return cUd3SbPar;
-  if (name == "test")
-    return testPar;
+  if (name == "fUd3Sa" ) return fUd3SaPar;
+  if (name == "mUd3Sa" ) return mUd3SaPar;
+  if (name == "mUd3Sb" ) return mUd3SbPar;
+  if (name == "cUd2Sa" ) return cUd2SaPar;
+  if (name == "cUd2SaL") return cUd2SaLPar;
+  if (name == "cUd3Sa" ) return cUd3SaPar;
+  if (name == "cUd3Sb" ) return cUd3SbPar;
+  if (name == "test"   ) return testPar;
 
   throw std::invalid_argument("unknown JLQCD ensemble preset '" + std::string(name) + "'");
 }
 
 inline void Jlqcd::addMixedPrecisionSolver(Application &app, const EnsembleParameters &par,
-                                           const std::string &solverName,
-                                           const std::string &gaugeName, const double mass,
-                                           const std::string &twist, const std::string &boundary,
-                                           const double residual)
+                                           const std::string &solverName, const std::string &gaugeName,
+                                           const double mass, const std::string &twist,
+                                           const std::string &boundary, const double residual)
 {
   const std::string prefix = solverName;
 
@@ -182,6 +248,7 @@ inline void Jlqcd::addNoFailSolver(Application &app, const EnsembleParameters &p
   app.createModule<MSolver::RBPrecCGNoFail>(solverName, solverPar);
 }
 
+// Light solvers (no deflation)
 inline void Jlqcd::addLightSolver(Application &app, const EnsembleParameters &par,
                                   const std::string &solverName, const std::string &gaugeName,
                                   const double mass, const std::string &twist,
@@ -198,12 +265,62 @@ inline void Jlqcd::addLightSolver(Application &app, const EnsembleParameters &pa
   addLightSolver(app, par, solverName, gaugeName, par.ml, twist, boundary, residual);
 }
 
+inline void Jlqcd::addfUd3SaLightSolver(Application &app, const std::string &solverName,
+                                        const std::string &gaugeName, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+    addLightSolver(app, fUd3SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addmUd3SaLightSolver(Application &app, const std::string &solverName,
+                                        const std::string &gaugeName, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+    addLightSolver(app, mUd3SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addmUd3SbLightSolver(Application &app, const std::string &solverName,
+                                        const std::string &gaugeName, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+    addLightSolver(app, mUd3SbPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd2SaLightSolver(Application &app, const std::string &solverName,
+                                        const std::string &gaugeName, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+    addLightSolver(app, cUd2SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd2SaLLightSolver(Application &app, const std::string &solverName,
+                                        const std::string &gaugeName, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+    addLightSolver(app, cUd2SaLPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd3SaLightSolver(Application &app, const std::string &solverName,
+                                        const std::string &gaugeName, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+    addLightSolver(app, cUd3SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd3SbLightSolver(Application &app, const std::string &solverName,
+                                        const std::string &gaugeName, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+    addLightSolver(app, cUd3SbPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+// Strange solvers
 inline void Jlqcd::addStrangeSolver(Application &app, const EnsembleParameters &par,
                                     const std::string &solverName, const std::string &gaugeName,
                                     const double mass, const std::string &twist,
                                     const std::string &boundary, const double residual)
 {
-  addMixedPrecisionSolver(app, par, solverName, gaugeName, mass, twist, boundary, residual);
+    addMixedPrecisionSolver(app, par, solverName, gaugeName, mass, twist, boundary, residual);
 }
 
 inline void Jlqcd::addStrangeSolver(Application &app, const EnsembleParameters &par,
@@ -211,9 +328,59 @@ inline void Jlqcd::addStrangeSolver(Application &app, const EnsembleParameters &
                                     const std::string &twist, const std::string &boundary,
                                     const double residual)
 {
-  addStrangeSolver(app, par, solverName, gaugeName, par.ms, twist, boundary, residual);
+    addStrangeSolver(app, par, solverName, gaugeName, par.ms, twist, boundary, residual);
 }
 
+inline void Jlqcd::addfUd3SaStrangeSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist,
+                                          const std::string &boundary, const double residual)
+{
+    addStrangeSolver(app, fUd3SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addmUd3SaStrangeSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist,
+                                          const std::string &boundary, const double residual)
+{
+    addStrangeSolver(app, mUd3SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addmUd3SbStrangeSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist,
+                                          const std::string &boundary, const double residual)
+{
+    addStrangeSolver(app, mUd3SbPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd2SaStrangeSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist,
+                                          const std::string &boundary, const double residual)
+{
+    addStrangeSolver(app, cUd2SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd2SaLStrangeSolver(Application &app, const std::string &solverName,
+                                           const std::string &gaugeName, const std::string &twist,
+                                           const std::string &boundary, const double residual)
+{
+    addStrangeSolver(app, cUd2SaLPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd3SaStrangeSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist,
+                                          const std::string &boundary, const double residual)
+{
+    addStrangeSolver(app, cUd3SaPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+inline void Jlqcd::addcUd3SbStrangeSolver(Application &app, const std::string &solverName,
+                                          const std::string &gaugeName, const std::string &twist,
+                                          const std::string &boundary, const double residual)
+{
+    addStrangeSolver(app, cUd3SbPar, solverName, gaugeName, twist, boundary, residual);
+}
+
+// Heavy solvers
 inline void Jlqcd::addHeavySolver(Application &app, const EnsembleParameters &par,
                                   const std::string &solverName, const std::string &gaugeName,
                                   const double mass, const std::string &twist,
@@ -221,4 +388,48 @@ inline void Jlqcd::addHeavySolver(Application &app, const EnsembleParameters &pa
 {
   addNoFailSolver(app, par, solverName, gaugeName, mass, twist, boundary, residual);
 }
+
+inline void Jlqcd::addfUd3SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                        const double mass, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+  addHeavySolver(app, fUd3SaPar, solverName, gaugeName, mass, twist, boundary, residual);
+}
+inline void Jlqcd::addmUd3SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                        const double mass, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+  addHeavySolver(app, mUd3SaPar, solverName, gaugeName, mass, twist, boundary, residual);
+}
+inline void Jlqcd::addmUd3SbHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                        const double mass, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+  addHeavySolver(app, mUd3SbPar, solverName, gaugeName, mass, twist, boundary, residual);
+}
+inline void Jlqcd::addcUd2SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                        const double mass, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+  addHeavySolver(app, cUd2SaPar, solverName, gaugeName, mass, twist, boundary, residual);
+}
+inline void Jlqcd::addcUd2SaLHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                         const double mass, const std::string &twist,
+                                         const std::string &boundary, const double residual)
+{
+  addHeavySolver(app, cUd2SaLPar, solverName, gaugeName, mass, twist, boundary, residual);
+}
+inline void Jlqcd::addcUd3SaHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                        const double mass, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+  addHeavySolver(app, cUd3SaPar, solverName, gaugeName, mass, twist, boundary, residual);
+}
+inline void Jlqcd::addcUd3SbHeavySolver(Application &app, const std::string &solverName, const std::string &gaugeName, 
+                                        const double mass, const std::string &twist,
+                                        const std::string &boundary, const double residual)
+{
+  addHeavySolver(app, cUd3SbPar, solverName, gaugeName, mass, twist, boundary, residual);
+}
+
 } // namespace hadpresets