anchor.diff:    total_rate = rate_release_s + rate_release_p
     enzyme_capture_radius = PTread_d(config, 'enzyme_capture_radius', loc=params_group)
     substrate_fraction = PTread_d(config, 'substrate_fraction', loc=params_group)
  +  product_relative_fraction = PTread_d(config, 'product_relative_fraction', loc=params_group)
   
     L = PTread_ivec(config, 'L', 3, loc=params_group)
     rho = PTread_i(config, 'rho', loc=params_group)
  
positive.code:      solvent% vel(3,k) = threefry_normal(state(1))*sqrt(T)
       if (threefry_double(state(1)) < substrate_fraction) then
          ! substrate
          solvent%species(k) = 1
       else
          ! neutral species
          solvent%species(k) = 3
  
negative.code: 
    T = PTread_d(config, 'T', loc=params_group)
    link_d = PTread_dvec(config, 'link_d', 2, loc=params_group)
    link_angles = PTread_dvec(config, 'link_angles', 2, loc=params_group)
    link_angle = link_angles(1)
    elastic_k = PTread_d(config, 'elastic_k', loc=params_group)
    angle_k = PTread_d(config, 'angle_k', loc=params_group)
  
  