anchor.diff:            n_effective = n_effective + 1
          end if
       end do
  +
  +    write(31,*) 'time', current_time, 'n_effective', n_effective, 'energy', energy
  +    if (n_effective == 0) then
  +       stop 'n_effective == 0 in add_energy_to_cell'
  +    end if
  +    if (n_effective == 1) then
  +       stop 'n_effective == 1 in add_energy_to_cell'
  +    end if
       com_v = com_v / n_effective
   
       kin = 0
  
positive.code:   subroutine bind_molecule(idx)
      integer, intent(in) :: idx
      double precision :: com_v(3), w_ab(3), excess_kinetic_energy, mu_ab
      integer :: p(3), cell_idx
  
      ! adjust velocity
  
negative.code:   double precision :: mass(2)
  
    double precision :: elastic_k, link_d(2)
    double precision :: angle_k, link_angle, link_angles(2)
  
    double precision :: e1, e2, e3
    double precision :: tau, dt , T
  