anchor.diff:      l = nc(i)
       do j = 1, l
         nb = nbs(i,j)
  -
         dx = px(i) - px(nb)
         dy = py(i) - py(nb)
         dz = pz(i) - pz(nb)
  
positive.code:   close(unit)
    deallocate(data)
  
    call mass_change
  
    ! Print out initial particle positions
    frame = 0
    call particle_write
  
negative.code:   ps_t = 0.0_WP
    tot_t = 0.0_WP
  
    call mass_change
  
    return
  end subroutine init
  
  