diff --git a/source/direct_ion.c b/source/direct_ion.c
index 58ed1970f..38696b7cd 100644
--- a/source/direct_ion.c
+++ b/source/direct_ion.c
@@ -89,8 +89,13 @@ compute_qrecomb_coeffs(T)
   int n, nvmin, ntmin;
   struct topbase_phot *xtop;
 
-  for (n = 0; n < nions; n++)
+  for (n = 0; n < nions; n++) //We need to generate data for the ions doing the recombining.
+    {
+      /* There is only any point doing this is we are not a neutral ion */
+
+      if (ion[n].istate > 1)  //There is only any point doing this is we are not a neutral ion
     {
+
       if (ion[n].dere_di_flag == 0)
   {
     //printf ("NO COEFFS\n");
@@ -100,24 +105,30 @@ compute_qrecomb_coeffs(T)
       else
   {
     /* we need to know about the bound-free jump, so we need the details
-       from the ground state cross-section for this ion */
+       from the ground state cross-section for for the ion below this one */
 
-    ntmin = ion[n].ntop_ground;  /* We only ever use the ground state cont_ptr. 
-                                       This is for topbase */
-    nvmin = ion[n].nxphot;
+    ntmin = ion[n-1].ntop_ground;  /* We only ever use the ground state cont_ptr. 
+                                      This is for topbase */
+    nvmin = ion[n-1].nxphot;
 
-    if (ion[n].phot_info > 0)  //topbase or hybrid VFKY (GS)+TB excited
+    if (ion[n-1].phot_info > 0)  //topbase or hybrid VFKY (GS)+TB excited
       { 
         xtop = &phot_top[ntmin];
       }
-    else if (ion[n].phot_info == 0)  // verner
+    else if (ion[n-1].phot_info == 0)  // verner
       {   //just the ground state ionization fraction.
         xtop = &phot_top[nvmin];
       }
-
+    else
+      {
+        Error ("compute_qrecomb_coeffs: no coll ionization data for ion %i\n",n-1);
+      }
+    
+    /* this will return 0 if there aren't any coeffs */
     qrecomb_coeffs[n] = q_recomb_dere(xtop, T);
   }
 
+    }   //End of if statement for neutral ions
     }   //End of loop over ions
 
   return (0);
@@ -348,15 +359,28 @@ q_recomb_dere (cont_ptr, electron_temperature)
     {
       root_etemp = sqrt(electron_temperature);     
       coeff = 2.07e-16 / (root_etemp * root_etemp * root_etemp);
-      coeff *= config[cont_ptr->nlev].g / config[cont_ptr->uplev].g;
+
+      /* the original way of getting mutilplicity doesn't work for non-matoms, 
+         because uplevel isn't identified */
+      //coeff *= config[cont_ptr->nlev].g / config[cont_ptr->uplev].g;
+
+      /* JM/NSH XXX -- This is the multiplicity of the ground states. 
+         Should be of order 1 so it may be better to just leave it out, 
+         since the collisional ionization cross section is doubtless 
+         averaged over upper states... */   
+      coeff *= ion[nion].g / ion[nion+1].g;   
+  
       coeff *= exp(u0);
       coeff *= q_ioniz_dere(nion, electron_temperature);
 
 
       /* do a sane check here, as there's an exponential which could blow up */
       if (sane_check(coeff))
-        Error("q_recomb is %8.4e for ion %i at temperature %8.4e\n",
+      {
+        Error("q_recomb is %8.4e for ion %i at temperature %8.4e, setting to zero\n",
                coeff, nion, electron_temperature);
+        coeff = 0.0;
+      }
     }
   else 
     coeff = 0.0;
@@ -374,7 +398,7 @@ JM 1301 -- Edited this to avoid need to call q_ioniz and exponential.
 Should improve speed and stability
 
 This equation comes from considering TE and getting the expression
-q_recomb = 2.07e-16 * gl/gu * exp(E/kT) * q_ioniz
+q_recomb = 2.07e-16 * gl/gu * exp(E/kT) * (T_e**-3/2) * q_ioniz
 then substituting the above expression for q_ioniz.
 */
 
diff --git a/source/para_update.c b/source/para_update.c
index 174fcc75d..52e27ac53 100644
--- a/source/para_update.c
+++ b/source/para_update.c
@@ -164,7 +164,7 @@ communicate_estimators_para ()
       iredhelper[mpi_i + 4 * NPLASMA] = plasmamain[mpi_i].ntot_wind;
       iredhelper[mpi_i + 5 * NPLASMA] = plasmamain[mpi_i].ntot_agn;
       iredhelper[mpi_i + 6 * NPLASMA] = plasmamain[mpi_i].nioniz;
-	  
+
       for (mpi_j = 0; mpi_j < NXBANDS; mpi_j++)
 	{
 	  iredhelper[mpi_i + (7 + mpi_j) * NPLASMA] = plasmamain[mpi_i].nxtot[mpi_j];
@@ -185,7 +185,8 @@ communicate_estimators_para ()
       plasmamain[mpi_i].ntot_disk = iredhelper2[mpi_i + 3 * NPLASMA];
       plasmamain[mpi_i].ntot_wind = iredhelper2[mpi_i + 4 * NPLASMA];
       plasmamain[mpi_i].ntot_agn = iredhelper2[mpi_i + 5 * NPLASMA];
-	  plasmamain[mpi_i].nioniz = iredhelper2[mpi_i + 6 * NPLASMA];
+	    plasmamain[mpi_i].nioniz = iredhelper2[mpi_i + 6 * NPLASMA];
+
       for (mpi_j = 0; mpi_j < NXBANDS; mpi_j++)
 	{
 	  plasmamain[mpi_i].nxtot[mpi_j] = iredhelper2[mpi_i + (7 + mpi_j) * NPLASMA];