First of all, we parametrize the reaction rate and the Q-value (for the back-reaction) in the corresponding subroutine in the main code file:

      subroutine rate_png(temp,den,fr,dfrdt,dfrdd,rr,drrdt,drrdd)
      include 'implno.dek'
      include 'tfactors.dek'
      include 'helper_deuteronrate.dek'
 
! declare the pass
      double precision temp,den,fr,dfrdt,dfrdd,rr,drrdt,drrdd
 
! locals
      double precision term,dtermdt,rev,drevdt,aa,daa
 
 
! p(n,g)d
! smith,kawano,malany 1992
 
       aa      = 1.0d0 - npgd_a1*t912 + npgd_a2*t9 &
                 - npgd_a3*t932 + npgd_a4*t92 &
                 - npgd_a5*t952
 
       daa     =  -0.5d0*npgd_a1*t9i12 + npgd_a2 &
                 - 1.5d0*npgd_a3*t912 + 2.0d0*npgd_a4*t9 &
                 - 2.5d0*npgd_a5*t932
 
       term    = npgd_a0 * aa
       dtermdt = npgd_a0 * daa
 
 
! rates
      fr    = den * term
      dfrdt = den * dtermdt * 1.0d-9
      dfrdd = term
 
      rev      = 4.71d+09 * t932 * exp(-npgd_Eb*t9i)
      drevdt   = rev*(1.5d0*t9i + npgd_Eb*t9i2)
 
      rr    = rev * term
      drrdt = (drevdt * term + rev * dtermdt) * 1.0d-9
      drrdd = 0.0d0
 
      return
      end

These parameters should be read in in another subroutine, so they have to be global variables. Hence we create a helper file helper_deuteronrate.dek:

! helper file to create global variables for the n(p,g)d reaction
 
	real*8 :: npgd_Eb, npgd_a0, npgd_a1, npgd_a2, npgd_a3, npgd_a4, npgd_a5, npgd_a6, npgd_a7 
	common /helpdeut/ npgd_Eb, npgd_a0, npgd_a1, npgd_a2, npgd_a3, npgd_a4, npgd_a5, npgd_a6, npgd_a7 

This helper file needs of course to be added to all subroutines working with these parameters.

We then read in these variables in the net_input subroutine after querying the cosmological variables:

      write(6,02) 'enter a0 to a7, mn and mp => '
      read(5,*) npgd_Eb, npgd_a0, npgd_a1, npgd_a2, npgd_a3, npgd_a4, npgd_a5, npgd_a6, npgd_a7, mn, mp

We then create a wrapper script to read in the parameters from a file, run the BBN code and save the final composition of each run to a new file:

#!/bin/bash 
 
# wrapper script to batch-feed mn, mp and the p(n,g)d rate coefficients into bbn code
 
### config section ###
prefix="test_"
endtime=1e6
inputfile="coefficients.in" #input file containing a column used as label as well as a0-a7, mn and mp for the p(n,g)d reaction rate
variablelabel="Lambda_QCD" #just a label used in the final.dat file
 
rm $prefix*
n=0
while read line; do
  if [ ${line:0:1} != "#" ]; then #skip comment lines
    prefixn=$prefix$n"_"
 
    var=`echo $line | awk '{ print $1; }'`
    params=`echo $line | cut -f2- -d ' '`
    echo -e "\n$params\n$endtime\n$prefixn\n1\n" | ./bigbang_deuteron.out
 
    #copy header for first file
    if [ $n -eq 0 ]; then
      head -n 2 $prefixn'0002_z1.dat' > $prefix"final.dat"
      echo " * variable varied: $variablelabel" >> $prefix"final.dat"
      awk '{$1=$2="";print;}' $prefixn'0003_z1.dat'  | paste $prefixn'0002_z1.dat' - | sed '3q;d' >> $prefix"final.dat"
    fi
 
    result1=`tail -n 1 $prefixn'0002_z1.dat'`
    result2=`tail -n 1 $prefixn'0003_z1.dat' | awk '{$1=$2="";print;}' # | cut -f3-`
    echo $var" "$result1" "$result2 >> $prefix"final.dat"
    n=$((n+1))
  fi
done < $inputfile

This script expects an $inputfile as follows (made-up data for now!)

#Lambda	Eb	a0	a1	a2	a3	a4	a5	a6	a7	mn	mp
217	25.82	4.742e4	0.8504	0.4895	0.09623	8.471e-3	2.80e-4	0	0	1.67492721184d-24	1.67262163783d-24
217	27	4.742e4	0.8504	0.4895	0.09623	8.471e-3	2.80e-4	0	0	1.67492721184d-24	1.67262163783d-24
220	27	7.342e4	0.8504	0.4895	0.09623	8.471e-3	2.80e-4	0	0	1.67492721184d-24	1.67492721184d-24

Now, “all” which remains is to figure out how the binding energy, the a's and the nucleon masses depend on $\Lambda_{QCD}$ and to calculate the above table.

The complete BBN code with our changes can be found here.

The rates and other input variables have been calculated outside the BBN code.