diff --git a/pyradex/base_class.py b/pyradex/base_class.py
index 3f026ed..1374ac0 100644
--- a/pyradex/base_class.py
+++ b/pyradex/base_class.py
@@ -257,6 +257,11 @@ def _validate_colliders(self):
                 bad_colliders.append(collider)
                 OK = False
 
+        if 'h2' in valid_colliders and density['H2'] == 0 and (density['oH2'] != 0
+                                                               or
+                                                               density['pH2'] != 0):
+            OK = True
+
         if not OK:
             raise ValueError("There are colliders with specified densities >0 "
                              "that do not have corresponding collision rates."
diff --git a/pyradex/core.py b/pyradex/core.py
index c5eed62..33ee8a5 100644
--- a/pyradex/core.py
+++ b/pyradex/core.py
@@ -24,6 +24,16 @@
 __all__ = ['pyradex', 'write_input', 'parse_outfile', 'call_radex', 'Radex',
            'density_distribution']
 
+
+def _init_radex(radex):
+    """
+    Initialize a RADEX object to blank values so it doesn't trigger any errors.
+    The goal of this method is to re-initialize the radex.so fortran objects
+    and remove any values stored in the fortran-wrapped structures, since these
+    values will persist between independent python objects in unclear and
+    unpredictable ways.
+    """
+    radex.cphys.density[:] = 0
         
 
 def pyradex(executable='radex', minfreq=100, maxfreq=130,
@@ -279,6 +289,11 @@ def __init__(self,
         from pyradex.radex import radex
         self.radex = radex
 
+        # the 'radex' object is stateful; it needs to be reset
+        # (this is bad, but a necessary hack to deal with the underlying
+        # fortran objects being wrapped)
+        _init_radex(radex)
+
         self.mu = mu
 
         if os.getenv('RADEX_DATAPATH') and datapath is None:
@@ -509,6 +524,8 @@ def density(self, collider_density):
                 raise ValueError('Collider %s is not one of the valid colliders: %s' %
                                  (k,self._all_valid_colliders))
 
+        valid_collider_lowercase = [x.lower() for x in self.valid_colliders]
+
         if (('OH2' in collider_densities and collider_densities['OH2'] !=0) or
             ('PH2' in collider_densities and collider_densities['PH2'] !=0)):
 
@@ -526,9 +543,10 @@ def density(self, collider_density):
                 self.radex.cphys.density[2] = collider_densities['OH2']
             self._use_thermal_opr = False
         elif 'H2' in collider_densities:
-            warnings.warn("Using a default ortho-to-para ratio (which "
-                          "will only affect species for which independent "
-                          "ortho & para collision rates are given)")
+            if 'h2' not in valid_collider_lowercase:
+                warnings.warn("Using a default ortho-to-para ratio (which "
+                              "will only affect species for which independent "
+                              "ortho & para collision rates are given)")
             self._use_thermal_opr = True
             #self.radex.cphys.density[0] = collider_densities['H2']
 
@@ -545,12 +563,11 @@ def density(self, collider_density):
 
         # RADEX relies on n(H2) = n(oH2) + n(pH2)
         # We have set n(oH2) and n(pH2) above
-        vc = [x.lower() for x in self.valid_colliders]
-        if 'h2' in vc:
+        if 'h2' in valid_collider_lowercase:
             self.radex.cphys.density[0] = self.radex.cphys.density[1:3].sum()
             self.radex.cphys.density[1] = 0
             self.radex.cphys.density[2] = 0
-        elif 'oh2' in vc or 'ph2' in vc:
+        elif 'oh2' in valid_collider_lowercase or 'ph2' in valid_collider_lowercase:
             self.radex.cphys.density[0] = 0
 
         self.radex.cphys.density[3] = collider_densities['E']
@@ -746,8 +763,17 @@ def temperature(self, tkin):
         if self._use_thermal_opr:
             # Reset the density to a thermal value
             lp = self._locked_parameter
-            self.density = (unitless(self.density['H2']) or
-                            unitless(self.density['oH2']+self.density['pH2']))
+            # In order to preserve secondary colliders (He, e-, etc), we first load the current
+            # density, then we set the total H2 density to be oH2+pH2
+            # Within the self.density setter, radex.readdata will be called, and it will re-compute
+            # the appropriate O/P density given the newly assigned temperature
+            # (we have to make the density modifiable first)
+            dens_to_set = dict(self.density)
+            if dens_to_set['H2'] == 0:
+                dens_to_set['H2'] = unitless(self.density['oH2']+self.density['pH2'])
+                dens_to_set['oH2'] = 0
+                dens_to_set['pH2'] = 0
+            self.density = dens_to_set
             self._locked_parameter = lp
 
     @property