diff --git a/exec/immersedIons/main_driver.cpp b/exec/immersedIons/main_driver.cpp
index 51dd39c8f..cb099f475 100644
--- a/exec/immersedIons/main_driver.cpp
+++ b/exec/immersedIons/main_driver.cpp
@@ -339,7 +339,7 @@ void main_driver(const char* argv)
     double simParticles = 0;
     double wetRad[nspecies];
     double dxAv = (dx[0] + dx[1] + dx[2])/3.0; //This is probably the wrong way to do this.
-    std::string line;
+    //std::string line;
 
     for(int j=0;j<nspecies;j++) {
        if (pkernel_fluid[j] == 3) {
@@ -838,11 +838,80 @@ void main_driver(const char* argv)
         efieldCC[d].setVal(0.);
     }
 
-    // external field
-    std::array< MultiFab, AMREX_SPACEDIM > external;
+    // external field cell-centered
+    std::array< MultiFab, AMREX_SPACEDIM > externalCC;
+    for (int d=0; d<AMREX_SPACEDIM; ++d) {
+        externalCC[d].define(bp, dmap, 1, ngp);
+    }
+
+    // external field face-centered
+    std::array< MultiFab, AMREX_SPACEDIM > externalFC;
+    for (int d=0; d<AMREX_SPACEDIM; ++d) {
+        externalFC[d].define(convert(bp,nodal_flag_dir[d]), dmap, 1, ngp);
+    }
+
+    // Set beta_es and permittivity_fc equal to permittivity
+    std::array< MultiFab, AMREX_SPACEDIM > beta_es;
+    std::array< MultiFab, AMREX_SPACEDIM > permittivity_fc;
+    for (int d=0; d<AMREX_SPACEDIM; ++d) {
+        beta_es[d].define(convert(ba,nodal_flag_dir[d]), dmap, 1, 0);
+        permittivity_fc[d].define(convert(ba,nodal_flag_dir[d]), dmap, 1, 0);
+    }
+
+    for (int d=0; d<AMREX_SPACEDIM; ++d) {
+        beta_es[d].setVal(permittivity);
+        permittivity_fc[d].setVal(permittivity);
+    }
+
+    // Read in space-varying permittivity info
+    int i_mf, j_mf, k_mf;
+    Real new_permittivity;
+    std::map< IntVect, Real > permittivity_map;
+    int epscount = 0;
+
+    std::string filename = "epsilon_mem.dat";
+    std::ifstream permittivityFile(filename); // permittivity file for non-solvent part
+    std::string line;
+    while (std::getline(permittivityFile, line)) {
+        std::istringstream ss(line);
+        ss >> i_mf;
+        ss >> j_mf;
+        ss >> k_mf;
+        ss >> new_permittivity;
+        IntVect mf_idx = IntVect(AMREX_D_DECL(i_mf,j_mf,k_mf));
+        permittivity_map.insert(make_pair(mf_idx, new_permittivity));
+        epscount++;
+    }
+    Print() << "Loaded " << epscount << " cells of different permittivity." << std::endl;
+    permittivityFile.close();
+
+    // Update non-uniform permittivity
     for (int d=0; d<AMREX_SPACEDIM; ++d) {
-        external[d].define(bp, dmap, 1, ngp);
+        for (MFIter mfi(beta_es[d]); mfi.isValid(); ++mfi) {
+	    auto& d_fab = beta_es[d][mfi];
+#ifdef AMREX_USE_GPU
+            FArrayBox h_fab(d_fab.box(), The_Pinned_Arena());
+            auto const& h_data = h_fab.array();
+#else
+            auto const& h_data = d_fab.array();
+#endif
+            amrex::LoopOnCpu(d_fab.box(), [=] (int i, int j, int k)
+            {
+                auto search = permittivity_map.find(IntVect(i,j,k));
+                if (search != permittivity_map.end()) {
+                    h_data(i,j,k) = search->second;
+                } else {
+                    h_data(i,j,k) = permittivity; // default value from input file
+                }
+            });
+#ifdef AMREX_USE_GPU
+            Gpu::htod_memcpy_async(d_fab.dataPtr(), h_fab.dataPtr(), h_fab.nBytes());
+            Gpu::streamSynchronize();
+#endif
+
+	}
     }
+
     
     ///////////////////////////////////////////
     // structure factor for charge-charge
@@ -1316,7 +1385,8 @@ void main_driver(const char* argv)
         //Most of these functions are sensitive to the order of execution. We can fix this, but for now leave them in this order.
 
         for (int d=0; d<AMREX_SPACEDIM; ++d) {
-            external[d].setVal(eamp[d]*cos(efreq[d]*time + ephase[d]));  // external field
+            externalCC[d].setVal(eamp[d]*cos(efreq[d]*time + ephase[d]));  // external field
+            externalFC[d].setVal(eamp[d]*cos(efreq[d]*time + ephase[d]));  // staggered external field
             source    [d].setVal(body_force_density[d]);      // reset source terms
             sourceTemp[d].setVal(0.0);      // reset source terms
             sourceRFD[d].setVal(0.0);      // reset source terms
@@ -1362,7 +1432,7 @@ void main_driver(const char* argv)
         
         // do Poisson solve using 'charge' for RHS, and put potential in 'potential'.
         // Then calculate gradient and put in 'efieldCC', then add 'external'.
-        esSolve(potential, charge, efieldCC, external, geomP);
+        esSolve(potential, charge, efieldCC, permittivity_fc, beta_es, externalCC, externalFC, geomP);
 
         if (es_tog==2) {
             // compute pairwise Coulomb force (currently hard-coded to work with y-wall).
diff --git a/exec/immersedIons/particleGen.cpp b/exec/immersedIons/particleGen.cpp
index 007ac4cd6..babb62863 100644
--- a/exec/immersedIons/particleGen.cpp
+++ b/exec/immersedIons/particleGen.cpp
@@ -31,9 +31,11 @@ void FhdParticleContainer::InitParticles(species* particleInfo, const Real* dxp)
 
             proc0_enter = false;
 
-            std::ifstream particleFile("particles.dat");
  //           Print() << "SPEC TOTAL: " << particleInfo[0].total << "\n";
             for(int i_spec=0; i_spec < nspecies; i_spec++) {
+	        std::string filename = Concatenate("particles_",i_spec);
+	        filename += ".dat";
+                std::ifstream particleFile(filename);
                 for (int i_part=0; i_part<particleInfo[i_spec].total;i_part++) {
                     ParticleType p;
                     p.id()  = ParticleType::NextID();
@@ -147,9 +149,10 @@ void FhdParticleContainer::InitParticles(species* particleInfo, const Real* dxp)
 
                     pcount++;
                 }
+
+                particleFile.close();
             }
 
-            particleFile.close();
         }
     }
 
diff --git a/src_common/ComputeDivAndGrad.cpp b/src_common/ComputeDivAndGrad.cpp
index e5236549c..6d69730cb 100644
--- a/src_common/ComputeDivAndGrad.cpp
+++ b/src_common/ComputeDivAndGrad.cpp
@@ -5,7 +5,7 @@
 void ComputeDiv(MultiFab& div,
                 const std::array<MultiFab, AMREX_SPACEDIM>& phi_fc,
                 int start_incomp, int start_outcomp, int ncomp,
-                const Geometry& geom, int increment)
+                const Geometry& geom, Real factor, int increment)
 {
     BL_PROFILE_VAR("ComputeDiv()",ComputeDiv);
 
@@ -27,10 +27,15 @@ void ComputeDiv(MultiFab& div,
 
         amrex::ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
         {
-            div_fab(i,j,k,start_outcomp+n) +=
+	    Real temp = 
                 AMREX_D_TERM(  (phix_fab(i+1,j,k,start_incomp+n) - phix_fab(i,j,k,start_incomp+n)) / dx[0],
                              + (phiy_fab(i,j+1,k,start_incomp+n) - phiy_fab(i,j,k,start_incomp+n)) / dx[1],
                              + (phiz_fab(i,j,k+1,start_incomp+n) - phiz_fab(i,j,k,start_incomp+n)) / dx[2]);;
+	    if (factor == 0.) {
+                div_fab(i,j,k,start_outcomp+n) += temp;
+	    } else {
+                div_fab(i,j,k,start_outcomp+n) += factor*temp;
+	    }
         });
     }
 }
@@ -183,7 +188,7 @@ void ComputeGrad(const MultiFab & phi_in, std::array<MultiFab, AMREX_SPACEDIM> &
 // Computes gradient at cell centres from cell centred data - ouputs to a three component mf.
 void ComputeCentredGrad(const MultiFab & phi,
                         std::array<MultiFab, AMREX_SPACEDIM> & gphi,
-                        const Geometry & geom)
+                        const Geometry & geom, Real factor)
 {
     BL_PROFILE_VAR("ComputeCentredGrad()",ComputeCentredGrad);
 
@@ -201,9 +206,9 @@ void ComputeCentredGrad(const MultiFab & phi,
 
         amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
-            AMREX_D_TERM(gphix_fab(i,j,k) = (phi_fab(i+1,j,k) - phi_fab(i-1,j,k) ) / (2.*dx[0]);,
-                         gphiy_fab(i,j,k) = (phi_fab(i,j+1,k) - phi_fab(i,j-1,k) ) / (2.*dx[1]);,
-                         gphiz_fab(i,j,k) = (phi_fab(i,j,k+1) - phi_fab(i,j,k-1) ) / (2.*dx[2]););
+            AMREX_D_TERM(gphix_fab(i,j,k) = factor*(phi_fab(i+1,j,k) - phi_fab(i-1,j,k) ) / (2.*dx[0]);,
+                         gphiy_fab(i,j,k) = factor*(phi_fab(i,j+1,k) - phi_fab(i,j-1,k) ) / (2.*dx[1]);,
+                         gphiz_fab(i,j,k) = factor*(phi_fab(i,j,k+1) - phi_fab(i,j,k-1) ) / (2.*dx[2]););
         });
     }
 }
diff --git a/src_common/common_functions.H b/src_common/common_functions.H
index a3d1df052..ba89882b0 100644
--- a/src_common/common_functions.H
+++ b/src_common/common_functions.H
@@ -65,14 +65,14 @@ void AverageCCToEdge(const MultiFab& cc_in, std::array<MultiFab, NUM_EDGE>& edge
 
 void ComputeDiv(MultiFab & div, const std::array<MultiFab, AMREX_SPACEDIM> & phi_fc, 
                 int start_incomp, int start_outcomp, int ncomp, 
-                const Geometry & geom, int increment=0);
+                const Geometry & geom, Real factor, int increment=0);
 
 void ComputeGrad(const MultiFab & phi_in, std::array<MultiFab, AMREX_SPACEDIM> & gphi, 
                  int start_incomp, int start_outcomp, int ncomp, int bccomp, const Geometry & geom,
                  int increment=0);
 
 void ComputeCentredGrad(const MultiFab& phi, std::array<MultiFab, AMREX_SPACEDIM>& gphi,
-                        const Geometry& geom);
+                        const Geometry& geom, Real factor);
 
 void ComputeCentredGradCompDir(const MultiFab & phi,
                                MultiFab& gphi,
diff --git a/src_electrostatic/electrostatic.H b/src_electrostatic/electrostatic.H
index 9a66e327c..3830562f0 100644
--- a/src_electrostatic/electrostatic.H
+++ b/src_electrostatic/electrostatic.H
@@ -10,7 +10,11 @@
 using namespace amrex;
 
 void esSolve(MultiFab& potential, MultiFab& charge, std::array< MultiFab, AMREX_SPACEDIM >& efieldCC,
-             const std::array< MultiFab, AMREX_SPACEDIM >& external, const Geometry geom);
+	     const std::array< MultiFab, AMREX_SPACEDIM >& permittivity_fc,
+             const std::array< MultiFab, AMREX_SPACEDIM >& beta_es,
+             const std::array< MultiFab, AMREX_SPACEDIM >& externalCC,
+             const std::array< MultiFab, AMREX_SPACEDIM >& externalFC, 
+	     const Geometry geom);
 
 void calculateField(MultiFab& potential, const Geometry geom);
 
diff --git a/src_electrostatic/electrostatic.cpp b/src_electrostatic/electrostatic.cpp
index db6009ae9..ed5c90e92 100644
--- a/src_electrostatic/electrostatic.cpp
+++ b/src_electrostatic/electrostatic.cpp
@@ -1,12 +1,17 @@
 #include "electrostatic.H"
 #include "common_functions.H"
 #include <AMReX_MLMG.H>
+#include <AMReX_MLABecLaplacian.H>
 
 using namespace amrex;
 
 void esSolve(MultiFab& potential, MultiFab& charge,
              std::array< MultiFab, AMREX_SPACEDIM >& efieldCC,
-             const std::array< MultiFab, AMREX_SPACEDIM >& external, const Geometry geom)
+	     const std::array< MultiFab, AMREX_SPACEDIM >& permittivity_fc,
+	     const std::array< MultiFab, AMREX_SPACEDIM >& beta_es,
+             const std::array< MultiFab, AMREX_SPACEDIM >& externalCC, 
+             const std::array< MultiFab, AMREX_SPACEDIM >& externalFC, 
+	     const Geometry geom)
 {
     BL_PROFILE_VAR("esSolve()",esSolve);
 
@@ -45,11 +50,35 @@ void esSolve(MultiFab& potential, MultiFab& charge,
             }
         }
 
+	potential.setVal(0.);
         const BoxArray& ba = charge.boxArray();
         const DistributionMapping& dmap = charge.DistributionMap();
 
+	// Copy charge in the RHS of Poisson Solver
+        MultiFab charge_unscaled(ba,dmap,1,charge.nGrow()); // create new charge (unscaled) MFab
+        MultiFab::Copy(charge_unscaled,charge,0,0,1,charge.nGrow()); // copy contents of charge/eps into charge
+        charge_unscaled.mult(permittivity,charge.nGrow()); // multiply the charge value by permittvity to get correct charge
+        MultiFab rhs(ba,dmap,1,0); // create RHS MFab
+        MultiFab::Copy(rhs,charge_unscaled,0,0,1,0); // Copy charge into rhs
+	
+	std::array< MultiFab, AMREX_SPACEDIM > epsE_ext;
+        for (int d=0; d<AMREX_SPACEDIM; ++d) {
+            epsE_ext[d].define(convert(ba,nodal_flag_dir[d]), dmap, 1, 0);
+            MultiFab::Copy(epsE_ext[d],permittivity_fc[d],0,0,1,permittivity_fc[d].nGrow()); // copy contents of charge/eps into charge
+        }
+	
+        // compute epsilon*external (where epsilon is stored in permittivity_fc)
+        for (int i=0; i<AMREX_SPACEDIM; ++i) {
+            MultiFab::Multiply(epsE_ext[i],externalFC[i],0,0,1,0);
+        }
+
+        // compute div (epsilon*E_ext) and subtract it from the solver rhs
+        // only needed for spatially varying epsilon or external field
+        ComputeDiv(rhs,epsE_ext,0,0,1,geom,-1.,0);
+
         //create solver opject
-        MLPoisson linop({geom}, {ba}, {dmap});
+        //MLPoisson linop({geom}, {ba}, {dmap});
+	MLABecLaplacian linop({geom}, {ba}, {dmap});
  
         //set BCs
         linop.setDomainBC({AMREX_D_DECL(lo_linop_bc[0],
@@ -71,6 +100,12 @@ void esSolve(MultiFab& potential, MultiFab& charge,
         // be correct or the Poisson solver won't converge
         linop.setEnforceSingularSolvable(false);
 
+	// set alpha=0, beta=1 (will overwrite beta with epsilon next)
+        linop.setScalars(0.0, 1.0);
+
+        // set beta=epsilon
+        linop.setBCoeffs(0, amrex::GetArrOfConstPtrs(beta_es));
+
         //Multi Level Multi Grid
         MLMG mlmg(linop);
 
@@ -80,7 +115,7 @@ void esSolve(MultiFab& potential, MultiFab& charge,
         mlmg.setBottomVerbose(poisson_bottom_verbose);
         
         //Do solve
-        mlmg.solve({&potential}, {&charge}, poisson_rel_tol, 0.0);
+        mlmg.solve({&potential}, {&rhs}, poisson_rel_tol, 0.0);
             
         potential.FillBoundary(geom.periodicity());
         // set ghost cell values so electric field is calculated properly
@@ -88,19 +123,17 @@ void esSolve(MultiFab& potential, MultiFab& charge,
         MultiFabPotentialBC(potential, geom); 
 
         //Find e field, gradient from cell centers to faces
-        ComputeCentredGrad(potential, efieldCC, geom);
+        ComputeCentredGrad(potential, efieldCC, geom, -1.);
         
 
     }
 
     //Add external field on top, then fill boundaries, then setup BCs for peskin interpolation
     for (int d=0; d<AMREX_SPACEDIM; ++d) {
-        efieldCC[d].mult(-1.0,efieldCC[d].nGrow());
-        MultiFab::Add(efieldCC[d], external[d], 0, 0, 1, efieldCC[d].nGrow());
+        //efieldCC[d].mult(-1.0,efieldCC[d].nGrow());
+        MultiFab::Add(efieldCC[d], externalCC[d], 0, 0, 1, efieldCC[d].nGrow());
         efieldCC[d].FillBoundary(geom.periodicity());
         MultiFabElectricBC(efieldCC[d], geom);
     }
 
 }
-
-