anchor.diff:       call compute_cell_wise_max_v
   
        call bulk_reac_timer%tic()
  -     call bulk_reaction_endothermic(solvent, solvent_cells, 2, 1, bulk_rate, tau, state, delta_u)
  +     call bulk_reaction_endothermic(solvent, solvent_cells, 2, 1, bulk_rate(2), tau, state, delta_u)
  +     call bulk_reaction_endothermic(solvent, solvent_cells, 1, 2, bulk_rate(1), tau, state, -delta_u)
        call bulk_reac_timer%tac()
   
        if (sampling) then
  
positive.code:   end do
    solvent% vel = solvent% vel - spread(sum(solvent% vel, dim=2)/solvent% Nmax, 2, solvent% Nmax)
    solvent% force = 0
    solvent% species = 1
  
    call n_solvent_el%create_time(hfile%observables, 'n_solvent', &
         n_solvent, H5MD_LINEAR, step=N_MD_steps, &
  
negative.code:   double precision :: quaternion_treshold
    character(len=144) :: fluid_wall
    double precision :: delta_u
    double precision :: catalytic_probability(N_species)
    double precision :: bulk_rate(N_species)
  
    integer, parameter :: block_length = 8
    type(axial_correlator_t) :: axial_cf
  