diff --git a/exec/multispec/AdvanceTimestepBousq.cpp b/exec/multispec/AdvanceTimestepBousq.cpp
index 8d9f3931..57aeff97 100644
--- a/exec/multispec/AdvanceTimestepBousq.cpp
+++ b/exec/multispec/AdvanceTimestepBousq.cpp
@@ -9,6 +9,7 @@
 
 #include <AMReX_ParallelDescriptor.H>
 #include <AMReX_MultiFabUtil.H>
+#include <AMReX_Print.H>
 
 
 // argv contains the name of the inputs file entered at the command line
@@ -33,6 +34,10 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
                           MultiFab& permittivity,
                           StochMassFlux& sMassFlux,
                           StochMomFlux& sMomFlux,
+                          MultiFab& phi_old,
+                          MultiFab& phi_new,
+                          MultiFab& phitot_old,
+                          MultiFab& phitot_new,
                           const Real& dt,
                           const Real& time,
                           const int& istep,
@@ -122,7 +127,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         }
     }
     
-    if (use_multiphase || use_flory_huggins) {
+    if (use_multiphase || use_flory_huggins || use_ice_nucleation) {
         for (int d=0; d<AMREX_SPACEDIM; ++d) {
             div_reversible_stress[d].define(convert(ba,nodal_flag_dir[d]), dmap, 1, 0);
         }
@@ -199,7 +204,8 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
     if (any_grav) {
         Abort("AdvanceTimestepBousq.cpp gravity not implemented");
     }
-
+    
+if (use_ice_nucleation ==0) {
     if (variance_coef_mass != 0.) {
         sMassFlux.fillMassStochastic();
     }
@@ -217,7 +223,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
     ComputeMassFluxdiv(rho_old,rhotot_old,Temp,diff_mass_fluxdiv,stoch_mass_fluxdiv,
                        diff_mass_flux,stoch_mass_flux,sMassFlux,0.5*dt,time,geom,weights_mass,
                        charge_old,grad_Epot_old,Epot,permittivity);
-    
+                       
     // here is a reasonable place to call something to compute in reversible stress term
     // in this case want to get divergence so it looks like a add to rhs for stokes solver
     if (use_multiphase) {
@@ -240,6 +246,18 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         }
 
     }
+ } else {
+
+    	  ComputeRhotot(phi_old,phitot_old);
+        // compute reversible stress tensor ---added term
+        ComputeDivFHReversibleStress(div_reversible_stress,phitot_old,phi_old,geom);
+
+        // add divergence of reversible stress to gmres_rhs_v
+        for (int d=0; d<AMREX_SPACEDIM; ++d) {
+            MultiFab::Saxpy(gmres_rhs_v[d],1.,div_reversible_stress[d],0,0,1,0);
+        }
+
+}
 
     if (use_charged_fluid) {
 
@@ -358,6 +376,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         eta_ed[d].mult(2.,0,1,0);
     }
 
+if (use_ice_nucleation ==0){
     //////////////////////////////////////////////
     // Step 2: compute reactions and mass fluxes at t^n
     //////////////////////////////////////////////
@@ -368,7 +387,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
 
     end if
     */
-
+    
     if (advection_type == 1 || advection_type == 2) {
 
       // bds advection
@@ -415,7 +434,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
       Abort("Invalid advection_type");
 
     }
-
+    
     //////////////////////////////////////////////
     /// Step 3: density prediction to t^{n+1/2}
     //////////////////////////////////////////////
@@ -440,7 +459,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         end if
         */
     }
-
+    
     // compute rhotot^{n+1/2} from rho^{n+1/2} in VALID REGION
     ComputeRhotot(rho_new,rhotot_new);
 
@@ -450,6 +469,17 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
     // average rho_i^{n+1/2} to faces
     AverageCCToFace(rho_new,rho_fc,0,nspecies,SPEC_BC_COMP,geom);
 
+} else {
+    // compute rhotot^{n+1/2} from rho^{n+1/2} in VALID REGION
+    ComputeRhotot(rho_old,rhotot_new);
+
+    // fill rho and rhotot ghost cells at t^{n+1/2}
+    FillRhoRhototGhost(rho_old,rhotot_new,geom);
+    
+    // average rho_i^{n+1/2} to faces
+    //AverageCCToFace(rho_old,rho_fc,0,nspecies,SPEC_BC_COMP,geom);
+}
+
     if (use_charged_fluid) {
         // compute total charge at t^{n+1/2}
         DotWithZ(rho_new,charge_new);
@@ -461,7 +491,11 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
     }
 
     // compute (eta,kappa)^{n+1/2}
+if (use_ice_nucleation ==0){
     ComputeEta(rho_new, rhotot_new, eta);
+} else {
+    ComputeEta(phi_old, phitot_old, eta);
+}
     if (AMREX_SPACEDIM == 2) {
         AverageCCToNode(eta,eta_ed[0],0,1,SPEC_BC_COMP,geom);
     }
@@ -469,6 +503,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         AverageCCToEdge(eta,eta_ed,0,1,SPEC_BC_COMP,geom);
     }
 
+if (use_ice_nucleation ==0){
     //////////////////////////////////////////////
     // Step 4: compute mass fluxes and reactions at t^{n+1/2}
     //////////////////////////////////////////////
@@ -525,7 +560,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
     } else if (use_flory_huggins ==1){
 
         // compute reversible stress tensor ---added term
-          ComputeDivFHReversibleStress(div_reversible_stress,rhotot_new,rho_new,geom);
+        ComputeDivFHReversibleStress(div_reversible_stress,rhotot_new,rho_new,geom);
 
     }
 
@@ -568,7 +603,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         MkAdvSFluxdiv(umac_old,rho_fc,adv_mass_fluxdiv,geom,0,nspecies,false);
         MkAdvSFluxdiv(umac    ,rho_fc,adv_mass_fluxdiv,geom,0,nspecies,true);
     }
-
+    
     //////////////////////////////////////////////
     // Step 5: density integration to t^{n+1}
     //////////////////////////////////////////////
@@ -606,6 +641,108 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
 
     // average rho^{n+1} to faces
     AverageCCToFace(rhotot_new,rhotot_fc_new,0,1,RHO_BC_COMP,geom);
+    	
+} else {
+    	AverageCCToFace(rhotot_old,rhotot_fc_new,0,1,RHO_BC_COMP,geom);
+   	
+   	MultiFab phi_prd(ba,dmap,1,2);
+   	
+	   FillRhoRhototGhost(phi_old,phitot_old,geom);
+	   
+	   // predict phi at t+1/2
+	   for (MFIter mfi(phi_prd); mfi.isValid(); ++mfi )
+	   {
+	       Box bx = mfi.tilebox();
+	       const Array4<Real>& phiOld 	= phi_old.array(mfi);
+	       const Array4<Real>& phiPrd 	= phi_prd.array(mfi);
+	      
+	       const amrex::Array4<amrex::Real>& u = umac[0].array(mfi);
+	       const amrex::Array4<amrex::Real>& v = umac[1].array(mfi);
+	       const amrex::Array4<amrex::Real>& u_old = umac_old[0].array(mfi);
+	       const amrex::Array4<amrex::Real>& v_old = umac_old[1].array(mfi);
+#if (AMREX_SPACEDIM == 3)
+	       const amrex::Array4<amrex::Real>& w = umac[2].array(mfi);
+	       const amrex::Array4<amrex::Real>& w_old = umac_old[2].array(mfi);
+#endif
+	       int n=0;
+	       amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+	       {
+	       	
+	           phiPrd(i,j,k) = phiOld(i,j,k,n)
+	           		+ dt/2 * Mobility * 2 * GradEnCoef *
+	               ( (phiOld(i+1,j,k,n) - 2.*phiOld(i,j,k,n) + phiOld(i-1,j,k,n)) / (dx[0]*dx[0])
+	                +(phiOld(i,j+1,k,n) - 2.*phiOld(i,j,k,n) + phiOld(i,j-1,k,n)) / (dx[1]*dx[1])
+	               )
+	               - dt/2 * ((u(i+1,j,k)+u_old(i+1,j,k))/2*(phiOld(i+1,j,k,n)+phiOld(i,j,k,n))/2-(u(i,j,k)+u_old(i,j,k))/2*(phiOld(i,j,k,n)+phiOld(i-1,j,k,n))/2)/dx[0]
+	               - dt/2 * ((v(i,j+1,k)+v_old(i,j+1,k))/2*(phiOld(i,j+1,k,n)+phiOld(i,j,k,n))/2-(v(i,j,k)+v_old(i,j,k))/2*(phiOld(i,j,k,n)+phiOld(i,j-1,k,n))/2)/dx[1]
+#if (AMREX_SPACEDIM == 3)
+	               + dt/2 * Mobility * 2 * GradEnCoef * (phiOld(i,j,k+1,n) - 2.*phiOld(i,j,k,n) + phiOld(i,j,k-1,n)) / (dx[2]*dx[2])
+	               - dt/2 * ((w(i,j,k+1)+w_old(i,j,k+1))/2*(phiOld(i,j,k+1,n)+phiOld(i,j,k,n))/2-(w(i,j,k)+w_old(i,j,k))/2*(phiOld(i,j,k,n)+phiOld(i,j,k-1,n))/2)/dx[2]
+               	+ variance_coef_mass*sqrt(2.0*k_B*T_init[0]*Mobility*dt/(dx[0]*dx[1]*dx[2]))*amrex::RandomNormal(0,1)/sqrt(2)
+#elif (AMREX_SPACEDIM == 2)
+               	+ variance_coef_mass*sqrt(2.0*k_B*T_init[0]*Mobility*dt/(dx[0]*dx[1]*cell_depth))*amrex::RandomNormal(0,1)/sqrt(2)
+#endif
+						- dt/2 * Mobility * EnScale*(2*phiOld(i,j,k,n)*(phiOld(i,j,k,n)-1)*(phiOld(i,j,k,n)-1)+2*phiOld(i,j,k,n)*phiOld(i,j,k,n)*(phiOld(i,j,k,n)-1) - PotWellDepr)
+	               ;
+	       });
+	   }
+	   
+	   phi_prd.FillBoundary(geom.periodicity());
+	   // advance phi to t+1
+	   for (MFIter mfi(phi_new); mfi.isValid(); ++mfi )
+	   {
+	       Box bx = mfi.tilebox();
+	       const Array4<Real>& phiNew 	= phi_new.array(mfi);
+	       const Array4<Real>& phiOld 	= phi_old.array(mfi);
+	       const Array4<Real>& phiPrd 	= phi_prd.array(mfi);
+	      
+	       const amrex::Array4<amrex::Real>& u = umac[0].array(mfi);
+	       const amrex::Array4<amrex::Real>& v = umac[1].array(mfi);
+	       const amrex::Array4<amrex::Real>& u_old = umac_old[0].array(mfi);
+	       const amrex::Array4<amrex::Real>& v_old = umac_old[1].array(mfi);
+#if (AMREX_SPACEDIM == 3)
+	       const amrex::Array4<amrex::Real>& w = umac[2].array(mfi);
+	       const amrex::Array4<amrex::Real>& w_old = umac_old[2].array(mfi);
+#endif
+	       int n=0;
+	       amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+	       {
+	       	
+	           phiNew(i,j,k,n) = phiOld(i,j,k,n)
+	           		+ dt * Mobility * 2 * GradEnCoef *
+	               ( (phiPrd(i+1,j,k) - 2.*phiPrd(i,j,k) + phiPrd(i-1,j,k)) / (dx[0]*dx[0])
+	                +(phiPrd(i,j+1,k) - 2.*phiPrd(i,j,k) + phiPrd(i,j-1,k)) / (dx[1]*dx[1])
+	               )
+	               - dt * ((u(i+1,j,k)+u_old(i+1,j,k))/2*(phiPrd(i+1,j,k)+phiPrd(i,j,k))/2-(u(i,j,k)+u_old(i,j,k))/2*(phiPrd(i,j,k)+phiPrd(i-1,j,k))/2)/dx[0]
+	               - dt * ((v(i,j+1,k)+v_old(i,j+1,k))/2*(phiPrd(i,j+1,k)+phiPrd(i,j,k))/2-(v(i,j,k)+v_old(i,j,k))/2*(phiPrd(i,j,k)+phiPrd(i,j-1,k))/2)/dx[1]
+#if (AMREX_SPACEDIM == 3)
+	               + dt * Mobility * 2 * GradEnCoef * (phiPrd(i,j,k+1) - 2.*phiPrd(i,j,k) + phiPrd(i,j,k-1)) / (dx[2]*dx[2])
+	               - dt * ((w(i,j,k+1)+w_old(i,j,k+1))/2*(phiPrd(i,j,k+1)+phiPrd(i,j,k))/2-(w(i,j,k)+w_old(i,j,k))/2*(phiPrd(i,j,k)+phiPrd(i,j,k-1))/2)/dx[2]
+               	+ variance_coef_mass*sqrt(2.0*k_B*T_init[0]*Mobility*dt/(dx[0]*dx[1]*dx[2]))*amrex::RandomNormal(0,1)
+#elif (AMREX_SPACEDIM == 2)
+               	+ variance_coef_mass*sqrt(2.0*k_B*T_init[0]*Mobility*dt/(dx[0]*dx[1]*cell_depth))*amrex::RandomNormal(0,1)
+#endif
+						- dt * Mobility * EnScale*(2*phiPrd(i,j,k)*(phiPrd(i,j,k)-1)*(phiPrd(i,j,k)-1)+2*phiPrd(i,j,k)*phiPrd(i,j,k)*(phiPrd(i,j,k)-1) - PotWellDepr)
+	               ;
+	            
+			      phiNew(i,j,k,n+1) = 1-phiNew(i,j,k,n);
+	       });
+	   }
+	
+    	ComputeRhotot(phi_new,phitot_new);
+      // compute reversible stress tensor ---added term
+      ComputeDivFHReversibleStress(div_reversible_stress,phitot_new,phi_new,geom);
+	   
+		MultiFab::Copy(phi_old,phi_new,0,0,nspecies,0);
+		
+		amrex::Print() << "Mobility = " << Mobility <<"\n";
+		
+		amrex::Print() << "Advanced phi (phi1_avg = " << phi_new.sum(0)/(n_cells[0]*n_cells[1]
+#if (AMREX_SPACEDIM == 3)
+																				*n_cells[2]
+#endif
+																		 					) <<")\n";
+}
 
     if (use_charged_fluid) {
         // compute total charge at t^{n+1}
@@ -617,14 +754,18 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         }
     }
     
+if (use_ice_nucleation ==0){
     ComputeEta(rho_new, rhotot_new, eta);
+} else {
+    ComputeEta(phi_new, phitot_new, eta);
+}
     if (AMREX_SPACEDIM == 2) {
         AverageCCToNode(eta,eta_ed[0],0,1,SPEC_BC_COMP,geom);
     }
     else {
         AverageCCToEdge(eta,eta_ed,0,1,SPEC_BC_COMP,geom);
     }
-
+    
     //////////////////////////////////////////////
     // Step 6: solve for v^{n+1} and pi^{n+1/2} using GMRES
     //////////////////////////////////////////////
@@ -676,7 +817,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         // call mk_grav_force_bousq(mla,gmres_rhs_v,.true.,rho_fc,the_bc_tower)
     }
 
-    if (use_multiphase || use_flory_huggins) {
+    if (use_multiphase || use_flory_huggins || use_ice_nucleation) {
         // add divergence of reversible stress to gmres_rhs_v
         for (int d=0; d<AMREX_SPACEDIM; ++d) {
             MultiFab::Saxpy(gmres_rhs_v[d],1.,div_reversible_stress[d],0,0,1,0);
@@ -783,6 +924,7 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         eta_ed[d].mult(2.,0,1,0);
     }
 
+if (use_ice_nucleation ==0) {
     // if writing a plotfile for electro-explicit option, compute Epot using new densities
     // use existing machinery to compute all mass fluxes - yes this is overkill
     // but better than writing a separate routine for now
@@ -800,5 +942,6 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
                                charge_old,grad_Epot_old,Epot,permittivity,0);
         }
     }
+ }
     
 }
diff --git a/exec/multispec/Checkpoint.cpp b/exec/multispec/Checkpoint.cpp
index 06d20503..4d46609c 100644
--- a/exec/multispec/Checkpoint.cpp
+++ b/exec/multispec/Checkpoint.cpp
@@ -23,6 +23,8 @@ void WriteCheckPoint(int step,
                      const amrex::Real dt,
                      const MultiFab& rho,
                      const MultiFab& rhotot,
+                     const MultiFab& phi,
+                     const MultiFab& phitot,
                      const MultiFab& pi,
                      std::array< MultiFab, AMREX_SPACEDIM >& umac,
                      const MultiFab& Epot,
@@ -145,6 +147,10 @@ void WriteCheckPoint(int step,
                  amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "rho"));
     VisMF::Write(rhotot,
                  amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "rhotot"));
+    VisMF::Write(phi,
+                 amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "phi"));
+    VisMF::Write(phitot,
+                 amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "phitot"));
     VisMF::Write(pi,
                  amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "pi"));
 
@@ -159,6 +165,8 @@ void ReadCheckPoint(int& step,
                     amrex::Real& dt,
                     MultiFab& rho,
                     MultiFab& rhotot,
+                    MultiFab& phi,
+                    MultiFab& phitot,
                     MultiFab& pi,
                     std::array< MultiFab, AMREX_SPACEDIM >& umac,
                     MultiFab& Epot,
@@ -237,6 +245,8 @@ void ReadCheckPoint(int& step,
         
         rho   .define(ba, dmap, nspecies, ng_s);
         rhotot.define(ba, dmap,        1, ng_s);
+        phi   .define(ba, dmap, nspecies, ng_s);
+        phitot.define(ba, dmap,        1, ng_s);
         pi    .define(ba, dmap,        1, 1);
 
         if (use_charged_fluid) {
@@ -327,6 +337,10 @@ void ReadCheckPoint(int& step,
                 amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "rho"));
     VisMF::Read(rhotot,
                 amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "rhotot"));
+    VisMF::Read(phi,
+                amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "phi"));
+    VisMF::Read(phitot,
+                amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "phitot"));
     VisMF::Read(pi,
                 amrex::MultiFabFileFullPrefix(0, checkpointname, "Level_", "pi"));
 
diff --git a/exec/multispec/GNUmakefile b/exec/multispec/GNUmakefile
index befc7567..f94107b3 100644
--- a/exec/multispec/GNUmakefile
+++ b/exec/multispec/GNUmakefile
@@ -13,7 +13,7 @@ MAX_SPEC      = 8
 # MAX_ELEM needs to be MAX_SPEC*(MAX_SPEC-1)/2
 MAX_ELEM      = 28
 
-TINY_PROFILE = TRUE
+TINY_PROFILE = FALSE
 USE_PARTICLES = FALSE
 
 include $(AMREX_HOME)/Tools/GNUMake/Make.defs
diff --git a/exec/multispec/Init.cpp b/exec/multispec/Init.cpp
index b91870e6..da4671d4 100644
--- a/exec/multispec/Init.cpp
+++ b/exec/multispec/Init.cpp
@@ -7,6 +7,7 @@ using namespace multispec;
 
 void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac,
                  MultiFab& rho_in,
+                 MultiFab& phi_in,
                  const Geometry& geom)
 {
     
@@ -35,7 +36,7 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         Box bx = mfi.tilebox();
 
         const Array4<Real>& c = conc.array(mfi);
-
+				
         if (prob_type == 1) {
 
             /*
@@ -78,7 +79,7 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac,
                 
             });
 
-	} else if (prob_type == 5) {
+			} else if (prob_type == 5) {
 
             /*
               bubble with radius = 1/4 of domain in x
@@ -493,6 +494,68 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac,
                }
             });
 
+        } else if (prob_type == 21) {
+
+            /*
+              cylinder with radius = 1/4 of domain in x
+              c=c_init_1(:) inside, c=c_init_2(:) outside
+              can be discontinous or smooth depending on smoothing_width
+            */
+            //Real rad = L[0] / 8.;
+            int nsub = 100;
+            Real factor = nsub;
+            Real dxsub = dx[0]/factor;
+            Real dysub = dx[1]/factor;
+            Real x,y,z;
+            amrex::Print() << "smoothing width " << smoothing_width << " film_thickness " << film_thickness << std::endl;
+
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+               for (int n=0; n<nspecies; ++n) {
+                   c(i,j,k,n) = 0.;
+               }
+            Real x,y,z;
+
+            for(int i1=0; i1<nsub; ++i1) {
+            for(int j1=0; j1<nsub; ++j1) {
+
+                AMREX_D_TERM(x = prob_lo[0] + i*dx[0] + (i1+0.5)*dxsub ;,
+                             y = prob_lo[1] + j*dx[1] + (j1+0.5)*dysub ;,
+                             z = prob_lo[2] + (k+0.5)*dx[2] - center[2];);
+
+//                Real height = film_thickness + 0.5*film_thickness*amrex::Math::cospi(2.*x/L[0]);
+                Real height = film_thickness + x;
+                //height = film_thickness ;
+
+//              amrex::Print{} << "x,thickness,height " << x << " " << film_thickness << " " << height << std::endl;
+
+                if (smoothing_width == 0.) {
+
+                    // discontinuous interface
+                    if (y < height) {
+                        for (int n=0; n<nspecies; ++n) {
+                            c(i,j,k,n) += c_init_1[n];
+                        }
+                    } else {
+                        for (int n=0; n<nspecies; ++n) {
+                            c(i,j,k,n) += c_init_2[n];
+                        }
+                    }
+
+                } else {
+                    // smooth interface
+                    for (int n=0; n<nspecies; ++n) {
+                        c(i,j,k,n) += c_init_1[n] + (c_init_2[n]-c_init_1[n]) *
+                            0.5*(1. + std::tanh((y-height)/(smoothing_width*dx[0])));
+                    }
+                }
+             }
+             }
+               for (int n=0; n<nspecies; ++n) {
+                   c(i,j,k,n) = c(i,j,k,n)/(factor*factor);
+               }
+            });
+
 
 
         } else if (prob_type == 3) {
@@ -618,12 +681,12 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac,
         } else if (prob_type == 4) {
 
             /*
-              bubble with radius = 1/4 of domain in x
+              slab with thickness = 2*radius_cyl
               c=c_init_1(:) inside, c=c_init_2(:) outside
               can be discontinous or smooth depending on smoothing_width
             */
-            Real l1 = L[1] / 4.;
-            Real l2 = 3.*l1;
+            Real l1 = L[1]*0.5-radius_cyl;
+            Real l2 = L[1]*0.5+radius_cyl;
 
             amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
@@ -654,7 +717,7 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac,
                     for (int n=0; n<nspecies; ++n) {
                         c(i,j,k,n) = c_init_1[n] +
                             (c_init_2[n] - c_init_1[n])*
-                            (1./(1.+std::exp(-smoothing_width*(y-l1))) - 1./(1.+std::exp(-smoothing_width*(y-l2))));
+                            (1./(1.+std::exp(-(y-l1)/(smoothing_width*dx[0]))) - 1./(1.+std::exp(-(y-l2)/(smoothing_width*dx[0]))));
                     }
                 }
             });
@@ -734,6 +797,7 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac,
 
         const Array4<Real>& c = conc.array(mfi);
         const Array4<Real>& rho = rho_in.array(mfi);
+        const Array4<Real>& p = phi_in.array(mfi);
 
         amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
@@ -758,10 +822,11 @@ void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac,
                 }
                 rho_total = 1./sumtot;
             }
-
-            // calculate rho_i
+        			
+            // calculate rho_i and make phi=conc
             for (int n=0; n<nspecies; ++n) {
                 rho(i,j,k,n) = rho_total*c(i,j,k,n);
+					 p(i,j,k,n) = c(i,j,k,n);
             }
                     
         });
diff --git a/exec/multispec/WritePlotFile.cpp b/exec/multispec/WritePlotFile.cpp
index 4d9d09bd..abdadc23 100644
--- a/exec/multispec/WritePlotFile.cpp
+++ b/exec/multispec/WritePlotFile.cpp
@@ -16,7 +16,9 @@ void WritePlotFile(int step,
 		   const MultiFab& rho,
 		   const MultiFab& pres,
                    const MultiFab& charge,
-                   const MultiFab& Epot)
+                   const MultiFab& Epot,
+		   const MultiFab& phitot,
+		   const MultiFab& phi)
 {
     
     BL_PROFILE_VAR("WritePlotFile()",WritePlotFile);
@@ -28,11 +30,13 @@ void WritePlotFile(int step,
 
     // rhotot        1
     // rho           nspecies
+    // phitot        1
+    // phi           nspecies
     // c             nspecies
     // averaged vel  AMREX_SPACEDIM
     // shifted  vel  AMREX_SPACEDIM
     // pres          1
-    int nPlot = 2*AMREX_SPACEDIM + 2*nspecies + 2;
+    int nPlot = 2*AMREX_SPACEDIM + 3*nspecies + 3;
 
     if (use_charged_fluid) {
         // charge
@@ -54,7 +58,15 @@ void WritePlotFile(int step,
         x += (49+i);
         varNames[cnt++] = x;
     }
+    
+    varNames[cnt++] = "phi";
 
+    for (int i=0; i<nspecies; ++i) {
+        std::string x = "phi";
+        x += (49+i);
+        varNames[cnt++] = x;
+    }
+    
     for (int i=0; i<nspecies; ++i) {
         std::string x = "c";
         x += (49+i);
@@ -93,6 +105,16 @@ void WritePlotFile(int step,
         cnt++;
     }
     
+    // copy phitot into plotfile
+    MultiFab::Copy(plotfile, phitot, 0, cnt, 1, 0);
+    cnt++;
+    
+    // copy phi's into plotfile
+    for (int i=0; i<nspecies; ++i) {
+        MultiFab::Copy(plotfile, phi, i, cnt, 1, 0);
+        cnt++;
+    }
+    
     // copy densities and convert to concentrations
     for (int i=0; i<nspecies; ++i) {
         MultiFab::Copy(plotfile, rho, i, cnt, 1, 0);
diff --git a/exec/multispec/_Examples/README b/exec/multispec/_Examples/README
new file mode 100644
index 00000000..cea5cd89
--- /dev/null
+++ b/exec/multispec/_Examples/README
@@ -0,0 +1,13 @@
+This directory contains the inputs files which can be used to study propagation of ice front in a supercooled liquid.
+
+General instructions:
+1. Remember to change the values of thermophysical properties according to the table given within.
+2. You can turn stochasticity on/off by setting variance_coef_* to 1/0.
+
+In ./
+inputs_slab_3d:         inputs to measure the propagation speed of ice/water interface at various temperatures in a pillar geometry.
+inputs_nucleus_3d:      inputs for simulationg growth of an ice nucleus at various temperatures in a cubic domain.
+
+In ./_Validation/
+inputs_laplace:         to measure laplace pressure inside an ice nucleus in a deterministic setting.
+inputs_spectrum:        to measure the spectral energy in a long slab and compare with theory.
diff --git a/exec/multispec/_Examples/_Validation/inputs_laplace b/exec/multispec/_Examples/_Validation/inputs_laplace
new file mode 100644
index 00000000..85cde3fb
--- /dev/null
+++ b/exec/multispec/_Examples/_Validation/inputs_laplace
@@ -0,0 +1,152 @@
+  #----------------------
+  # Problem specification
+  #----------------------
+  prob_lo               = 0.      0.      # physical lo coordinate
+  prob_hi               = 48e-7   48e-7   # physical hi coordinate
+
+  cell_depth            = 1.0e-7
+
+  # refer to Init.cpp
+  prob_type             = 8
+  chk_base_name         = chk
+  plot_base_name        = plt
+  
+  # number of cells in domain. IMP: 0.5 nm or less for deterministic simulations
+  n_cells               = 192 192
+  # max number of cells in a box
+  max_grid_size         = 96 96
+
+  # Time-step control
+  fixed_dt              = 1.e-13
+
+  # Controls for number of steps between actions
+  max_step              = 1000000
+  plot_int              = 10000
+  chk_int               = 100000
+  restart               = -1
+
+  use_ice_nucleation    = 1
+
+  #----------------------
+  # Thermodynamic and transport properties:
+  #------------------------------------------------------------------------------------------------------------------------------------
+  # Temperature (K)	    | 230           235         240         245         250         255         260         265         270
+  #------------------------------------------------------------------------------------------------------------------------------------
+  # c0_water            | 0.1099 	    0.0957 	    0.082  	    0.0687 	    0.0558 	    0.0431 	    0.0308 	    0.0186 	    0.0066 
+  # c0_ice              | 1.0637 	    1.0591 	    1.0539 	    1.048       1.0415 	    1.0341 	    1.0259 	    1.0167 	    1.0064 
+  # mu_water (Poise)    | 0.9124  	    0.2588      0.1279      0.0782      0.0536      0.0394      0.0304      0.0243      0.0200 
+  # fh_chi              | 2.6194        2.5637      2.5103      2.4591      2.4099      2.3626      2.3172      2.2735      2.2314
+  # fh_kappa            | 2.563e-15	    2.508e-15   2.456e-15   2.406e-15   2.358e-15	2.311e-15	2.267e-15	2.224e-15	2.183e-15
+  # rho (g/cm^3)        | 0.9453	    0.9676	    0.9789	    0.9862	    0.9912	    0.9946	    0.997	    0.9986	    0.9995
+  #------------------------------------------------------------------------------------------------------------------------------------
+
+  k_B                   = 1.3806488e-16                     # Boltzmann's constant [units: cm2*g*s-2*K-1]
+  T_init                = 255.0  	                        # [units: K] IMP: Change visc_coef, rho, fh_chi, fh_kappa, c_init_1, c_init_2 accordingly.
+
+  # Viscous friction L phi operator
+  # if abs(visc_type) = 1, L = div beta grad
+  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
+  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
+  # positive = assume constant coefficients
+  # negative = assume spatially-varying coefficients
+  visc_coef             = 0.0394                            # [units: g*cm-1*s-1] dynamic (shear) viscosity of water at T_init
+  visc_type             = 1
+
+  nspecies              = 2
+  #molmass              = 2.98e-23 2.98e-23                 # molecular masses for nspecies (mass per molecule, *not* molar mass)
+  use_flory_huggins     = 1
+
+  fh_monomers           = 1. 1.
+  monomer_mass          = 2.98e-23
+  rho0                  = 0.9946
+  rhobar                = 0.9946 0.9946                     # pure component densities for all species. Same for ice and water for now.
+  fh_chi                = 0. 2.3626     2.3626      0.      # = 2.79e09/(rhobar*k_B*T_init/monomer_mass) ~ 602.478/T_init. Off-diagonal terms should be identical.
+  fh_kappa              = 0. 2.311e-15  2.311e-15   0.      # = 2.73e-6/(rhobar*k_B*T_init/monomer_mass) ~ 5.895e-13/T_init. Off-diagonal terms should be identical.
+
+  is_ideal_mixture      = 0
+
+  # initial values for c
+  c_init_1              = 1.0341 0.0431                     # ice   ~ -1.605e-05*T_init^2+0.006603*T_init+0.3938
+  c_init_2              = 0.0431 1.0341                     # water ~ -0.002575*T_init+0.7005
+  smoothing_width       = 0.0
+  radius_cyl            = 6.e-7
+
+  # These are lower-triangules of symmetric matrices represented as vectors
+  # Number of elements is (nspecies*(nspecies-1)/2)
+  # The values are red row by row starting from top going down
+  # (this allows easy addition/deletion of new species/rows)
+  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
+  Dbar                  = 5.e-5 # Maxwell-Stefan diffusion constant  
+
+  use_charged_fluid     = 0
+  dielectric_const      = 9.2e-20
+  charge_per_mass       = 0. 0.
+
+  mg_verbose            = 0                  # multigrid verbosity
+  algorithm_type        = 6
+  
+  #----------------------
+  # Stochastic parameters
+  #----------------------
+  seed                  = 1
+  variance_coef_mom     = 0.
+  variance_coef_mass    = 0.
+  initial_variance_mom  = 0.
+
+  struct_fact_int       = -1
+  n_steps_skip          = 0
+
+  #----------------------
+  # Boundary conditions
+  # ----------------------
+  # BC specifications:
+  # -1 = periodic
+  #  1 = slip
+  #  2 = no-slip
+  bc_vel_lo             = -1 -1 -1
+  bc_vel_hi             = -1 -1 -1
+
+  # -1 = periodic
+  #  1 = wall
+  #  2 = reservoir (Dirichlet values must be suppled by other means)
+  bc_mass_lo            = -1 -1 -1
+  bc_mass_hi            = -1 -1 -1
+
+  # -1 = periodic
+  #  1 = Dirichlet
+  #  2 = Neumann
+  bc_es_lo              = -1 -1 -1
+  bc_es_hi              = -1 -1 -1
+
+  potential_lo          = 0. 0. 0.
+  potential_hi          = 0. 0. 0.
+
+  #----------------------
+  # Staggered multigrid solver parameters
+  #----------------------
+  stag_mg_verbosity     = 0          # verbosity
+  stag_mg_max_vcycles   = 1         # max number of v-cycles
+  stag_mg_minwidth      = 2            # length of box at coarsest multigrid level
+  stag_mg_bottom_solver = 0       # bottom solver type
+  # 0 = smooths only, controlled by mg_nsmooths_bottom
+  # 4 = Fancy bottom solve that coarsens additionally
+  #     and then applies stag_mg_nsmooths_bottom smooths
+  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
+  stag_mg_nsmooths_up   = 2    # number of smooths at each level on the way up
+  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
+  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
+  stag_mg_omega         = 1.e0            # weighted-jacobi omega coefficient
+  stag_mg_smoother      = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
+  stag_mg_rel_tol       = 1.e-9         # relative tolerance stopping criteria
+
+
+  #----------------------
+  # GMRES solver parameters
+  #----------------------
+  gmres_rel_tol         = 1.e-12                # relative tolerance stopping criteria
+  gmres_abs_tol         = 0                     # absolute tolerance stopping criteria
+  gmres_verbose         = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
+  gmres_max_outer       = 20                  # max number of outer iterations
+  gmres_max_inner       = 5                   # max number of inner iterations, or restart number
+  gmres_max_iter        = 100                  # max number of gmres iterations
+  gmres_min_iter        = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_Examples/_Validation/inputs_spectrum b/exec/multispec/_Examples/_Validation/inputs_spectrum
new file mode 100644
index 00000000..01f41ede
--- /dev/null
+++ b/exec/multispec/_Examples/_Validation/inputs_spectrum
@@ -0,0 +1,152 @@
+  #----------------------
+  # Problem specification
+  #----------------------
+  prob_lo               = 0.        0.      # physical lo coordinate
+  prob_hi               = 256e-7    64e-7   # physical hi coordinate
+
+  cell_depth            = 1.0e-7
+
+  # refer to Init.cpp
+  prob_type             = 16
+  chk_base_name         = chk
+  plot_base_name        = plt
+  
+  # number of cells in domain. IMP: 0.5 nm for deterministic simulations
+  n_cells               = 256 64
+  # max number of cells in a box
+  max_grid_size         = 32 64
+
+  # Time-step control
+  fixed_dt              = 1.e-13
+
+  # Controls for number of steps between actions
+  max_step              = 1000000
+  plot_int              = 1000
+  chk_int               = 100000
+  restart               = -1
+
+  use_ice_nucleation    = 1
+
+  #----------------------
+  # Thermodynamic and transport properties:
+  #------------------------------------------------------------------------------------------------------------------------------------
+  # Temperature (K)	    | 230           235         240         245         250         255         260         265         270
+  #------------------------------------------------------------------------------------------------------------------------------------
+  # c0_water            | 0.1099 	    0.0957 	    0.082  	    0.0687 	    0.0558 	    0.0431 	    0.0308 	    0.0186 	    0.0066 
+  # c0_ice              | 1.0637 	    1.0591 	    1.0539 	    1.048       1.0415 	    1.0341 	    1.0259 	    1.0167 	    1.0064 
+  # mu_water (Poise)    | 0.9124  	    0.2588      0.1279      0.0782      0.0536      0.0394      0.0304      0.0243      0.0200 
+  # fh_chi              | 2.6194        2.5637      2.5103      2.4591      2.4099      2.3626      2.3172      2.2735      2.2314
+  # fh_kappa            | 2.563e-15	    2.508e-15   2.456e-15   2.406e-15   2.358e-15	2.311e-15	2.267e-15	2.224e-15	2.183e-15
+  # rho (g/cm^3)        | 0.9453	    0.9676	    0.9789	    0.9862	    0.9912	    0.9946	    0.997	    0.9986	    0.9995
+  #------------------------------------------------------------------------------------------------------------------------------------
+
+  k_B                   = 1.3806488e-16                     # Boltzmann's constant [units: cm2*g*s-2*K-1]
+  T_init                = 255.0  	                        # [units: K] IMP: Change visc_coef, rho, fh_chi, fh_kappa, c_init_1, c_init_2 accordingly.
+
+  # Viscous friction L phi operator
+  # if abs(visc_type) = 1, L = div beta grad
+  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
+  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
+  # positive = assume constant coefficients
+  # negative = assume spatially-varying coefficients
+  visc_coef             = 0.0394                            # [units: g*cm-1*s-1] dynamic (shear) viscosity of water at T_init
+  visc_type             = 1
+
+  nspecies              = 2
+  #molmass              = 2.98e-23 2.98e-23                 # molecular masses for nspecies (mass per molecule, *not* molar mass)
+  use_flory_huggins     = 1
+
+  fh_monomers           = 1. 1.
+  monomer_mass          = 2.98e-23
+  rho0                  = 0.9946
+  rhobar                = 0.9946 0.9946                     # pure component densities for all species. Same for ice and water for now.
+  fh_chi                = 0. 2.3626     2.3626      0.      # = 2.79e09/(rhobar*k_B*T_init/monomer_mass) ~ 602.478/T_init. Off-diagonal terms should be identical.
+  fh_kappa              = 0. 2.311e-15  2.311e-15   0.      # = 2.73e-6/(rhobar*k_B*T_init/monomer_mass) ~ 5.895e-13/T_init. Off-diagonal terms should be identical.
+
+  is_ideal_mixture      = 0
+
+  # initial values for c
+  c_init_1              = 1.0341 0.0431                     # ice   ~ -1.605e-05*T_init^2+0.006603*T_init+0.3938
+  c_init_2              = 0.0431 1.0341                     # water ~ -0.002575*T_init+0.7005
+  smoothing_width       = 0.75
+  film_thickness        = 32.e-7
+
+  # These are lower-triangules of symmetric matrices represented as vectors
+  # Number of elements is (nspecies*(nspecies-1)/2)
+  # The values are red row by row starting from top going down
+  # (this allows easy addition/deletion of new species/rows)
+  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
+  Dbar                  = 5.e-5 # Maxwell-Stefan diffusion constant  
+
+  use_charged_fluid     = 0
+  dielectric_const      = 9.2e-20
+  charge_per_mass       = 0. 0.
+
+  mg_verbose            = 0                  # multigrid verbosity
+  algorithm_type        = 6
+  
+  #----------------------
+  # Stochastic parameters
+  #----------------------
+  seed                  = 1
+  variance_coef_mom     = 1.
+  variance_coef_mass    = 1.
+  initial_variance_mom  = 0.
+
+  struct_fact_int       = -1
+  n_steps_skip          = 0
+
+  #----------------------
+  # Boundary conditions
+  # ----------------------
+  # BC specifications:
+  # -1 = periodic
+  #  1 = slip
+  #  2 = no-slip
+  bc_vel_lo             = -1 -1 -1
+  bc_vel_hi             = -1 -1 -1
+
+  # -1 = periodic
+  #  1 = wall
+  #  2 = reservoir (Dirichlet values must be suppled by other means)
+  bc_mass_lo            = -1 -1 -1
+  bc_mass_hi            = -1 -1 -1
+
+  # -1 = periodic
+  #  1 = Dirichlet
+  #  2 = Neumann
+  bc_es_lo              = -1 -1 -1
+  bc_es_hi              = -1 -1 -1
+
+  potential_lo          = 0. 0. 0.
+  potential_hi          = 0. 0. 0.
+
+  #----------------------
+  # Staggered multigrid solver parameters
+  #----------------------
+  stag_mg_verbosity     = 0          # verbosity
+  stag_mg_max_vcycles   = 1         # max number of v-cycles
+  stag_mg_minwidth      = 2            # length of box at coarsest multigrid level
+  stag_mg_bottom_solver = 0       # bottom solver type
+  # 0 = smooths only, controlled by mg_nsmooths_bottom
+  # 4 = Fancy bottom solve that coarsens additionally
+  #     and then applies stag_mg_nsmooths_bottom smooths
+  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
+  stag_mg_nsmooths_up   = 2    # number of smooths at each level on the way up
+  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
+  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
+  stag_mg_omega         = 1.e0            # weighted-jacobi omega coefficient
+  stag_mg_smoother      = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
+  stag_mg_rel_tol       = 1.e-9         # relative tolerance stopping criteria
+
+
+  #----------------------
+  # GMRES solver parameters
+  #----------------------
+  gmres_rel_tol         = 1.e-12                # relative tolerance stopping criteria
+  gmres_abs_tol         = 0                     # absolute tolerance stopping criteria
+  gmres_verbose         = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
+  gmres_max_outer       = 20                  # max number of outer iterations
+  gmres_max_inner       = 5                   # max number of inner iterations, or restart number
+  gmres_max_iter        = 100                  # max number of gmres iterations
+  gmres_min_iter        = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_Examples/inputs_nucleus_3d b/exec/multispec/_Examples/inputs_nucleus_3d
new file mode 100644
index 00000000..ef828f47
--- /dev/null
+++ b/exec/multispec/_Examples/inputs_nucleus_3d
@@ -0,0 +1,129 @@
+  #----------------------
+  # Problem specification
+  #----------------------
+  prob_lo               = 0.      0.      0.                # physical lo coordinate
+  prob_hi               = 32e-7   32e-7   32e-7             # physical hi coordinate [units: cm]
+  n_cells               = 32 32 32                          # number of cells in domain. IMP: 0.5 nm for deterministic simulations
+  max_grid_size         = 8 16 16                          # max number of cells in a box
+  fixed_dt              = 1.e-13                            # Time-step control
+
+  prob_type             = 1                                 # refer to Init.cpp
+  use_ice_nucleation    = 1                                 # usual multispec if set to 0
+
+  # Controls for number of steps between actions
+  max_step              = 1000000
+  plot_int              = 1000
+  chk_int               = 100000
+  restart               = -1
+
+  #----------------------
+  # Thermodynamic and transport properties:
+  #------------------------------------------------------------------------------------------------------------------------------------
+  # Temperature (K)	    | 230           235         240         245         250         255         260         265         270
+  #------------------------------------------------------------------------------------------------------------------------------------
+  # c0_water            | 0.1099 	    0.0957 	    0.082  	    0.0687 	    0.0558 	    0.0431 	    0.0308 	    0.0186 	    0.0066 
+  # c0_ice              | 1.0637 	    1.0591 	    1.0539 	    1.048       1.0415 	    1.0341 	    1.0259 	    1.0167 	    1.0064 
+  # mu_water (Poise)    | 0.9124  	    0.2588      0.1279      0.0782      0.0536      0.0394      0.0304      0.0243      0.0200 
+  # fh_chi              | 2.6194        2.5637      2.5103      2.4591      2.4099      2.3626      2.3172      2.2735      2.2314
+  # fh_kappa            | 2.563e-15	    2.508e-15   2.456e-15   2.406e-15   2.358e-15	2.311e-15	2.267e-15	2.224e-15	2.183e-15
+  # rho (g/cm^3)        | 0.9453	    0.9676	    0.9789	    0.9862	    0.9912	    0.9946	    0.997	    0.9986	    0.9995
+  #------------------------------------------------------------------------------------------------------------------------------------
+
+  k_B                   = 1.3806488e-16                     # Boltzmann's constant [units: erg/K]
+  T_init                = 255.0  	                        # [units: K] IMP: Change visc_coef, rho, fh_chi, fh_kappa, c_init_1, c_init_2 accordingly.
+
+  # Viscous friction L phi operator
+  # if abs(visc_type) = 1, L = div beta grad
+  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
+  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
+  # positive = assume constant coefficients
+  # negative = assume spatially-varying coefficients
+  visc_coef             = 0.0394                           # [units: g*cm-1*s-1] dynamic (shear) viscosity of water at T_init
+  visc_type             = 1
+
+  nspecies              = 2
+  #molmass              = 2.98e-23 2.98e-23                 # molecular masses for nspecies (mass per molecule, *not* molar mass)
+
+  fh_monomers           = 1. 1.
+  monomer_mass          = 2.98e-23                          # [units: g]
+  rho0                  = 0.9946                            # [units: g/cm^3]
+  rhobar                = 0.9946 0.9946                     # pure component densities for all species. Same for ice and water for now.
+  fh_chi                = 0. 2.3626     2.3626      0.      # = 2.79e08/(rhobar*k_B*T_init/monomer_mass)    ~ 602.478/T_init.     Off-diagonal terms should be identical.
+  fh_kappa              = 0. 2.311e-15  2.311e-15   0.      # = 2.73e-11/(rhobar*k_B*T_init/monomer_mass)   ~ 5.895e-13/T_init.   Off-diagonal terms should be identical.
+
+  is_ideal_mixture      = 0
+
+  # initial values for c
+  c_init_1              = 1.0341 0.0431                     # ice   ~ -1.605e-05*T_init^2+0.006603*T_init+0.3938
+  c_init_2              = 0.0431 1.0341                     # water ~ -0.002575*T_init+0.7005
+  smoothing_width       = 0.75
+  radius_cyl            = 4.e-7                             # [units: cm]
+
+  mg_verbose            = 0                                 # multigrid verbosity
+  algorithm_type        = 6
+  
+  #----------------------
+  # Stochastic parameters
+  #----------------------
+  seed                  = 1
+  variance_coef_mom     = 1.
+  variance_coef_mass    = 1.
+  initial_variance_mom  = 0.
+
+  struct_fact_int       = -1
+  n_steps_skip          = 0
+
+  #----------------------
+  # Boundary conditions
+  # ----------------------
+  # BC specifications:
+  # -1 = periodic
+  #  1 = slip
+  #  2 = no-slip
+  bc_vel_lo             = -1 -1 -1
+  bc_vel_hi             = -1 -1 -1
+
+  # -1 = periodic
+  #  1 = wall
+  #  2 = reservoir (Dirichlet values must be suppled by other means)
+  bc_mass_lo            = -1 -1 -1
+  bc_mass_hi            = -1 -1 -1
+
+  # -1 = periodic
+  #  1 = Dirichlet
+  #  2 = Neumann
+  bc_es_lo              = -1 -1 -1
+  bc_es_hi              = -1 -1 -1
+
+  potential_lo          = 0. 0. 0.
+  potential_hi          = 0. 0. 0.
+
+  #----------------------
+  # Staggered multigrid solver parameters
+  #----------------------
+  stag_mg_verbosity     = 0          # verbosity
+  stag_mg_max_vcycles   = 1         # max number of v-cycles
+  stag_mg_minwidth      = 2            # length of box at coarsest multigrid level
+  stag_mg_bottom_solver = 0       # bottom solver type
+  # 0 = smooths only, controlled by mg_nsmooths_bottom
+  # 4 = Fancy bottom solve that coarsens additionally
+  #     and then applies stag_mg_nsmooths_bottom smooths
+  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
+  stag_mg_nsmooths_up   = 2    # number of smooths at each level on the way up
+  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
+  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
+  stag_mg_omega         = 1.e0            # weighted-jacobi omega coefficient
+  stag_mg_smoother      = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
+  stag_mg_rel_tol       = 1.e-9         # relative tolerance stopping criteria
+
+
+  #----------------------
+  # GMRES solver parameters
+  #----------------------
+  gmres_rel_tol         = 1.e-12                # relative tolerance stopping criteria
+  gmres_abs_tol         = 0                     # absolute tolerance stopping criteria
+  gmres_verbose         = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
+  gmres_max_outer       = 20                  # max number of outer iterations
+  gmres_max_inner       = 5                   # max number of inner iterations, or restart number
+  gmres_max_iter        = 100                  # max number of gmres iterations
+  gmres_min_iter        = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_Examples/inputs_slab_3d b/exec/multispec/_Examples/inputs_slab_3d
new file mode 100644
index 00000000..8a31d998
--- /dev/null
+++ b/exec/multispec/_Examples/inputs_slab_3d
@@ -0,0 +1,129 @@
+  #----------------------
+  # Problem specification
+  #----------------------
+  prob_lo               = 0.      0.      0.                # physical lo coordinate
+  prob_hi               = 16e-7   32e-7   16e-7             # physical hi coordinate [units: cm]
+  n_cells               = 16 32 16                          # number of cells in domain. IMP: 0.5 nm for deterministic simulations
+  max_grid_size         = 16 16 16                          # max number of cells in a box
+  fixed_dt              = 1.e-13                            # Time-step control
+
+  prob_type             = 4                                 # refer to Init.cpp. 4: slab with thickness = 0.2*length_in_y
+  use_ice_nucleation    = 1                                 # usual multispec if set to 0
+
+  # Controls for number of steps between actions
+  max_step              = 1000000
+  plot_int              = 1000
+  chk_int               = 100000
+  restart               = -1
+
+  #----------------------
+  # Thermodynamic and transport properties:
+  #------------------------------------------------------------------------------------------------------------------------------------
+  # Temperature (K)	    | 230           235         240         245         250         255         260         265         270
+  #------------------------------------------------------------------------------------------------------------------------------------
+  # c0_water            | 0.1099 	    0.0957 	    0.082  	    0.0687 	    0.0558 	    0.0431 	    0.0308 	    0.0186 	    0.0066 
+  # c0_ice              | 1.0637 	    1.0591 	    1.0539 	    1.048       1.0415 	    1.0341 	    1.0259 	    1.0167 	    1.0064 
+  # mu_water (Poise)    | 0.9124  	    0.2588      0.1279      0.0782      0.0536      0.0394      0.0304      0.0243      0.0200 
+  # fh_chi              | 2.6194        2.5637      2.5103      2.4591      2.4099      2.3626      2.3172      2.2735      2.2314
+  # fh_kappa            | 2.563e-15	    2.508e-15   2.456e-15   2.406e-15   2.358e-15	2.311e-15	2.267e-15	2.224e-15	2.183e-15
+  # rho (g/cm^3)        | 0.9453	    0.9676	    0.9789	    0.9862	    0.9912	    0.9946	    0.997	    0.9986	    0.9995
+  #------------------------------------------------------------------------------------------------------------------------------------
+
+  k_B                   = 1.3806488e-16                     # Boltzmann's constant [units: erg/K]
+  T_init                = 245.0  	                        # [units: K] IMP: Change visc_coef, rho, fh_chi, fh_kappa, c_init_1, c_init_2 accordingly.
+
+  # Viscous friction L phi operator
+  # if abs(visc_type) = 1, L = div beta grad
+  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
+  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
+  # positive = assume constant coefficients
+  # negative = assume spatially-varying coefficients
+  visc_coef             = 0.0782                           # [units: g*cm-1*s-1] dynamic (shear) viscosity of water at T_init
+  visc_type             = 1
+
+  nspecies              = 2
+  #molmass              = 2.98e-23 2.98e-23                 # molecular masses for nspecies (mass per molecule, *not* molar mass)
+
+  fh_monomers           = 1. 1.
+  monomer_mass          = 2.98e-23                          # [units: g]
+  rho0                  = 0.9862                            # [units: g/cm^3]
+  rhobar                = 0.9862 0.9862                     # pure component densities for all species. Same for ice and water for now.
+  fh_chi                = 0. 2.4591     2.4591      0.      # = 2.79e08/(rhobar*k_B*T_init/monomer_mass)    ~ 602.478/T_init.     Off-diagonal terms should be identical.
+  fh_kappa              = 0. 2.406e-15  2.406e-15   0.      # = 2.73e-11/(rhobar*k_B*T_init/monomer_mass)   ~ 5.895e-13/T_init.   Off-diagonal terms should be identical.
+
+  is_ideal_mixture      = 0
+
+  # initial values for c
+  c_init_1              = 0.0687 1.048                     # waterPhi ~ -0.002575*T_init+0.7005
+  c_init_2              = 1.0048 0.0687                     # icePhi   ~ -1.605e-05*T_init^2+0.006603*T_init+0.3938
+  smoothing_width       = 0.25
+  radius_cyl            = 3.e-7                             # [units: cm]
+
+  mg_verbose            = 0                                 # multigrid verbosity
+  algorithm_type        = 6
+  
+  #----------------------
+  # Stochastic parameters
+  #----------------------
+  seed                  = 1
+  variance_coef_mom     = 1.
+  variance_coef_mass    = 1.
+  initial_variance_mom  = 0.
+
+  struct_fact_int       = -1
+  n_steps_skip          = 0
+
+  #----------------------
+  # Boundary conditions
+  # ----------------------
+  # BC specifications:
+  # -1 = periodic
+  #  1 = slip
+  #  2 = no-slip
+  bc_vel_lo             = -1 -1 -1
+  bc_vel_hi             = -1 -1 -1
+
+  # -1 = periodic
+  #  1 = wall
+  #  2 = reservoir (Dirichlet values must be suppled by other means)
+  bc_mass_lo            = -1 -1 -1
+  bc_mass_hi            = -1 -1 -1
+
+  # -1 = periodic
+  #  1 = Dirichlet
+  #  2 = Neumann
+  bc_es_lo              = -1 -1 -1
+  bc_es_hi              = -1 -1 -1
+
+  potential_lo          = 0. 0. 0.
+  potential_hi          = 0. 0. 0.
+
+  #----------------------
+  # Staggered multigrid solver parameters
+  #----------------------
+  stag_mg_verbosity     = 0          # verbosity
+  stag_mg_max_vcycles   = 1         # max number of v-cycles
+  stag_mg_minwidth      = 2            # length of box at coarsest multigrid level
+  stag_mg_bottom_solver = 0       # bottom solver type
+  # 0 = smooths only, controlled by mg_nsmooths_bottom
+  # 4 = Fancy bottom solve that coarsens additionally
+  #     and then applies stag_mg_nsmooths_bottom smooths
+  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
+  stag_mg_nsmooths_up   = 2    # number of smooths at each level on the way up
+  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
+  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
+  stag_mg_omega         = 1.e0            # weighted-jacobi omega coefficient
+  stag_mg_smoother      = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
+  stag_mg_rel_tol       = 1.e-9         # relative tolerance stopping criteria
+
+
+  #----------------------
+  # GMRES solver parameters
+  #----------------------
+  gmres_rel_tol         = 1.e-12                # relative tolerance stopping criteria
+  gmres_abs_tol         = 0                     # absolute tolerance stopping criteria
+  gmres_verbose         = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
+  gmres_max_outer       = 20                  # max number of outer iterations
+  gmres_max_inner       = 5                   # max number of inner iterations, or restart number
+  gmres_max_iter        = 100                  # max number of gmres iterations
+  gmres_min_iter        = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/README b/exec/multispec/_PNAS_inputs/README
deleted file mode 100644
index f7da5079..00000000
--- a/exec/multispec/_PNAS_inputs/README
+++ /dev/null
@@ -1,15 +0,0 @@
-This directory contains the inputs files used for the results in "Fluctuating Hydrodynamics and
-the Rayleigh-Plateau Instability
-
-_thickness directory:  inputs to measure interface thickness at different chi
-
-_laplace directory:    inputs to measure interfacial tension from Laplace pressure
-
-inputs_spectrum:       inputs to measure capillary wave spectrum
-
-inputs_RayleighPlateau:  inputs for default Rayleigh Plateau simulations
-                         run with different seeds to create an ensemble
-                         for smaller runs reduce length and numoer of mesh points in z direction
-                             keeping mesh size the same
-
-inputs_high_schmidt:     inputs for default configure but with higher viscosity and Schmidt number
diff --git a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_3.0 b/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_3.0
deleted file mode 100644
index d9d15746..00000000
--- a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_3.0
+++ /dev/null
@@ -1,141 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 8
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 96 96 
-  # max number of cells in a box
-  max_grid_size = 48 48
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 3. 3. 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-  c_init_1 = 0.070720 0.929280
-  c_init_2 = 0.929280 0.070720
-  smoothing_width = .0
-  radius_cyl = 6.e-7
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_3.25 b/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_3.25
deleted file mode 100644
index bc29fde3..00000000
--- a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_3.25
+++ /dev/null
@@ -1,143 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 8
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 96 96
-  # max number of cells in a box
-  max_grid_size = 48 48
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 3.25 3.25 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-  #c_init_1 = 0.034811 0.965188
-  #c_init_2 = 0.965188 0.034811
-  c_init_1 = 0.051361 0.948639
-  c_init_2 = 0.948639 0.051361
-  smoothing_width = .0
-  radius_cyl = 6.e-7
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_3.5 b/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_3.5
deleted file mode 100644
index fb4a8dc4..00000000
--- a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_3.5
+++ /dev/null
@@ -1,143 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 8
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 96 96
-  # max number of cells in a box
-  max_grid_size = 48 48
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 3.5 3.5 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-  #c_init_1 = 0.034811 0.965188
-  #c_init_2 = 0.965188 0.034811
-  c_init_1 = 0.037874 0.962126
-  c_init_2 = 0.962126 0.037874
-  smoothing_width = .0
-  radius_cyl = 6.e-7
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_3.75 b/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_3.75
deleted file mode 100644
index d8cc793c..00000000
--- a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_3.75
+++ /dev/null
@@ -1,141 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 8
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 96 96
-  # max number of cells in a box
-  max_grid_size = 48 48
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 3.75 3.75 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-  c_init_1 = 0.028246 0.971754
-  c_init_2 = 0.971754 0.028246
-  smoothing_width = .0
-  radius_cyl = 6.e-7
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_4.0 b/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_4.0
deleted file mode 100644
index 88fe8d9a..00000000
--- a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_4.0
+++ /dev/null
@@ -1,141 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 8
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 96 96
-  # max number of cells in a box
-  max_grid_size = 48 48
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 4. 4. 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-  c_init_1 = 0.021248 0.978752
-  c_init_2 = 0.978752 0.021248
-  smoothing_width = .0
-  radius_cyl = 6.e-7
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_4.25 b/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_4.25
deleted file mode 100644
index 33309b35..00000000
--- a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_4.25
+++ /dev/null
@@ -1,141 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 8
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 96 96
-  # max number of cells in a box
-  max_grid_size = 48 48
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 4.25 4.25 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-  c_init_1 = 0.016092 0.983908
-  c_init_2 = 0.983908 0.016092
-  smoothing_width = .0
-  radius_cyl = 6.e-7
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_4.5 b/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_4.5
deleted file mode 100644
index adf036b7..00000000
--- a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_4.5
+++ /dev/null
@@ -1,141 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 8
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 96 96
-  # max number of cells in a box
-  max_grid_size = 48 48
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 4.5 4.5 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-  c_init_1 = 0.012252 0.987748
-  c_init_2 = 0.987748 0.012252
-  smoothing_width = .0
-  radius_cyl = 6.e-7
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_4.75 b/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_4.75
deleted file mode 100644
index 37018cd0..00000000
--- a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_4.75
+++ /dev/null
@@ -1,141 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 8
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 96 96
-  # max number of cells in a box
-  max_grid_size = 48 48
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 4.75 4.75 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-  c_init_1 = 0.009368 0.990632
-  c_init_2 = 0.990632 0.009368
-  smoothing_width = .0
-  radius_cyl = 6.e-7
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_5.0 b/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_5.0
deleted file mode 100644
index a29a5278..00000000
--- a/exec/multispec/_PNAS_inputs/_laplace_pressure/inputs_bubble_chi_5.0
+++ /dev/null
@@ -1,141 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 8
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 96 96
-  # max number of cells in a box
-  max_grid_size = 48 48
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 5. 5. 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-  c_init_1 = 0.007188 0.992812
-  c_init_2 = 0.992812 0.007188
-  smoothing_width = .0
-  radius_cyl = 6.e-7
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_3.0 b/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_3.0
deleted file mode 100644
index 231eb978..00000000
--- a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_3.0
+++ /dev/null
@@ -1,143 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  #prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-  prob_hi = 3.2e-6 6.4e-6 
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 16
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 32 64
-  # max number of cells in a box
-  max_grid_size = 32 32
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 3. 3. 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-  c_init_1 = 0.070720 0.929280
-  c_init_2 = 0.929280 0.070720
-  smoothing_width = .0
-  #radius_cyl = 6.e-7
-  film_thickness = 3.2e-6
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_3.25 b/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_3.25
deleted file mode 100644
index 709d839a..00000000
--- a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_3.25
+++ /dev/null
@@ -1,150 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  #prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-  prob_hi = 3.2e-6 6.4e-6 
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 16
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 32 64
-  # max number of cells in a box
-  max_grid_size = 32 32
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 3.25 3.25 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-
-  c_init_1 = 0.051361 0.948639
-  c_init_2 = 0.948639 0.051361
-
-
-
-
-
-
-  smoothing_width = .0
-  #radius_cyl = 6.e-7
-  film_thickness = 3.2e-6
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_3.5 b/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_3.5
deleted file mode 100644
index 61b2b8b0..00000000
--- a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_3.5
+++ /dev/null
@@ -1,150 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  #prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-  prob_hi = 3.2e-6 6.4e-6 
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 16
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 32 64
-  # max number of cells in a box
-  max_grid_size = 32 32
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 3.5 3.5 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-
-  c_init_1 = 0.037874 0.962126
-  c_init_2 = 0.962126 0.037874
-
-
-
-
-
-
-  smoothing_width = .0
-  #radius_cyl = 6.e-7
-  film_thickness = 3.2e-6
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_3.75 b/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_3.75
deleted file mode 100644
index a26ff437..00000000
--- a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_3.75
+++ /dev/null
@@ -1,151 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  #prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-  prob_hi = 3.2e-6 6.4e-6 
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 16
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 32 64
-  # max number of cells in a box
-  max_grid_size = 32 32
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 3.75 3.75 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-
-  c_init_1 = 0.028246 0.971754
-  c_init_2 = 0.971754 0.028246
-
-
-
-
-
-
-
-  smoothing_width = .0
-  #radius_cyl = 6.e-7
-  film_thickness = 3.2e-6
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_4.0 b/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_4.0
deleted file mode 100644
index 97026f56..00000000
--- a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_4.0
+++ /dev/null
@@ -1,152 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  #prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-  prob_hi = 3.2e-6 6.4e-6 
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 16
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 32 64
-  # max number of cells in a box
-  max_grid_size = 32 32
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 4. 4. 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-
-
-  c_init_1 = 0.021248 0.978752
-  c_init_2 = 0.978752 0.021248
-
-
-
-
-
-
-
-  smoothing_width = .0
-  #radius_cyl = 6.e-7
-  film_thickness = 3.2e-6
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_4.25 b/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_4.25
deleted file mode 100644
index f30a3e3e..00000000
--- a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_4.25
+++ /dev/null
@@ -1,153 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  #prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-  prob_hi = 3.2e-6 6.4e-6 
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 16
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 32 64
-  # max number of cells in a box
-  max_grid_size = 32 32
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 4.25 4.25 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-
-
-  c_init_1 = 0.016092 0.983908
-  c_init_2 = 0.983908 0.016092
-
-
-
-
-
-
-
-
-  smoothing_width = .0
-  #radius_cyl = 6.e-7
-  film_thickness = 3.2e-6
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_4.5 b/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_4.5
deleted file mode 100644
index 5ec45aa6..00000000
--- a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_4.5
+++ /dev/null
@@ -1,153 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  #prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-  prob_hi = 3.2e-6 6.4e-6 
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 16
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 32 64
-  # max number of cells in a box
-  max_grid_size = 32 32
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 4.5 4.5 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-
-
-  c_init_1 = 0.012252 0.987748
-  c_init_2 = 0.987748 0.012252
-
-
-
-
-
-
-
-
-  smoothing_width = .0
-  #radius_cyl = 6.e-7
-  film_thickness = 3.2e-6
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_4.75 b/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_4.75
deleted file mode 100644
index c05c5b28..00000000
--- a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_4.75
+++ /dev/null
@@ -1,153 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  #prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-  prob_hi = 3.2e-6 6.4e-6 
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 16
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 32 64
-  # max number of cells in a box
-  max_grid_size = 32 32
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 4.75 4.75 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-
-
-
-  c_init_1 = 0.009368 0.990632
-  c_init_2 = 0.990632 0.009368
-
-
-
-
-
-
-
-  smoothing_width = .0
-  #radius_cyl = 6.e-7
-  film_thickness = 3.2e-6
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_5.0 b/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_5.0
deleted file mode 100644
index f9af1451..00000000
--- a/exec/multispec/_PNAS_inputs/_thickness/inputs_slab_chi_5.0
+++ /dev/null
@@ -1,149 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  #prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-  prob_hi = 3.2e-6 6.4e-6 
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 16
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 32 64
-  # max number of cells in a box
-  max_grid_size = 32 32
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 1000
-  chk_int = 10000
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 0.
-  variance_coef_mass = 0.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  #fh_chi = 0. 3.571 3.571 0.
-  fh_chi = 0. 5. 5. 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-
-  c_init_1 = 0.007188 0.992812
-  c_init_2 = 0.992812 0.007188
-
-
-
-
-
-  smoothing_width = .0
-  #radius_cyl = 6.e-7
-  film_thickness = 3.2e-6
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/inputs_RayleighPlateau b/exec/multispec/_PNAS_inputs/inputs_RayleighPlateau
deleted file mode 100644
index 74d2c881..00000000
--- a/exec/multispec/_PNAS_inputs/inputs_RayleighPlateau
+++ /dev/null
@@ -1,141 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  #prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-  prob_hi = 4.8e-6 4.8e-6 360.e-7  # physical hi coordinate
-
-  # cell_depth = 0.5e-7
-
-  # refer to Init.cpp
-  prob_type = 8
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 48 48 360
-  # max number of cells in a box
-  max_grid_size = 48 48 24
-
-  # Time-step control
-  fixed_dt = 1.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 100
-  chk_int = 100
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.0 2464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 1.
-  variance_coef_mass = 1.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0 	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = .6e-22
-  use_flory_huggins = 1
-  fh_chi = 0. 3.571 3.571 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-  c_init_1 = 0.0348115 0.9651885
-  c_init_2 = 0.9651885 0.0348115
-  smoothing_width = 1.
-  radius_cyl = 6.e-7
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/inputs_high_schmidt b/exec/multispec/_PNAS_inputs/inputs_high_schmidt
deleted file mode 100644
index 3b4ec590..00000000
--- a/exec/multispec/_PNAS_inputs/inputs_high_schmidt
+++ /dev/null
@@ -1,141 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  #prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-  prob_hi = 4.8e-6 4.8e-6 360.e-7  # physical hi coordinate
-
-  # cell_depth = 0.5e-7
-
-  # refer to Init.cpp
-  prob_type = 8
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 48 48 360
-  # max number of cells in a box
-  max_grid_size = 48 48 24
-
-  # Time-step control
-  fixed_dt = 1.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 100
-  chk_int = 100
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.02464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 1.
-  variance_coef_mass = 1.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0 	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = .6e-22
-  use_flory_huggins = 1
-  fh_chi = 0. 3.571 3.571 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-  c_init_1 = 0.0348115 0.9651885
-  c_init_2 = 0.9651885 0.0348115
-  smoothing_width = 1.
-  radius_cyl = 6.e-7
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/inputs_spectrum b/exec/multispec/_PNAS_inputs/inputs_spectrum
deleted file mode 100644
index 8b327a14..00000000
--- a/exec/multispec/_PNAS_inputs/inputs_spectrum
+++ /dev/null
@@ -1,142 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  #prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-  #prob_hi = 4.8e-6 4.8e-6 4.8e-6  # physical hi coordinate
-  prob_hi =  25.6e-6 6.4e-6   # physical hi coordinate
-
-  cell_depth = 1.0e-7
-
-  # refer to Init.cpp
-  prob_type = 16
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 256 64 
-  # max number of cells in a box
-  max_grid_size = 64 64
-
-  # Time-step control
-  fixed_dt = 1.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 100
-  chk_int = 100
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 1.
-  variance_coef_mass = 1.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0 	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 2 -1
-  bc_vel_hi = -1 2 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 1 -1
-  bc_mass_hi = -1 1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = .6e-22
-  use_flory_huggins = 1
-  fh_chi = 0. 3.571 3.571 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-  c_init_1 = 0.0348115 0.9651885
-  c_init_2 = 0.9651885 0.0348115
-  smoothing_width = .75
-  film_thickness = 6.e-7
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 5.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/_PNAS_inputs/inputs_torus b/exec/multispec/_PNAS_inputs/inputs_torus
deleted file mode 100644
index c0147b99..00000000
--- a/exec/multispec/_PNAS_inputs/inputs_torus
+++ /dev/null
@@ -1,141 +0,0 @@
-
-
-  # Problem specification
-  prob_lo = 0.      0.       0. # physical lo coordinate
-  #prob_hi = 9.6e-7 9.6e-7 9.6e-7  # physical hi coordinate
-  prob_hi = 2.56e-5 2.56e-5 6.4e-6
-
-  # cell_depth = 0.5e-7
-
-  # refer to Init.cpp
-  prob_type = 9
-  chk_base_name = chk
-  plot_base_name= plt
-  
-  # number of cells in domain
-  n_cells = 256 256 64
-  # max number of cells in a box
-  max_grid_size = 64 64 64
-
-  # Time-step control
-  fixed_dt = 4.e-13
-
-  # Controls for number of steps between actions
-  max_step = 1000000
-  plot_int = 100
-  chk_int = 100
-  restart = -1
-
-  # Viscous friction L phi operator
-  # if abs(visc_type) = 1, L = div beta grad
-  # if abs(visc_type) = 2, L = div [ beta (grad + grad^T) ]
-  # if abs(visc_type) = 3, L = div [ beta (grad + grad^T) + I (gamma - (2/3)*beta) div ]
-  # positive = assume constant coefficients
-  # negative = assume spatially-varying coefficients
-  visc_coef = 0.002464   # [units: g*cm-1*s-1] dynamic (shear) viscosity of water
-  visc_type = 1
-  
-  # Stochastic parameters
-  seed = 1
-  variance_coef_mom = 1.
-  variance_coef_mass = 1.
-  initial_variance_mom = 0.
-
-  k_B = 1.3806488e-16   # Boltzmann's constant [units: cm2*g*s-2*K-1]
-  T_init = 84.0  	# [units: K]
-
-  algorithm_type = 6
-
-  struct_fact_int = -1
-  n_steps_skip = 0
-
-  # Boundary conditions
-  # ----------------------
-  # BC specifications:
-  # -1 = periodic
-  #  1 = slip
-  #  2 = no-slip
-  bc_vel_lo = -1 -1 -1
-  bc_vel_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = wall
-  #  2 = reservoir (Dirichlet values must be suppled by other means)
-  bc_mass_lo = -1 -1 -1
-  bc_mass_hi = -1 -1 -1
-
-  # -1 = periodic
-  #  1 = Dirichlet
-  #  2 = Neumann
-  bc_es_lo = -1 -1 -1
-  bc_es_hi = -1 -1 -1
-
-  potential_lo = 0. 0. 0.
-  potential_hi = 0. 0. 0.
-
-  # Thermodynamic and transport properties:
-  #----------------------
-  nspecies = 2
-  #molmass = 2.4e-22 2.4e-22  # molecular masses for nspecies (mass per molecule, *not* molar mass)
-  fh_monomers = 1. 1.
-  monomer_mass = 6.e-23
-  use_flory_huggins = 1
-  fh_chi = 0. 3.571 3.571 0.
-  fh_kappa = 0. 2.7e-14 2.7e-14 0.
-
-  rhobar = 1.4 1.4     # pure component densities for all species
-  rho0 = 1.4
-
-
-
-  is_ideal_mixture = 0
-
-  # initial values for c
-  c_init_1 = 0.034811 0.965188
-  c_init_2 = 0.965188 0.034811
-  smoothing_width = 1.
-  radius_cyl = 6e-7
-
-  # These are lower-triangules of symmetric matrices represented as vectors
-  # Number of elements is (nspecies*(nspecies-1)/2)
-  # The values are red row by row starting from top going down
-  # (this allows easy addition/deletion of new species/rows)
-  # So D_12; D_13, D_23; D_14, D_24, D_34; ...
-  Dbar = 1.e-5 # Maxwell-Stefan diffusion constant  
-
-  use_charged_fluid = 0
-  dielectric_const = 9.2e-20
-  charge_per_mass = 0. 0.
-
-  #use_multiphase = 1
-  #alpha_gex = 4.
-  #n_gex = 1
-  #kc_tension = 4.e-12
-
-  mg_verbose = 0                  # multigrid verbosity
-
-  # Staggered multigrid solver parameters
-  stag_mg_verbosity = 0          # verbosity
-  stag_mg_max_vcycles = 1         # max number of v-cycles
-  stag_mg_minwidth = 2            # length of box at coarsest multigrid level
-  stag_mg_bottom_solver = 0       # bottom solver type
-  # 0 = smooths only, controlled by mg_nsmooths_bottom
-  # 4 = Fancy bottom solve that coarsens additionally
-  #     and then applies stag_mg_nsmooths_bottom smooths
-  stag_mg_nsmooths_down = 2  # number of smooths at each level on the way down
-  stag_mg_nsmooths_up = 2    # number of smooths at each level on the way up
-  stag_mg_nsmooths_bottom = 8     # number of smooths at the bottom
-  stag_mg_max_bottom_nlevels = 10 # for stag_mg_bottom_solver 4, number of additional levels of multigrid
-  stag_mg_omega = 1.e0            # weighted-jacobi omega coefficient
-  stag_mg_smoother = 1            # 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
-  stag_mg_rel_tol = 1.e-9         # relative tolerance stopping criteria
-
-
-  # GMRES solver parameters
-  gmres_rel_tol = 1.e-12                # relative tolerance stopping criteria
-  gmres_abs_tol = 0                     # absolute tolerance stopping criteria
-  gmres_verbose = 1                     # gmres verbosity; if greater than 1, more residuals will be printed out
-  gmres_max_outer = 20                  # max number of outer iterations
-  gmres_max_inner = 5                   # max number of inner iterations, or restart number
-  gmres_max_iter = 100                  # max number of gmres iterations
-  gmres_min_iter = 1                    # min number of gmres iterations
diff --git a/exec/multispec/main_driver.cpp b/exec/multispec/main_driver.cpp
index 60ee3f38..be949736 100644
--- a/exec/multispec/main_driver.cpp
+++ b/exec/multispec/main_driver.cpp
@@ -129,9 +129,13 @@ void main_driver(const char* argv)
     MultiFab Epot;
     std::array< MultiFab, AMREX_SPACEDIM > grad_Epot_old;
     
+    /***/
+    MultiFab phi_old;
+    MultiFab phitot_old;
+    
     if (restart > 0) {
         ReadCheckPoint(init_step,time,dt,
-                       rho_old,rhotot_old,pi,umac,Epot,grad_Epot_old,
+                       rho_old,rhotot_old,phi_old,phitot_old,pi,umac,Epot,grad_Epot_old,
                        ba,dmap);
     }
     else {
@@ -171,6 +175,10 @@ void main_driver(const char* argv)
         if (use_charged_fluid) {
             Epot.define(ba, dmap, 1, 1);
         }
+        
+    	  /***/
+        phi_old.define 		(ba, dmap, nspecies, ng_s);
+        phitot_old.define	(ba, dmap, 1       , ng_s);
     }
 
     // moved this to here so can change dt from value in checkpoint
@@ -191,8 +199,9 @@ void main_driver(const char* argv)
 
     if (restart < 0) {
     
+    	  /***/
         // initialize rho and umac in valid region only
-        InitRhoUmac(umac,rho_old,geom);
+        InitRhoUmac(umac,rho_old,phi_old,geom);
 
         // initialize pi, including ghost cells
         pi.setVal(0.);
@@ -207,6 +216,9 @@ void main_driver(const char* argv)
     // pressure ghost cells
     pi.FillBoundary(geom.periodicity());
     MultiFabPhysBC(pi,geom,0,1,PRES_BC_COMP);
+    	  /***/
+    ComputeRhotot(phi_old,phitot_old);
+    FillRhoRhototGhost(phi_old,phitot_old,geom);
 
     //=======================================================
     // Build multifabs for all the variables
@@ -219,6 +231,11 @@ void main_driver(const char* argv)
     MultiFab eta              (ba, dmap, 1       , 1);
     MultiFab kappa            (ba, dmap, 1       , 1);
     
+    	  /***/
+    MultiFab phi_new       	(ba, dmap, nspecies, ng_s);
+    MultiFab phitot_new       (ba, dmap, 1       , ng_s);
+    //MultiFab local_MobScale   (ba, dmap, 1		 , 0);
+    //local_MobScale.setVal(1.0);	// Change this near to wall BC
     /////////////////////////////////////////
 
     // eta and Temp on nodes (2d) or edges (3d)
@@ -506,7 +523,7 @@ void main_driver(const char* argv)
         
         // write initial plotfile and structure factor
         if (plot_int > 0) {
-            WritePlotFile(0,0.,geom,umac,rhotot_old,rho_old,pi,charge_old,Epot);
+            WritePlotFile(0,0.,geom,umac,rhotot_old,rho_old,pi,charge_old,Epot,phitot_old,phi_old);
             if (n_steps_skip == 0 && struct_fact_int > 0) {
                 structFact.WritePlotFile(0,0.,geom,"plt_SF");
             }
@@ -548,12 +565,14 @@ void main_driver(const char* argv)
         }
         else if (algorithm_type == 6) {
             // boussinesq
+    	  /***/
             AdvanceTimestepBousq(umac,rho_old,rho_new,rhotot_old,rhotot_new,
                                  pi,eta,eta_ed,kappa,Temp,Temp_ed,
                                  diff_mass_fluxdiv,stoch_mass_fluxdiv,stoch_mass_flux,
                                  grad_Epot_old,grad_Epot_new,
                                  charge_old,charge_new,Epot,permittivity,
                                  sMassFlux,sMomFlux,
+                                 phi_old,phi_new,phitot_old,phitot_new,
                                  dt,time,istep,geom);
         }
         else {
@@ -591,7 +610,7 @@ void main_driver(const char* argv)
 
         // write plotfile at specific intervals
         if (plot_int > 0 && istep%plot_int == 0) {
-            WritePlotFile(istep,time,geom,umac,rhotot_new,rho_new,pi,charge_new,Epot);
+            WritePlotFile(istep,time,geom,umac,rhotot_new,rho_new,pi,charge_new,Epot,phitot_new,phi_new);
             if (istep > n_steps_skip && struct_fact_int > 0) {
                 structFact.WritePlotFile(istep,time,geom,"plt_SF");
             }
@@ -599,7 +618,7 @@ void main_driver(const char* argv)
 
         // write checkpoint at specific intervals
         if (chk_int > 0 && istep%chk_int == 0) {
-            WriteCheckPoint(istep,time,dt,rho_new,rhotot_new,pi,umac,Epot,grad_Epot_new);
+            WriteCheckPoint(istep,time,dt,rho_new,rhotot_new,phi_new,phitot_new,pi,umac,Epot,grad_Epot_new);
         }
 
         // set old state to new state
diff --git a/exec/multispec/multispec_test_functions.H b/exec/multispec/multispec_test_functions.H
index 2c0cb2b8..0dac6e50 100644
--- a/exec/multispec/multispec_test_functions.H
+++ b/exec/multispec/multispec_test_functions.H
@@ -60,6 +60,10 @@ void AdvanceTimestepBousq(std::array< MultiFab, AMREX_SPACEDIM >& umac,
                           MultiFab& permittivity,
                           StochMassFlux& sMassFlux,
                           StochMomFlux& sMomFlux,
+                          MultiFab& phi_old,
+                          MultiFab& phi_new,
+                          MultiFab& phitot_old,
+                          MultiFab& phitot_new,
                           const Real& dt,
                           const Real& time,
                           const int& istep,
@@ -73,6 +77,8 @@ void WriteCheckPoint(int step,
                      const amrex::Real dt,
                      const MultiFab& rho,
                      const MultiFab& rhotot,
+                     const MultiFab& phi,
+                     const MultiFab& phitot,
                      const MultiFab& pi,
                      std::array< MultiFab, AMREX_SPACEDIM >& umac,
                      const MultiFab& Epot,
@@ -83,6 +89,8 @@ void ReadCheckPoint(int& step,
                     amrex::Real& dt,
                     MultiFab& rho,
                     MultiFab& rhotot,
+                    MultiFab& phi,
+                    MultiFab& phitot,
                     MultiFab& pi,
                     std::array< MultiFab, AMREX_SPACEDIM >& umac,
                     MultiFab& Epot,
@@ -98,6 +106,7 @@ void ReadFile(const std::string& filename, Vector<char>& charBuf,
 
 void InitRhoUmac(std::array< MultiFab, AMREX_SPACEDIM >& umac,
                  MultiFab& rho_in,
+                 MultiFab& phi_old,
                  const Geometry& geom);
 
 ///////////////////////////
@@ -108,10 +117,12 @@ void WritePlotFile(const int step,
                    const Geometry geom,
                    std::array< MultiFab, AMREX_SPACEDIM >& umac,
                    const MultiFab& rhotot,
-		   const MultiFab& rho,
+		   			 const MultiFab& rho,
                    const MultiFab& pres,
                    const MultiFab& charge,
-                   const MultiFab& Epot);
+                   const MultiFab& Epot,
+                   const MultiFab& phitot,
+		   			 const MultiFab& phi);
 
 ///////////////////////////
 
diff --git a/src_multispec/InitialProjection.cpp b/src_multispec/InitialProjection.cpp
index 4aa1a869..7f7567df 100644
--- a/src_multispec/InitialProjection.cpp
+++ b/src_multispec/InitialProjection.cpp
@@ -95,11 +95,12 @@ void InitialProjection(std::array< MultiFab, AMREX_SPACEDIM >& umac,
     if (variance_coef_mass != 0. && algorithm_type != 6) {
         sMassFlux.fillMassStochastic();
     }
-        
+
+if (use_ice_nucleation ==0) {        
     ComputeMassFluxdiv(rho,rhotot,Temp,diff_mass_fluxdiv,stoch_mass_fluxdiv,
                        diff_mass_flux,stoch_mass_flux,sMassFlux,dt_eff,time,geom,weights,
                        charge_old,grad_Epot_old,Epot,permittivity);
-
+}
     // assumble total fluxes to be used in reservoirs
     //
     //
diff --git a/src_multispec/multispec_functions.cpp b/src_multispec/multispec_functions.cpp
index 436d051b..d7d18d45 100644
--- a/src_multispec/multispec_functions.cpp
+++ b/src_multispec/multispec_functions.cpp
@@ -18,11 +18,20 @@ AMREX_GPU_MANAGED int                                       multispec::is_ideal_
 int                                                         multispec::use_lapack;
 AMREX_GPU_MANAGED int                                       multispec::use_multiphase;
 AMREX_GPU_MANAGED int                                       multispec::use_flory_huggins;
+AMREX_GPU_MANAGED int                                       multispec::use_ice_nucleation;
+AMREX_GPU_MANAGED int                                       multispec::use_log_potential;
+AMREX_GPU_MANAGED amrex::Real                               multispec::EnScale;
+AMREX_GPU_MANAGED amrex::Real                               multispec::GradEnCoef;
+AMREX_GPU_MANAGED amrex::Real                               multispec::PotWellDepr;
+AMREX_GPU_MANAGED amrex::Real                               multispec::Mobility;
 AMREX_GPU_MANAGED amrex::Real                               multispec::kc_tension;
 AMREX_GPU_MANAGED amrex::Real                               multispec::alpha_gex;
 AMREX_GPU_MANAGED amrex::Array2D<Real,0,MAX_SPECIES-1,0,MAX_SPECIES-1> multispec::fh_kappa;
 AMREX_GPU_MANAGED amrex::Array2D<Real,0,MAX_SPECIES-1,0,MAX_SPECIES-1> multispec::fh_chi;
 AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, MAX_SPECIES> multispec::fh_monomers;
+AMREX_GPU_MANAGED amrex::Real                               multispec::fh_tension;
+AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>      multispec::contact_angle_lo;
+AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>      multispec::contact_angle_hi;
 AMREX_GPU_MANAGED amrex::Real                               multispec::monomer_mass;
 AMREX_GPU_MANAGED int                                       multispec::n_gex;
 AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, MAX_SPECIES> multispec::c_init_1;
@@ -80,10 +89,22 @@ void InitializeMultispecNamespace() {
     use_lapack = 0;         // Use LAPACK or iterative method for diffusion matrix (recommend False)
     use_multiphase = 0;     // for RTIL
     use_flory_huggins = 0;   // for flory huggins
+    use_ice_nucleation 	= 0;   																			// for ice nucleation simulations
+    use_log_potential 	= 0;																				// use quartic potential by default
+    EnScale 				= 2.79e09;
+    GradEnCoef 			= 2.73e-6;
+    PotWellDepr 			= -2.687e-05*T_init[0]*T_init[0]+0.009943*T_init[0]-0.7128;		// For polynomial potential
+    Mobility 				= 0.0001808*T_init[0]*T_init[0]-0.08481*T_init[0]+9.968;			// AC Mobility
     kc_tension = 0;         // for RTIL
     alpha_gex = 0;          // for RTIL
     n_gex = 1;              // for RTIL
     chi_iterations = 10;    // number of iterations used in Dbar2chi_iterative
+    fh_tension = 0.;
+
+    for (int i=0; i<AMREX_SPACEDIM; ++i) {
+        contact_angle_lo[i] = 90.;        // spacedim-vector specifying external E field
+        contact_angle_hi[i] = 90.;        // spacedim-vector specifying external E field
+    }
 
     // Initial and boundary conditions 
     //----------------------
@@ -183,8 +204,12 @@ void InitializeMultispecNamespace() {
     pp.query("use_lapack",use_lapack);
     pp.query("use_multiphase",use_multiphase);
     pp.query("use_flory_huggins",use_flory_huggins);
+    pp.query("use_ice_nucleation",use_ice_nucleation);
+    pp.query("use_log_potential",use_log_potential);
     pp.query("monomer_mass",monomer_mass);
     pp.query("kc_tension",kc_tension);
+    pp.query("Mobility",Mobility);
+    pp.query("fh_tension",fh_tension);
     pp.query("alpha_gex",alpha_gex);
     pp.query("n_gex",n_gex);
     pp.query("chi_iterations",chi_iterations);
@@ -203,6 +228,17 @@ void InitializeMultispecNamespace() {
         }
         }
     }
+   // contact angles in degrees.  converty to radians
+    if(pp.queryarr("contact_angle_lo",temp)) {
+        for (int i=0; i<AMREX_SPACEDIM; ++i) {
+            contact_angle_lo[i] = temp[i]*M_PI/180.;
+        }
+    }
+    if(pp.queryarr("contact_angle_hi",temp)) {
+        for (int i=0; i<AMREX_SPACEDIM; ++i) {
+            contact_angle_hi[i] = temp[i]*M_PI/180.;
+        }
+    }
     if(pp.queryarr("fh_monomers",temp)) {
         for (int i=0; i<nspecies; ++i) {
             fh_monomers[i] = temp[i];
diff --git a/src_multispec/multispec_namespace.H b/src_multispec/multispec_namespace.H
index 6d89ebd5..3d01c3fe 100644
--- a/src_multispec/multispec_namespace.H
+++ b/src_multispec/multispec_namespace.H
@@ -24,6 +24,9 @@ namespace multispec {
     extern AMREX_GPU_MANAGED amrex::Array2D<amrex::Real,0,MAX_SPECIES-1,0,MAX_SPECIES-1>  fh_kappa;
     extern AMREX_GPU_MANAGED amrex::Array2D<amrex::Real,0,MAX_SPECIES-1,0,MAX_SPECIES-1>  fh_chi;
     extern AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, MAX_SPECIES> fh_monomers;
+    extern AMREX_GPU_MANAGED amrex::Real fh_tension;
+    extern AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>          contact_angle_lo;
+    extern AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>          contact_angle_hi;
     extern AMREX_GPU_MANAGED amrex::Real monomer_mass;
     extern AMREX_GPU_MANAGED amrex::Real kc_tension;
     extern AMREX_GPU_MANAGED amrex::Real alpha_gex;
@@ -57,5 +60,13 @@ namespace multispec {
     extern amrex::Real                L_pos;
     extern amrex::Real                L_trans;
     extern amrex::Real                L_zero;
+    
+    // ice-nucleation
+    extern AMREX_GPU_MANAGED int      		use_ice_nucleation;
+    extern AMREX_GPU_MANAGED int      		use_log_potential;
+    extern AMREX_GPU_MANAGED amrex::Real 	EnScale;
+    extern AMREX_GPU_MANAGED amrex::Real 	GradEnCoef;
+    extern AMREX_GPU_MANAGED amrex::Real 	PotWellDepr;
+    extern AMREX_GPU_MANAGED amrex::Real 	Mobility;
 }