anchor.diff:          end do
        end do
     end do
  -  
  +
     call n_solvent_el%create_time(hfile%observables, 'n_solvent', &
          n_solvent, ior(H5MD_LINEAR, H5MD_STORE_TIME), step=N_MD_steps, &
          time=N_MD_steps*dt)
  
positive.code:     bound_molecule_id(enzyme_idx) = solvent%id(idx)
      solvent%species(idx) = 0
  
      en1 = compute_bead_energy(enzyme_idx)
  
      link_angle(enzyme_idx) = link_angles(2)
  
      en2 = compute_bead_energy(enzyme_idx)
  
      call add_energy_to_cell(cell_idx, excess_kinetic_energy + en1 - en2)
  
negative.code: 
  
    !> Bind solvent molecule idx to enzyme
    subroutine bind_molecule(enzyme_idx, idx)
      integer, intent(in) :: enzyme_idx
      integer, intent(in) :: idx
  