diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md
index 66a8f443d..8d6003da1 100644
--- a/doc/docs/Python_User_Interface.md
+++ b/doc/docs/Python_User_Interface.md
@@ -1072,7 +1072,7 @@ As described in the tutorial, you normally use `add_flux` via statements like:
 
 **`transmission = sim.add_flux(...)`**
 —
-to store the flux object in a variable. `add_flux` initializes the fields if necessary, just like calling `run`, so you should only call it *after* initializing your `Simulation` object which includes specifying `geometry`, `sources`, `boundary_layers`, etcetera. You can create as many flux objects as you want, e.g. to look at powers flowing in different regions or in different frequency ranges. Note, however, that Meep has to store (and update at every time step) a number of Fourier components equal to the number of grid points intersecting the flux region multiplied by the number of electric and magnetic field components required to get the Poynting vector multiplied by `nfreq`, so this can get quite expensive (in both memory and time) if you want a lot of frequency points over large regions of space.
+to store the flux object in a variable.  You can create as many flux objects as you want, e.g. to look at powers flowing in different regions or in different frequency ranges. Note, however, that Meep has to store (and update at every time step) a number of Fourier components equal to the number of grid points intersecting the flux region multiplied by the number of electric and magnetic field components required to get the Poynting vector multiplied by `nfreq`, so this can get quite expensive (in both memory and time) if you want a lot of frequency points over large regions of space.
 
 Once you have called `add_flux`, the Fourier transforms of the fields are accumulated automatically during time-stepping by the [run functions](#run-functions). At any time, you can ask for Meep to print out the current flux spectrum via:
 
@@ -1205,7 +1205,7 @@ As for energy regions, you normally use `add_energy` via statements like:
 En = sim.add_energy(...)
 ```
 
-to store the energy object in a variable. `add_energy` initializes the fields if necessary, just like calling `run`, so you should only call it *after* setting up your `geometry`, `sources`, `boundary_layers`, etcetera. You can create as many energy objects as you want, e.g. to look at the energy densities in different objects or in different frequency ranges. Note, however, that Meep has to store (and update at every time step) a number of Fourier components equal to the number of grid points intersecting the energy region multiplied by `nfreq`, so this can get quite expensive (in both memory and time) if you want a lot of frequency points over large regions of space.
+to store the energy object in a variable.  You can create as many energy objects as you want, e.g. to look at the energy densities in different objects or in different frequency ranges. Note, however, that Meep has to store (and update at every time step) a number of Fourier components equal to the number of grid points intersecting the energy region multiplied by `nfreq`, so this can get quite expensive (in both memory and time) if you want a lot of frequency points over large regions of space.
 
 Once you have called `add_energy`, the Fourier transforms of the fields are accumulated automatically during time-stepping by the `run` functions. At any time, you can ask for Meep to print out the current energy density spectrum via:
 
@@ -1287,7 +1287,7 @@ As for force regions, you normally use `add_force` via statements like:
 Fx = sim.add_force(...)
 ```
 
-to store the force object in a variable. `add_force` initializes the fields if necessary, just like calling `run`, so you should only call it *after* initializing your `Simulation` object which includes specifying `geometry`, `sources`, `boundary_layers`, etcetera. You can create as many force objects as you want, e.g. to look at forces on different objects, in different directions, or in different frequency ranges. Note, however, that Meep has to store (and update at every time step) a number of Fourier components equal to the number of grid points intersecting the force region, multiplied by the number of electric and magnetic field components required to get the stress vector, multiplied by `nfreq`, so this can get quite expensive (in both memory and time) if you want a lot of frequency points over large regions of space.
+to store the force object in a variable.  You can create as many force objects as you want, e.g. to look at forces on different objects, in different directions, or in different frequency ranges. Note, however, that Meep has to store (and update at every time step) a number of Fourier components equal to the number of grid points intersecting the force region, multiplied by the number of electric and magnetic field components required to get the stress vector, multiplied by `nfreq`, so this can get quite expensive (in both memory and time) if you want a lot of frequency points over large regions of space.
 
 Once you have called `add_force`, the Fourier transforms of the fields are accumulated automatically during time-stepping by the `run` functions. At any time, you can ask for Meep to print out the current force spectrum via:
 
diff --git a/doc/docs/Scheme_User_Interface.md b/doc/docs/Scheme_User_Interface.md
index 8aca81d1b..03df55aca 100644
--- a/doc/docs/Scheme_User_Interface.md
+++ b/doc/docs/Scheme_User_Interface.md
@@ -904,7 +904,7 @@ As described in the tutorial, you normally use `add-flux` via statements like:
 
 **`(define transmission (add-flux ...))`**
 —
-to store the flux object in a variable. `add-flux` initializes the fields if necessary, just like calling `run`, so you should only call it *after* setting up your `geometry`, `sources`, `pml-layers`, etcetera. You can create as many flux objects as you want, e.g. to look at powers flowing in different regions or in different frequency ranges. Note, however, that Meep has to store (and update at every time step) a number of Fourier components equal to the number of grid points intersecting the flux region multiplied by the number of electric and magnetic field components required to get the Poynting vector multiplied by `nfreq`, so this can get quite expensive (in both memory and time) if you want a lot of frequency points over large regions of space.
+to store the flux object in a variable. `add-flux` initializes the fields if necessary, just like calling `run`, so you should only call it *after* setting up your `geometry`, `sources`, `pml-layers`, `k-point`, etcetera. You can create as many flux objects as you want, e.g. to look at powers flowing in different regions or in different frequency ranges. Note, however, that Meep has to store (and update at every time step) a number of Fourier components equal to the number of grid points intersecting the flux region multiplied by the number of electric and magnetic field components required to get the Poynting vector multiplied by `nfreq`, so this can get quite expensive (in both memory and time) if you want a lot of frequency points over large regions of space.
 
 Once you have called `add-flux`, the Fourier transforms of the fields are accumulated automatically during time-stepping by the [run functions](#run-functions). At any time, you can ask for Meep to print out the current flux spectrum via:
 
@@ -1005,7 +1005,7 @@ As for energy regions, you normally use `add-energy` via statements like:
 (define En (add-energy ...))
 ```
 
-to store the energy object in a variable. `add-energy` initializes the fields if necessary, just like calling `run`, so you should only call it *after* setting up your `geometry`, `sources`, `pml-layers`, etcetera. You can create as many energy objects as you want, e.g. to look at the energy densities in different objects or in different frequency ranges. Note, however, that Meep has to store (and update at every time step) a number of Fourier components equal to the number of grid points intersecting the energy region multiplied by `nfreq`, so this can get quite expensive (in both memory and time) if you want a lot of frequency points over large regions of space.
+to store the energy object in a variable. `add-energy` initializes the fields if necessary, just like calling `run`, so you should only call it *after* setting up your `geometry`, `sources`, `pml-layers`, `k-point`, etcetera. You can create as many energy objects as you want, e.g. to look at the energy densities in different objects or in different frequency ranges. Note, however, that Meep has to store (and update at every time step) a number of Fourier components equal to the number of grid points intersecting the energy region multiplied by `nfreq`, so this can get quite expensive (in both memory and time) if you want a lot of frequency points over large regions of space.
 
 Once you have called `add-energy`, the Fourier transforms of the fields are accumulated automatically during time-stepping by the `run` functions. At any time, you can ask for Meep to print out the current energy density spectrum via:
 
diff --git a/python/simulation.py b/python/simulation.py
index 83eb13053..fa14f3d10 100644
--- a/python/simulation.py
+++ b/python/simulation.py
@@ -1255,11 +1255,6 @@ def run_k_points(self, t, k_points):
 
         for k in k_points:
             k_index += 1
-
-            if k_index == 1:
-                self.init_sim()
-                output_epsilon(self)
-
             harminv = self.run_k_point(t, k)
             freqs = [complex(m.freq, m.decay) for m in harminv.modes]
 
diff --git a/scheme/meep.scm.in b/scheme/meep.scm.in
index fee50d199..4fbb632cc 100644
--- a/scheme/meep.scm.in
+++ b/scheme/meep.scm.in
@@ -1353,7 +1353,7 @@
   (define all-freqs '())
   (map (lambda (k)
 	 (set! k-index (+ k-index 1))
-	 (if (= k-index 1) (begin (init-fields) (output-epsilon)))
+	 (if (= k-index 1) (begin (change-k-point! k) (if (null? fields) (init-fields)) (output-epsilon)))
 	 (let ((freqs (run-k-point T k)))
 	   (print "freqs:, " k-index
 		  ", " (vector3-x k) ", " (vector3-y k) ", " (vector3-z k))