diff --git a/libpympb/pympb.cpp b/libpympb/pympb.cpp
index 4f83726bf..40ac8dafa 100644
--- a/libpympb/pympb.cpp
+++ b/libpympb/pympb.cpp
@@ -185,14 +185,11 @@ static void deflation_constraint(evectmatrix X, void *data) {
 /******* mode_solver *******/
 
 mode_solver::mode_solver(int num_bands,
-                         int parity,
                          double resolution[3],
                          lattice lat,
                          double tolerance,
                          int mesh_size,
                          meep_geom::material_data *_default_material,
-                         geometric_object_list geom,
-                         bool reset_fields,
                          bool deterministic,
                          double target_freq,
                          int dims,
@@ -208,7 +205,6 @@ mode_solver::mode_solver(int num_bands,
                          int eigensolver_nwork,
                          int eigensolver_block_size):
   num_bands(num_bands),
-  parity(parity),
   target_freq(target_freq),
   tolerance(tolerance),
   mesh_size(mesh_size),
@@ -218,6 +214,7 @@ mode_solver::mode_solver(int num_bands,
   force_mu(force_mu),
   use_simple_preconditioner(use_simple_preconditioner),
   grid_size(grid_size),
+  nwork_alloc(0),
   eigensolver_nwork(eigensolver_nwork),
   eigensolver_block_size(eigensolver_block_size),
   last_parity(-2),
@@ -235,11 +232,17 @@ mode_solver::mode_solver(int num_bands,
   curfield_type('-'),
   eps(true) {
 
-  this->lat = lat;
-
   geometry_lattice = lat;
   dimensions = dims;
   ensure_periodicity = periodicity;
+  geometry_tree = NULL;
+  H.data = NULL;
+  Hblock.data = NULL;
+  muinvH.data = NULL;
+
+  for (int i = 0; i < MAX_NWORK; i++) {
+    W[i].data = NULL;
+  }
 
   for (int i = 0; i < 3; ++i) {
     this->resolution[i] = resolution[i];
@@ -249,21 +252,14 @@ mode_solver::mode_solver(int num_bands,
     }
   }
 
-  default_material = _default_material;
-
 #ifndef WITH_HERMITIAN_EPSILON
   meep_geom::medium_struct *m;
-  if (meep_geom::is_medium(default_material, &m)) {
+  if (meep_geom::is_medium(_default_material, &m)) {
     meep_geom::check_offdiag(m);
   }
+#else
+  (void)_default_material;
 #endif
-
-  geometry = geom;
-
-  // `init` is called in the constructor to avoid the need to copy the
-  // geometric objects and default material. They can then be safely freed by
-  // typemaps once the mode_solver constructor call returns to python.
-  init(parity, reset_fields);
 }
 
 mode_solver::~mode_solver() {
@@ -716,9 +712,12 @@ bool mode_solver::using_mu() {
   return mdata && mdata->mu_inv != NULL;
 }
 
-void mode_solver::init(int p, bool reset_fields) {
+void mode_solver::init(int p, bool reset_fields, geometric_object_list geometry,
+                       meep_geom::material_data *_default_material) {
   int have_old_fields = 0;
 
+  default_material = _default_material;
+
   n[0] = grid_size.x;
   n[1] = grid_size.y;
   n[2] = grid_size.z;
@@ -759,6 +758,36 @@ void mode_solver::init(int p, bool reset_fields) {
     block_size = num_bands;
   }
 
+  if (mdata) {
+    if (n[0] == mdata->nx && n[1] == mdata->ny && n[2] == mdata->nz &&
+        block_size == Hblock.alloc_p && num_bands == H.p &&
+        eigensolver_nwork + (mdata->mu_inv != NULL) == nwork_alloc) {
+
+      have_old_fields = 1;
+    }
+    else {
+      destroy_evectmatrix(H);
+      for (int i = 0; i < nwork_alloc; ++i) {
+        destroy_evectmatrix(W[i]);
+      }
+
+      if (Hblock.data != H.data) {
+        destroy_evectmatrix(Hblock);
+      }
+      if (muinvH.data != H.data) {
+        destroy_evectmatrix(muinvH);
+      }
+    }
+    destroy_maxwell_target_data(mtdata);
+    mtdata = NULL;
+    destroy_maxwell_data(mdata);
+    mdata = NULL;
+    curfield_reset();
+  }
+  else {
+    srand(time(NULL));
+  }
+
   if (deterministic) {
     // seed should be the same for each run, although
     // it should be different for each process.
@@ -774,7 +803,7 @@ void mode_solver::init(int p, bool reset_fields) {
     mtdata = create_maxwell_target_data(mdata, target_freq);
   }
 
-  init_epsilon();
+  init_epsilon(&geometry);
 
   if (check_maxwell_dielectric(mdata, 0)) {
     meep::abort("invalid dielectric function for MPB");
@@ -816,7 +845,7 @@ void mode_solver::init(int p, bool reset_fields) {
   evectmatrix_flops = eigensolver_flops;
 }
 
-void mode_solver::init_epsilon() {
+void mode_solver::init_epsilon(geometric_object_list *geometry) {
   mpb_real no_size_x = geometry_lattice.size.x == 0 ? 1 : geometry_lattice.size.x;
   mpb_real no_size_y = geometry_lattice.size.y == 0 ? 1 : geometry_lattice.size.y;
   mpb_real no_size_z = geometry_lattice.size.z == 0 ? 1 : geometry_lattice.size.z;
@@ -844,20 +873,20 @@ void mode_solver::init_epsilon() {
   matrix3x3_to_arr(R, Rm);
   matrix3x3_to_arr(G, Gm);
 
-  geom_fix_objects0(geometry);
+  geom_fix_objects0(*geometry);
 
   meep::master_printf("Geometric objects:\n");
   if (meep::am_master()) {
-    for (int i = 0; i < geometry.num_items; ++i) {
+    for (int i = 0; i < geometry->num_items; ++i) {
 
 #ifndef WITH_HERMITIAN_EPSILON
       meep_geom::medium_struct *mm;
-      if (meep_geom::is_medium(geometry.items[i].material, &mm)) {
+      if (meep_geom::is_medium(geometry->items[i].material, &mm)) {
         meep_geom::check_offdiag(mm);
       }
 #endif
 
-      display_geometric_object_info(5, geometry.items[i]);
+      display_geometric_object_info(5, geometry->items[i]);
 
       // meep_geom::medium_struct *mm;
       // if (meep_geom::is_medium(geometry.items[i].material, &mm)) {
@@ -884,7 +913,7 @@ void mode_solver::init_epsilon() {
     b0.high.x += tmp_size.x / mdata->nx;
     b0.high.y += tmp_size.y / mdata->ny;
     b0.high.z += tmp_size.z / mdata->nz;
-    geometry_tree = create_geom_box_tree0(geometry, b0);
+    geometry_tree = create_geom_box_tree0(*geometry, b0);
   }
 
   if (verbose && meep::am_master()) {
@@ -896,14 +925,14 @@ void mode_solver::init_epsilon() {
   int tree_nobjects;
   geom_box_tree_stats(geometry_tree, &tree_depth, &tree_nobjects);
   meep::master_printf("Geometric object tree has depth %d and %d object nodes"
-                      " (vs. %d actual objects)\n", tree_depth, tree_nobjects, geometry.num_items);
+                      " (vs. %d actual objects)\n", tree_depth, tree_nobjects, geometry->num_items);
 
   // restricted_tree = geometry_tree;
 
-  reset_epsilon();
+  reset_epsilon(geometry);
 }
 
-void mode_solver::reset_epsilon() {
+void mode_solver::reset_epsilon(geometric_object_list *geometry) {
   int mesh[3] = {
     mesh_size,
     (dimensions > 1) ? mesh_size : 1,
@@ -922,7 +951,7 @@ void mode_solver::reset_epsilon() {
   set_maxwell_dielectric(mdata, mesh, R, G, dielectric_function, mean_epsilon_func,
                          static_cast<void *>(this));
 
-  if (has_mu()) {
+  if (has_mu(geometry)) {
     meep::master_printf("Initializing mu function...\n");
     eps = false;
     set_maxwell_mu(mdata, mesh, R, G, dielectric_function, mean_epsilon_func,
@@ -931,14 +960,14 @@ void mode_solver::reset_epsilon() {
   }
 }
 
-bool mode_solver::has_mu() {
+bool mode_solver::has_mu(geometric_object_list *geometry) {
   // TODO: mu_file_func
   if (material_has_mu(default_material) || force_mu) {
     return true;
   }
 
-  for (int i = 0; i < geometry.num_items; ++i) {
-    if (material_has_mu(geometry.items[i].material)) {
+  for (int i = 0; i < geometry->num_items; ++i) {
+    if (material_has_mu(geometry->items[i].material)) {
       return true;
     }
   }
@@ -1001,6 +1030,11 @@ void mode_solver::set_parity(integer p) {
   set_kpoint_index(0);  /* reset index */
 }
 
+void mode_solver::set_num_bands(int nb) {
+  num_bands = nb;
+  freqs.resize(nb);
+}
+
 int mode_solver::get_kpoint_index() {
   return kpoint_index;
 }
@@ -1586,6 +1620,32 @@ std::vector<int> mode_solver::get_dims() {
   return dims;
 }
 
+void mode_solver::set_grid_size(vector3 gs) {
+  grid_size.x = gs.x;
+  grid_size.y = gs.y;
+  grid_size.z = gs.z;
+}
+
+int mode_solver::get_libctl_dimensions() {
+  return dimensions;
+}
+
+void mode_solver::set_libctl_dimensions(int val) {
+  dimensions = val;
+}
+
+bool mode_solver::get_libctl_ensure_periodicity() {
+  return ensure_periodicity;
+}
+
+void mode_solver::set_libctl_ensure_periodicity(bool val) {
+  ensure_periodicity = val;
+}
+
+void mode_solver::set_libctl_geometry_lattice(lattice val) {
+  geometry_lattice = val;
+}
+
 void mode_solver::get_curfield(double *data, int size) {
   mpb_real *p = (mpb_real *)curfield;
 
diff --git a/libpympb/pympb.hpp b/libpympb/pympb.hpp
index 86473ad4e..c32c8877d 100644
--- a/libpympb/pympb.hpp
+++ b/libpympb/pympb.hpp
@@ -29,7 +29,6 @@ struct mode_solver {
   static const int NUM_FFT_BANDS = 20;
 
   int num_bands;
-  int parity;
   double resolution[3];
   double target_freq;
   lattice lat;
@@ -57,10 +56,7 @@ struct mode_solver {
   int iterations;
   double eigensolver_flops;
 
-  geometric_object_list geometry;
-
   geom_box_tree geometry_tree;
-  // geom_box_tree restricted_tree;
 
   mpb_real vol;
 
@@ -81,25 +77,39 @@ struct mode_solver {
   evectmatrix muinvH;
   evectmatrix W[MAX_NWORK];
 
-  meep_geom::material_data *default_md;
-
   std::vector<mpb_real> freqs;
 
   bool verbose;
   bool deterministic;
 
-  mode_solver(int num_bands, int parity, double resolution[3], lattice lat, double tolerance,
-              int mesh_size, meep_geom::material_data *_default_material, geometric_object_list geom,
-              bool reset_fields, bool deterministic, double target_freq, int dims, bool verbose,
-              bool periodicity, double flops, bool negative_epsilon_ok, std::string epsilon_input_file,
-              std::string mu_input_file, bool force_mu, bool use_simple_preconditioner, vector3 grid_size,
-              int eigensolver_nwork, int eigensolver_block_size);
+  mode_solver(int num_bands,
+              double resolution[3],
+              lattice lat,
+              double tolerance,
+              int mesh_size,
+              meep_geom::material_data *_default_material,
+              bool deterministic,
+              double target_freq,
+              int dims,
+              bool verbose,
+              bool periodicity,
+              double flops,
+              bool negative_epsilon_ok,
+              std::string epsilon_input_file,
+              std::string mu_input_file,
+              bool force_mu,
+              bool use_simple_preconditioner,
+              vector3 grid_size,
+              int eigensolver_nwork,
+              int eigensolver_block_size);
+
   ~mode_solver();
 
-  void init(int p, bool reset_fields);
+  void init(int p, bool reset_fields, geometric_object_list geometry, meep_geom::material_data *_default_material);
   void solve_kpoint(vector3 kpoint);
   bool using_mu();
   void set_parity(int p);
+  void set_num_bands(int nb);
   int get_kpoint_index();
   void set_kpoint_index(int i);
   void get_epsilon();
@@ -112,9 +122,9 @@ struct mode_solver {
                     mpb_real d1, mpb_real d2, mpb_real d3, mpb_real tol, const mpb_real r[3]);
 
   void randomize_fields();
-  void init_epsilon();
-  void reset_epsilon();
-  bool has_mu();
+  void init_epsilon(geometric_object_list *geometry);
+  void reset_epsilon(geometric_object_list *geometry);
+  bool has_mu(geometric_object_list *geometry);
   bool material_has_mu(void *mt);
   void curfield_reset();
 
@@ -146,6 +156,12 @@ struct mode_solver {
   void set_curfield_type(char t);
   std::string get_parity_string();
   std::vector<int> get_dims();
+  void set_grid_size(vector3 gs);
+  int get_libctl_dimensions();
+  void set_libctl_dimensions(int val);
+  bool get_libctl_ensure_periodicity();
+  void set_libctl_ensure_periodicity(bool val);
+  void set_libctl_geometry_lattice(lattice val);
   std::vector<mpb_real> get_output_k();
 
   mpb_real get_val(int ix, int iy, int iz, int nx, int ny, int nz, int last_dim_size,
diff --git a/python/geom.py b/python/geom.py
index cea1619a5..d2c34d04b 100644
--- a/python/geom.py
+++ b/python/geom.py
@@ -752,8 +752,7 @@ def newton(x, a, b, dx):
         a_prime = x if f < 0 else a
         b_prime = x if f > 0 else b
 
-        cond = f_memo(lazy(a_prime))[0] * f_memo(lazy(b_prime))[0] > 0
-        if dx != x_max - x_min and dx * (f / df) < 0 and cond:
+        if dx != x_max - x_min and dx * (f / df) < 0 and f_memo(lazy(a_prime))[0] * f_memo(lazy(b_prime))[0] > 0:
             raise ValueError("failed to bracket the root in find_root_deriv")
 
         if isinstance(a, numbers.Number) and isinstance(b, numbers.Number):
diff --git a/python/solver.py b/python/solver.py
index 97cce21e1..0fcf760c9 100644
--- a/python/solver.py
+++ b/python/solver.py
@@ -98,46 +98,63 @@ def __init__(self,
                  eigensolver_nwork=3,
                  eigensolver_block_size=-11):
 
+        self.mode_solver = None
         self.resolution = resolution
-        self.is_negative_epsilon_ok = is_negative_epsilon_ok
-        self.eigensolver_flops = eigensolver_flops
-        self.eigensolver_nwork = eigensolver_nwork
-        self.eigensolver_block_size = eigensolver_block_size
         self.eigensolver_flags = eigensolver_flags
-        self.use_simple_preconditioner = use_simple_preconditioner
-        self.force_mu = force_mu
-        self.mu_input_file = mu_input_file
-        self.epsilon_input_file = epsilon_input_file
-        self.mesh_size = mesh_size
-        self.target_freq = target_freq
-        self.tolerance = tolerance
-        self.num_bands = num_bands
         self.k_points = k_points
-        self.ensure_periodicity = ensure_periodicity
         self.geometry = geometry
         self.geometry_lattice = geometry_lattice
         self.geometry_center = geometry_center
         self.default_material = default_material
-        self.dimensions = dimensions
         self.random_fields = random_fields
         self.filename_prefix = filename_prefix
-        self.deterministic = deterministic
-        self.verbose = verbose
         self.optimize_grid_size = optimize_grid_size
-        self.eigensolver_nwork = eigensolver_nwork
-        self.eigensolver_block_size = eigensolver_block_size
         self.parity = ''
         self.iterations = 0
         self.all_freqs = None
         self.freqs = []
         self.band_range_data = []
-        self.eigensolver_flops = 0
         self.total_run_time = 0
         self.current_k = mp.Vector3()
         self.k_split_num = 1
         self.k_split_index = 0
         self.eigensolver_iters = []
-        self.mode_solver = None
+
+        grid_size = self._adjust_grid_size()
+
+        if type(self.default_material) is not mp.Medium and callable(self.default_material):
+            self.default_material.eps = False
+
+        self.mode_solver = mode_solver(
+            num_bands,
+            self.resolution,
+            self.geometry_lattice,
+            tolerance,
+            mesh_size,
+            self.default_material,
+            deterministic,
+            target_freq,
+            dimensions,
+            verbose,
+            ensure_periodicity,
+            eigensolver_flops,
+            is_negative_epsilon_ok,
+            epsilon_input_file,
+            mu_input_file,
+            force_mu,
+            use_simple_preconditioner,
+            grid_size,
+            eigensolver_nwork,
+            eigensolver_block_size,
+        )
+
+    @property
+    def num_bands(self):
+        return self.mode_solver.num_bands
+
+    @num_bands.setter
+    def num_bands(self, val):
+        self.mode_solver.set_num_bands(val)
 
     @property
     def resolution(self):
@@ -153,6 +170,151 @@ def resolution(self, val):
             t = type(val)
             raise TypeError("resolution must be a number or a Vector3: Got {}".format(t))
 
+        if self.mode_solver:
+            self.mode_solver.resolution = self._resolution
+            grid_size = self._adjust_grid_size()
+            self.mode_solver.set_grid_size(grid_size)
+
+    @property
+    def geometry_lattice(self):
+        return self._geometry_lattice
+
+    @geometry_lattice.setter
+    def geometry_lattice(self, val):
+        self._geometry_lattice = val
+        if self.mode_solver:
+            self.mode_solver.set_libctl_geometry_lattice(val)
+            grid_size = self._adjust_grid_size()
+            self.mode_solver.set_grid_size(grid_size)
+
+    @property
+    def tolerance(self):
+        return self.mode_solver.tolerance
+
+    @tolerance.setter
+    def tolerance(self, val):
+        self.mode_solver.tolerance = val
+
+    @property
+    def mesh_size(self):
+        return self.mode_solver.mesh_size
+
+    @mesh_size.setter
+    def mesh_size(self, val):
+        self.mode_solver.mesh_size = val
+
+    @property
+    def deterministic(self):
+        return self.mode_solver.deterministic
+
+    @deterministic.setter
+    def deterministic(self, val):
+        self.mode_solver.deterministic = val
+
+    @property
+    def target_freq(self):
+        return self.mode_solver.target_freq
+
+    @target_freq.setter
+    def target_freq(self, val):
+        self.mode_solver.target_freq = val
+
+    @property
+    def dimensions(self):
+        return self.mode_solver.get_libctl_dimensions()
+
+    @dimensions.setter
+    def dimensions(self, val):
+        self.mode_solver.set_libctl_dimensions(val)
+
+    @property
+    def verbose(self):
+        return self.mode_solver.verbose
+
+    @verbose.setter
+    def verbose(self, val):
+        self.mode_solver.verbose = val
+
+    @property
+    def ensure_periodicity(self):
+        return self.mode_solver.get_libctl_ensure_periodicity()
+
+    @ensure_periodicity.setter
+    def ensure_periodicity(self, val):
+        self.mode_solver.set_libctl_ensure_periodicity(val)
+
+    @property
+    def eigensolver_flops(self):
+        return self.mode_solver.eigensolver_flops
+
+    @eigensolver_flops.setter
+    def eigensolver_flops(self, val):
+        self.mode_solver.eigensolver_flops = val
+
+    @property
+    def is_negative_epsilon_ok(self):
+        return self.mode_solver.negative_epsilon_ok
+
+    @is_negative_epsilon_ok.setter
+    def is_negative_epsilon_ok(self, val):
+        self.mode_solver.negative_epsilon_ok = val
+
+    @property
+    def epsilon_input_file(self):
+        return self.mode_solver.epsilon_input_file
+
+    @epsilon_input_file.setter
+    def epsilon_input_file(self, val):
+        self.mode_solver.epsilon_input_file = val
+
+    @property
+    def mu_input_file(self):
+        return self.mode_solver.mu_input_file
+
+    @mu_input_file.setter
+    def mu_input_file(self, val):
+        self.mode_solver.mu_input_file = val
+
+    @property
+    def force_mu(self):
+        return self.mode_solver.force_mu
+
+    @force_mu.setter
+    def force_mu(self, val):
+        self.mode_solver.force_mu = val
+
+    @property
+    def use_simple_preconditioner(self):
+        return self.mode_solver.use_simple_preconditioner
+
+    @use_simple_preconditioner.setter
+    def use_simple_preconditioner(self, val):
+        self.mode_solver.use_simple_preconditioner = val
+
+    @property
+    def eigensolver_nwork(self):
+        return self.mode_solver.eigensolver_nwork
+
+    @eigensolver_nwork.setter
+    def eigensolver_nwork(self, val):
+        self.mode_solver.eigensolver_nwork = val
+
+    @property
+    def eigensolver_block_size(self):
+        return self.mode_solver.eigensolver_block_size
+
+    @eigensolver_block_size.setter
+    def eigensolver_block_size(self, val):
+        self.mode_solver.eigensolver_block_size = val
+
+    def _adjust_grid_size(self):
+        grid_size = self._get_grid_size()
+
+        if self.optimize_grid_size:
+            grid_size = self._optimize_grid_size(grid_size)
+
+        return grid_size
+
     def allow_negative_epsilon(self):
         self.is_negative_epsilon_ok = True
         self.target_freq = 1 / mp.inf
@@ -679,10 +841,11 @@ def _get_grid_size(self):
 
         return grid_size
 
-    def _optimize_grid_size(self):
-        self.grid_size.x = self.next_factor2357(self.grid_size.x)
-        self.grid_size.y = self.next_factor2357(self.grid_size.y)
-        self.grid_size.z = self.next_factor2357(self.grid_size.z)
+    def _optimize_grid_size(self, grid_size):
+        grid_size.x = self.next_factor2357(grid_size.x)
+        grid_size.y = self.next_factor2357(grid_size.y)
+        grid_size.z = self.next_factor2357(grid_size.z)
+        return grid_size
 
     def next_factor2357(self, n):
 
@@ -699,36 +862,7 @@ def divby(n, p):
         return self.next_factor2357(n + 1)
 
     def init_params(self, p, reset_fields):
-        self.grid_size = self._get_grid_size()
-
-        if self.optimize_grid_size:
-            self._optimize_grid_size()
-
-        self.mode_solver = mode_solver(
-            self.num_bands,
-            p,
-            self.resolution,
-            self.geometry_lattice,
-            self.tolerance,
-            self.mesh_size,
-            self.default_material,
-            self.geometry,
-            True if reset_fields else False,
-            self.deterministic,
-            self.target_freq,
-            self.dimensions,
-            self.verbose,
-            self.ensure_periodicity,
-            self.eigensolver_flops,
-            self.is_negative_epsilon_ok,
-            self.epsilon_input_file,
-            self.mu_input_file,
-            self.force_mu,
-            self.use_simple_preconditioner,
-            self.grid_size,
-            self.eigensolver_nwork,
-            self.eigensolver_block_size,
-        )
+        self.mode_solver.init(p, reset_fields, self.geometry, self.default_material)
 
     def set_parity(self, p):
         self.mode_solver.set_parity(p)
@@ -750,10 +884,6 @@ def run_parity(self, p, reset_fields, *band_functions):
         print("Initializing eigensolver data")
         print("Computing {} bands with {} tolerance".format(self.num_bands, self.tolerance))
 
-        if type(self.default_material) is not mp.Medium and callable(self.default_material):
-            # TODO: Support epsilon_function user materials like meep?
-            self.default_material.eps = False
-
         self.init_params(p, reset_fields)
 
         if isinstance(reset_fields, basestring):
@@ -846,6 +976,7 @@ def run_yodd_zodd(self, *band_functions):
 
     def find_k(self, p, omega, band_min, band_max, korig_and_kdir, tol,
                kmag_guess, kmag_min, kmag_max, *band_funcs):
+
         num_bands_save = self.num_bands
         kpoints_save = self.k_points
         nb = band_max - band_min + 1
diff --git a/python/tests/mpb.py b/python/tests/mpb.py
index 21a062287..78d62c903 100644
--- a/python/tests/mpb.py
+++ b/python/tests/mpb.py
@@ -69,6 +69,58 @@ def init_solver(self, geom=True):
             tolerance=1e-12
         )
 
+    def test_attribute_accessors(self):
+        ms = mpb.ModeSolver()
+        self.assertEqual(ms.mode_solver.num_bands, 1)
+        self.assertEqual(ms.mode_solver.deterministic, False)
+        self.assertEqual(ms.mode_solver.tolerance, 1.0e-7)
+        self.assertEqual(ms.mode_solver.mesh_size, 3)
+        self.assertEqual(ms.mode_solver.target_freq, 0)
+        self.assertEqual(ms.mode_solver.get_libctl_dimensions(), 3)
+        self.assertEqual(ms.mode_solver.verbose, False)
+        self.assertEqual(ms.mode_solver.get_libctl_ensure_periodicity(), True)
+        self.assertEqual(ms.mode_solver.eigensolver_flops, 0)
+        self.assertEqual(ms.mode_solver.negative_epsilon_ok, False)
+        self.assertEqual(ms.mode_solver.epsilon_input_file, '')
+        self.assertEqual(ms.mode_solver.mu_input_file, '')
+        self.assertEqual(ms.mode_solver.force_mu, False)
+        self.assertEqual(ms.mode_solver.use_simple_preconditioner, False)
+        self.assertEqual(ms.mode_solver.eigensolver_nwork, 3)
+        self.assertEqual(ms.mode_solver.eigensolver_block_size, -11)
+
+        ms.num_bands = 2
+        ms.deterministic = True
+        ms.tolerance = 1.0e-12
+        ms.mesh_size = 2
+        ms.target_freq = 1
+        ms.dimensions = 2
+        ms.verbose = True
+        ms.ensure_periodicity = False
+        ms.eigensolver_flops = 20
+        ms.is_negative_epsilon_ok = True
+        ms.epsilon_input_file = 'test'
+        ms.mu_input_file = 'test'
+        ms.force_mu = True
+        ms.use_simple_preconditioner = True
+        ms.eigensolver_nwork = 4
+        ms.eigensolver_block_size = 11
+        self.assertEqual(ms.mode_solver.num_bands, 2)
+        self.assertEqual(ms.mode_solver.deterministic, True)
+        self.assertEqual(ms.mode_solver.tolerance, 1.0e-12)
+        self.assertEqual(ms.mode_solver.mesh_size, 2)
+        self.assertEqual(ms.mode_solver.target_freq, 1)
+        self.assertEqual(ms.mode_solver.get_libctl_dimensions(), 2)
+        self.assertEqual(ms.mode_solver.verbose, True)
+        self.assertEqual(ms.mode_solver.get_libctl_ensure_periodicity(), False)
+        self.assertEqual(ms.mode_solver.eigensolver_flops, 20)
+        self.assertEqual(ms.mode_solver.negative_epsilon_ok, True)
+        self.assertEqual(ms.mode_solver.epsilon_input_file, 'test')
+        self.assertEqual(ms.mode_solver.mu_input_file, 'test')
+        self.assertEqual(ms.mode_solver.force_mu, True)
+        self.assertEqual(ms.mode_solver.use_simple_preconditioner, True)
+        self.assertEqual(ms.mode_solver.eigensolver_nwork, 4)
+        self.assertEqual(ms.mode_solver.eigensolver_block_size, 11)
+
     def test_resolution(self):
         ms = self.init_solver()
         self.assertEqual([32, 32, 32], ms.resolution)
@@ -1349,7 +1401,7 @@ def test_poynting(self):
     def test_fractional_lattice(self):
 
         ms = self.init_solver()
-        ms.geometry_lattice.size = mp.Vector3(0.1, 0.1)
+        ms.geometry_lattice = mp.Lattice(size=mp.Vector3(0.1, 0.1))
         ms.run_te()
 
         expected_brd = [