anchor.diff:      ! Use c.o.m. velocity so that no kinetic energy exchange must take place
       solvent%vel(:,i) = colloids%vel(:,enz_2)
   
  +    ! compute conformational energy difference when changing the angle
       en1 = compute_bead_energy(enzyme_idx)
  -
       link_angle(enzyme_idx) = link_angles(1)
  -
       en2 = compute_bead_energy(enzyme_idx)
   
       p = solvent_cells%cartesian_indices(solvent%pos(:, i))
  
positive.code:         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)
  
negative.code:         do k = 1, colloids%Nmax
             call compute_radial_histogram(radial_hist(k), colloids%pos(:,k), solvent_cells%edges, solvent)
          end do
          n_solvent = 0
          !$omp parallel do private(k,j) reduction(+:n_solvent)
          do k = 1, solvent%Nmax
             j = solvent%species(k)
             if (j <= 0) cycle
  