diff --git a/docs/index.rst b/docs/index.rst
index 92a7444..6139688 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -6,8 +6,7 @@ Population Synthesis (FSPS) Fortran library
 <https://github.com/cconroy20/fsps>`_.
 Python-FSPS makes it easy to generate spectra and magnitudes for arbitrary
 stellar populations.
-These bindings are updated in conjunction with FSPS, and are known to work with
-FSPS v3.0.
+These bindings are updated in conjunction with FSPS.
 
 
 User Guide
diff --git a/docs/installation.rst b/docs/installation.rst
index 5748883..6fa4bf6 100644
--- a/docs/installation.rst
+++ b/docs/installation.rst
@@ -21,7 +21,7 @@ fail to import if it is set incorrectly.
 
 Python-FSPS is built against specific versions of the FSPS Fortran API and data
 files, so it is important that you have a recent version of FSPS through git.
-Currently Python-FSPS is built against FSPS v3.2.
+Currently Python-FSPS is built against FSPS v4.0.
 
 Installing stable version
 -------------------------
diff --git a/src/fsps/fsps.py b/src/fsps/fsps.py
index 27fbeaa..74c2f96 100644
--- a/src/fsps/fsps.py
+++ b/src/fsps/fsps.py
@@ -46,20 +46,11 @@ class StellarPopulation(object):
         populations (SSPs) is performed before computing composite stellar
         population (CSP) models:
 
-        * 0: No interpolation, use the metallicity index specified by ``zmet``.
-        * 1: The SSPs are interpolated to the value of ``logzsol`` before the
-          spectra and magnitudes are computed, and the value of ``zmet`` is
-          ignored.
-        * 2: The SSPs are convolved with a metallicity distribution function
-          specified by the ``logzsol`` and ``pmetals`` parameters. The value of
-          ``zmet`` is ignored.
-        * 3: Use all available SSP metallicities when computing the composite
-          model, for use exclusively with tabular SFHs where the metallicity
-          evolution as function of age is given (see `set_tabular_sfh()`).  The
-          values of ``zmet`` and ``logzsol`` are ignored.  Furthermore
-          ``add_neb_emission`` must be set to False.
-
-        Can only be changed during initialization.
+        * 0: No interpolation, use the metallicity index specified by ``zmet``
+          (and ``afeindx``).
+        * 1: The SSPs are interpolated to the value of ``logzsol`` and ``afe``
+          before the spectra and magnitudes are computed, and the values of
+          ``zmet`` and ``afeindx`` are ignored.
 
     :param add_agb_dust_model: (default: True)
         Switch to turn on/off the AGB circumstellar dust model presented in
@@ -1409,17 +1400,22 @@ def check_params(self):
         NZ = driver.get_nz()
         assert self._params["zmet"] in range(
             1, NZ + 1
-        ), "zmet={0} out of range [1, {1}]".format(self._params["zmet"], NZ)
+        ), f"zmet={self._params['zmet']} out of range [1, {NZ}]"
+        NAFE = driver.get_nafe()
+        assert self._params["afeindx"] in range(
+            1, NAFE + 1
+        ), f"afeindx={self._params['afeindx']} out of range [1, {NAFE}]"
         assert self._params["dust_type"] in range(
             7
-        ), "dust_type={0} out of range [0, 6]".format(self._params["dust_type"])
+        ), f"dust_type={self._params['dust_type']} out of range [0, 6]"
         assert self._params["imf_type"] in range(
             6
-        ), "imf_type={0} out of range [0, 5]".format(self._params["imf_type"])
+        ), f"imf_type={self._params['imf_type']} out of range [0, 5]"
         assert (self._params["tage"] <= 0) | (
             self._params["tage"] > self._params["sf_start"]
-        ), "sf_start={0} is greater than tage={1}".format(
-            self._params["sf_start"], self._params["tage"]
+        ), (
+            f"sf_start={self._params['sf_start']} is greater than "
+            f"tage={self._params['tage']}"
         )
         assert (
             self._params["const"] + self._params["fburst"]
diff --git a/tests/tests.py b/tests/tests.py
index 6242248..10168e9 100644
--- a/tests/tests.py
+++ b/tests/tests.py
@@ -91,6 +91,18 @@ def test_param_checks(pop_and_params):
     pop.params["sf_start"] = 0.1
     w, s = pop.get_spectrum(tage=pop.params["tage"])
 
+    pop.params["zmet"] = 100
+    with pytest.raises(AssertionError):
+        w, s = pop.get_spectrum(tage=pop.params["tage"])
+    pop.params["zmet"] = 1
+    w, s = pop.get_spectrum(tage=pop.params["tage"])
+
+    pop.params["afeindx"] = pop.n_afe + 1
+    with pytest.raises(AssertionError):
+        w, s = pop.get_spectrum(tage=pop.params["tage"])
+    pop.params["afeindx"] = 1
+    w, s = pop.get_spectrum(tage=pop.params["tage"])
+
 
 def test_smooth_lsf(pop_and_params):
     # recomputes SSPs