Fix rate tables for r_li7_pa_he4

Works on #497

The problem occurs when using a custom rate table for for r_li7_pa_he4
in that we could not find the rate id in the reaclib data. During
normal operation we dont look up the rate in reaclib as it has a hard coded
rate id and we use a function call in set_raw_rate to get the rate value, rather
than use reaclib.

But when we load the rate from a table we do need to do a lookup on reaclib
to know the forward and reverse rate id. But the code for working out the
handle for a rate was assuming this rate was still r_li7_p_to_he4_he4 as we had
two identical outputs, rather than marking it as an _ap reaction. Thus we
could never find the rate in reaclib.

Reported-by Hannah Brinkman

Author: rjfarmer

rates/private/raw_rates.f90

@@ -909,7 +909,10 @@ subroutine eval_table_reverse(ir, rir, tf, temp, fr, rr, ierr)
          dlambda_dlnT = 0d0
 
          call do_reaclib_indices_for_reaction(reaction_name(rir), reaclib_rates, lo, hi, ierr)
-         if(ierr/=0) return 
+         if(ierr/=0) then
+            write(*,*) "Error: Could not get reaclib index for rate ",rir, ' ',trim(reaction_name(rir))
+            return 
+         end if
 
          inv_lambda = 0d0
          dinv_lambda_dlnT = 0d0

rates/private/reaclib_support.f90

@@ -636,8 +636,22 @@ end function two_one
 
          logical function two_two()
             two_two = .true.
-            if (in1 == 0 .or. in2 == 0 .or. out1 == 0 .or. out2 == 0 .or. &
-                in1 == in2 .or. out1 == out2) return
+
+            if (in1 == 0 .or. in2 == 0 .or. out1 == 0 .or. out2 == 0) return
+            
+            ! Special case r_li7_pa_he4, this must come first otherwise the out1==out2
+            ! check will label this rate as a _to_ reaction instead of an _ap reaction
+            if (nuclides% Z(in1) == 1 .and. nuclides% N(in1) == 0 .and. &
+               nuclides% Z(out1) == 2 .and. nuclides% N(out1) == 2 .and. &
+               nuclides% Z(out2) == 2 .and.  nuclides% N(out2) == 2 .and. &
+               nuclides% Z(in2) == 3 .and. nuclides% N(in2) == 4) then
+               handle = trim(handle) // '_pa'
+               two_two = .false.
+               return
+            end if
+
+            if(in1==in2 .or. out1==out2) return
+            
             two_two = .false.
             if (nuclides% Z(in1) == 2 .and. nuclides% N(in1) == 2 .and. &
                      nuclides% Z(out1) == 1 .and. nuclides% N(out1) == 0 .and. &