<commit_before>program single_janus_pbc
  use rmpcdmd_module
  use hdf5
  use h5md_module
  use threefry_module
  use ParseText
  use iso_c_binding
  use omp_lib
  implicit none

  type(cell_system_t) :: solvent_cells
  type(particle_system_t) :: solvent
  type(particle_system_t) :: colloids
  type(neighbor_list_t) :: neigh
  type(lj_params_t) :: solvent_colloid_lj
  type(lj_params_t) :: colloid_lj

  integer, parameter :: N_species = 2

  integer :: rho
  integer :: N
  integer :: error

  double precision :: max_cut
  double precision :: epsilon(2,2)
  double precision :: sigma, sigma_v(2,2), sigma_cut(2,2)
  double precision :: mass(2)

  double precision :: e1, e2
  double precision :: tau, dt , T
  double precision :: prob
  double precision :: bulk_rate
  double precision :: skin, co_max, so_max
  integer :: N_MD_steps, N_loop
  integer :: n_extra_sorting

  type(PTo) :: config

  integer :: i, L(3)
  integer :: j, k

  double precision :: total_time
  type(timer_t), target :: varia, main, time_flag, time_refuel, time_change
  type(timer_list_t) :: timer_list
  integer(HID_T) :: timers_group

  type(threefry_rng_t), allocatable :: state(:)
  integer :: n_threads
  type(h5md_file_t) :: hfile
  type(h5md_element_t) :: dummy_element
  integer(HID_T) :: fields_group
  integer(HID_T) :: connectivity_group
  integer(HID_T) :: box_group
  integer(HID_T) :: tmp_id
  type(thermo_t) :: thermo_data
  double precision :: temperature, kin_e
  double precision :: v_com(3), x(3)
  double precision :: tmp_mass
  type(particle_system_io_t) :: janus_io
  type(particle_system_io_t) :: solvent_io

  integer, dimension(N_species) :: n_solvent
  type(h5md_element_t) :: n_solvent_el

  type(histogram_t) :: radial_hist
  type(h5md_element_t) :: radial_hist_el

  integer, allocatable :: links(:,:)
  double precision, allocatable :: links_d(:)
  double precision :: link_treshold
  integer :: i_link
  double precision :: dist, rattle_pos_tolerance, rattle_vel_tolerance

  type(args_t) :: args
  character(len=144) :: data_filename

  args = get_input_args()
  call PTparse(config, args%input_file, 11)

  n_threads = omp_get_max_threads()
  allocate(state(n_threads))
  call threefry_rng_init(state, args%seed)

  call main%init('main')
  call varia%init('varia')
  call time_flag%init('flag')
  call time_refuel%init('refuel')
  call time_change%init('change')

  call timer_list%init(16)
  call timer_list%append(varia)
  call timer_list%append(time_flag)
  call timer_list%append(time_refuel)
  call timer_list%append(time_change)

  call h5open_f(error)

  prob = PTread_d(config,'probability')
  bulk_rate = PTread_d(config,'bulk_rate')

  L = PTread_ivec(config, 'L', 3)
  rho = PTread_i(config, 'rho')
  N = rho *L(1)*L(2)*L(3)

  T = PTread_d(config, 'T')
  
  tau = PTread_d(config, 'tau')
  N_MD_steps = PTread_i(config, 'N_MD')
  dt = tau / N_MD_steps
  N_loop = PTread_i(config, 'N_loop')

  rattle_pos_tolerance = PTread_d(config, 'rattle_pos_tolerance')
  rattle_vel_tolerance = PTread_d(config, 'rattle_vel_tolerance')

  sigma = PTread_d(config, 'sigma_colloid')
  sigma_v = sigma
  sigma_cut = sigma_v*2**(1.d0/6.d0)
  mass = rho * sigma**3 * 4 * 3.14159265/3

  epsilon = PTread_d(config, 'epsilon_colloid')

  call colloid_lj% init(epsilon, sigma_v, sigma_cut)

  ! solvent index first, colloid index second, in solvent_colloid_lj
  sigma = PTread_d(config, 'sigma')
  sigma_v = sigma
  sigma_cut = sigma_v*2**(1.d0/6.d0)
  max_cut = maxval(sigma_cut)

  epsilon(:,1) = PTread_dvec(config, 'epsilon_C', 2)
  epsilon(:,2) = PTread_dvec(config, 'epsilon_N', 2)

  call solvent_colloid_lj% init(epsilon, sigma_v, sigma_cut)

  call solvent% init(N, N_species, system_name='solvent') !there will be 2 species of solvent particles

  call hfile%create(args%output_file, 'RMPCDMD::single_janus_pbc', &
       'N/A', '<NAME>')
  call thermo_data%init(hfile, n_buffer=50, step=N_MD_steps, time=N_MD_steps*dt)

  data_filename = PTread_s(config, 'data_filename')
  link_treshold = PTread_d(config, 'link_treshold')

  call PTkill(config)

  call colloids%init_from_file(data_filename, 'janus', H5MD_FIXED)
  colloids%image = 0
  colloids%vel = 0
  colloids%force = 0
  colloids%mass = mass

  i_link = 0
  do i = 1, colloids%Nmax
     do j = i+1, colloids%Nmax
        x = rel_pos(colloids%pos(:,i), colloids%pos(:,j), solvent_cells%edges)
        dist = norm2(x)
        if (dist < link_treshold) then
           ! count link
           i_link = i_link + 1
        end if
     end do
  end do
  allocate(links(2,i_link), links_d(i_link))

  i_link = 0
  do i = 1, colloids%Nmax
     do j = i+1, colloids%Nmax
        x = rel_pos(colloids%pos(:,i), colloids%pos(:,j), solvent_cells%edges)
        dist = norm2(x)
        if (dist < link_treshold) then
           ! add link
           i_link = i_link + 1
           links(1, i_link) = i
           links(2, i_link) = j
           links_d(i_link) = dist
        end if
     end do
  end do

  write(*,*) 'number of links:', i_link

  janus_io%force_info%store = .false.
  janus_io%id_info%store = .false.
  janus_io%position_info%store = .true.
  janus_io%position_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  janus_io%position_info%step = N_MD_steps
  janus_io%position_info%time = N_MD_steps*dt
  janus_io%image_info%store = .true.
  janus_io%image_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  janus_io%image_info%step = N_MD_steps
  janus_io%image_info%time = N_MD_steps*dt
  janus_io%velocity_info%store = .true.
  janus_io%velocity_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  janus_io%velocity_info%step = N_MD_steps
  janus_io%velocity_info%time = N_MD_steps*dt
  janus_io%species_info%store = .true.
  janus_io%species_info%mode = H5MD_FIXED
  call janus_io%init(hfile, 'janus', colloids)

  solvent_io%force_info%store = .false.
  solvent_io%id_info%store = .false.
  solvent_io%position_info%store = .true.
  solvent_io%position_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  solvent_io%position_info%step = N_loop*N_MD_steps
  solvent_io%position_info%time = N_loop*N_MD_steps*dt
  solvent_io%image_info%store = .true.
  solvent_io%image_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  solvent_io%image_info%step = N_loop*N_MD_steps
  solvent_io%image_info%time = N_loop*N_MD_steps*dt
  solvent_io%velocity_info%store = .true.
  solvent_io%velocity_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  solvent_io%velocity_info%step = N_loop*N_MD_steps
  solvent_io%velocity_info%time = N_loop*N_MD_steps*dt
  solvent_io%species_info%store = .true.
  solvent_io%species_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  solvent_io%species_info%step = N_loop*N_MD_steps
  solvent_io%species_info%time = N_loop*N_MD_steps*dt
  call solvent_io%init(hfile, 'solvent', solvent)

  do k = 1, solvent%Nmax
     solvent% vel(1,k) = threefry_normal(state(1))*sqrt(T)
     solvent% vel(2,k) = threefry_normal(state(1))*sqrt(T)
     solvent% vel(3,k) = threefry_normal(state(1))*sqrt(T)
  end do
  solvent% vel = solvent% vel - spread(sum(solvent% vel, dim=2)/solvent% Nmax, 2, solvent% Nmax)
  solvent% force = 0
  solvent% species = 1

  call solvent_cells%init(L, 1.d0)

  call n_solvent_el%create_time(hfile%observables, 'n_solvent', &
       n_solvent, H5MD_LINEAR, step=N_MD_steps, &
       time=N_MD_steps*dt)

  call h5gcreate_f(janus_io%group, 'box', box_group, error)
  call h5md_write_attribute(box_group, 'dimension', 3)
  call dummy_element%create_fixed(box_group, 'edges', solvent_cells%edges)
  call h5gclose_f(box_group, error)

  call h5gcreate_f(solvent_io%group, 'box', box_group, error)
  call h5md_write_attribute(box_group, 'dimension', 3)
  call dummy_element%create_fixed(box_group, 'edges', solvent_cells%edges)
  call h5gclose_f(box_group, error)

  call h5gcreate_f(hfile%id, 'fields', fields_group, error)
  call radial_hist%init(0.d0, 5*sigma, 100, 2)
  call radial_hist_el%create_time(fields_group, 'radial_histogram', radial_hist%data, &
       H5MD_LINEAR, step=N_MD_steps, time=N_MD_steps*dt)
  call h5md_write_attribute(radial_hist_el%id, 'xmin', radial_hist%xmin)
  call h5md_write_attribute(radial_hist_el%id, 'dx', radial_hist%dx)
  call h5gclose_f(fields_group, error)

  call h5gcreate_f(hfile%id, 'connectivity', connectivity_group, error)
  call h5md_write_dataset(connectivity_group, 'janus_links', links-1)
  call h5dopen_f(connectivity_group, 'janus_links', tmp_id, error)
  call h5md_write_attribute(tmp_id, 'particles_group', 'janus')
  call h5dclose_f(tmp_id, error)
  call h5gclose_f(connectivity_group, error)

  call solvent% random_placement(solvent_cells% edges, colloids, solvent_colloid_lj, state(1))

  call solvent% sort(solvent_cells)

  call neigh% init(colloids% Nmax, int(300*sigma**3))

  skin = 1.5
  n_extra_sorting = 0

  call neigh% make_stencil(solvent_cells, max_cut+skin)

  call neigh% update_list(colloids, solvent, max_cut+skin, solvent_cells, solvent_colloid_lj)

  ! insert polar histogram

  e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
  e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)
  solvent% force_old = solvent% force
  colloids% force_old = colloids% force

  write(*,*) 'Running for', N_loop, 'loops'
  call main%tic()
  do i = 1, N_loop
     if (modulo(i,5) == 0) write(*,'(i05)',advance='no') i
     md_loop: do j = 1, N_MD_steps
        call md_pos(solvent, dt)

        call varia%tic()
        ! Extra copy for rattle
        colloids% pos_rattle = colloids% pos
        do k=1, colloids% Nmax
           colloids% pos(:,k) = colloids% pos(:,k) + dt * colloids% vel(:,k) + &
                dt**2 * colloids% force(:,k) / (2 * colloids% mass(colloids%species(k)))
        end do
        call varia%tac()

        call rattle_body_pos(colloids, links, links_d, dt, solvent_cells% edges, rattle_pos_tolerance)

        so_max = solvent% maximum_displacement()
        co_max = colloids% maximum_displacement()

        if ( (co_max >= skin/2) .or. (so_max >= skin/2) ) then
           call varia%tic()
           call apply_pbc(solvent, solvent_cells% edges)
           call apply_pbc(colloids, solvent_cells% edges)
           call varia%tac()
           call solvent% sort(solvent_cells)
           call neigh% update_list(colloids, solvent, max_cut + skin, solvent_cells, solvent_colloid_lj)
           call varia%tic()
           solvent% pos_old = solvent% pos
           colloids% pos_old = colloids% pos
           n_extra_sorting = n_extra_sorting + 1
           call varia%tac()
        end if

        call varia%tic()
        call switch(solvent% force, solvent% force_old)
        call switch(colloids% force, colloids% force_old)<commit_msg>parallelize pos_old storage for solvent<commit_after>program single_janus_pbc
  use rmpcdmd_module
  use hdf5
  use h5md_module
  use threefry_module
  use ParseText
  use iso_c_binding
  use omp_lib
  implicit none

  type(cell_system_t) :: solvent_cells
  type(particle_system_t) :: solvent
  type(particle_system_t) :: colloids
  type(neighbor_list_t) :: neigh
  type(lj_params_t) :: solvent_colloid_lj
  type(lj_params_t) :: colloid_lj

  integer, parameter :: N_species = 2

  integer :: rho
  integer :: N
  integer :: error

  double precision :: max_cut
  double precision :: epsilon(2,2)
  double precision :: sigma, sigma_v(2,2), sigma_cut(2,2)
  double precision :: mass(2)

  double precision :: e1, e2
  double precision :: tau, dt , T
  double precision :: prob
  double precision :: bulk_rate
  double precision :: skin, co_max, so_max
  integer :: N_MD_steps, N_loop
  integer :: n_extra_sorting

  type(PTo) :: config

  integer :: i, L(3)
  integer :: j, k

  double precision :: total_time
  type(timer_t), target :: varia, main, time_flag, time_refuel, time_change
  type(timer_list_t) :: timer_list
  integer(HID_T) :: timers_group

  type(threefry_rng_t), allocatable :: state(:)
  integer :: n_threads
  type(h5md_file_t) :: hfile
  type(h5md_element_t) :: dummy_element
  integer(HID_T) :: fields_group
  integer(HID_T) :: connectivity_group
  integer(HID_T) :: box_group
  integer(HID_T) :: tmp_id
  type(thermo_t) :: thermo_data
  double precision :: temperature, kin_e
  double precision :: v_com(3), x(3)
  double precision :: tmp_mass
  type(particle_system_io_t) :: janus_io
  type(particle_system_io_t) :: solvent_io

  integer, dimension(N_species) :: n_solvent
  type(h5md_element_t) :: n_solvent_el

  type(histogram_t) :: radial_hist
  type(h5md_element_t) :: radial_hist_el

  integer, allocatable :: links(:,:)
  double precision, allocatable :: links_d(:)
  double precision :: link_treshold
  integer :: i_link
  double precision :: dist, rattle_pos_tolerance, rattle_vel_tolerance

  type(args_t) :: args
  character(len=144) :: data_filename

  args = get_input_args()
  call PTparse(config, args%input_file, 11)

  n_threads = omp_get_max_threads()
  allocate(state(n_threads))
  call threefry_rng_init(state, args%seed)

  call main%init('main')
  call varia%init('varia')
  call time_flag%init('flag')
  call time_refuel%init('refuel')
  call time_change%init('change')

  call timer_list%init(16)
  call timer_list%append(varia)
  call timer_list%append(time_flag)
  call timer_list%append(time_refuel)
  call timer_list%append(time_change)

  call h5open_f(error)

  prob = PTread_d(config,'probability')
  bulk_rate = PTread_d(config,'bulk_rate')

  L = PTread_ivec(config, 'L', 3)
  rho = PTread_i(config, 'rho')
  N = rho *L(1)*L(2)*L(3)

  T = PTread_d(config, 'T')
  
  tau = PTread_d(config, 'tau')
  N_MD_steps = PTread_i(config, 'N_MD')
  dt = tau / N_MD_steps
  N_loop = PTread_i(config, 'N_loop')

  rattle_pos_tolerance = PTread_d(config, 'rattle_pos_tolerance')
  rattle_vel_tolerance = PTread_d(config, 'rattle_vel_tolerance')

  sigma = PTread_d(config, 'sigma_colloid')
  sigma_v = sigma
  sigma_cut = sigma_v*2**(1.d0/6.d0)
  mass = rho * sigma**3 * 4 * 3.14159265/3

  epsilon = PTread_d(config, 'epsilon_colloid')

  call colloid_lj% init(epsilon, sigma_v, sigma_cut)

  ! solvent index first, colloid index second, in solvent_colloid_lj
  sigma = PTread_d(config, 'sigma')
  sigma_v = sigma
  sigma_cut = sigma_v*2**(1.d0/6.d0)
  max_cut = maxval(sigma_cut)

  epsilon(:,1) = PTread_dvec(config, 'epsilon_C', 2)
  epsilon(:,2) = PTread_dvec(config, 'epsilon_N', 2)

  call solvent_colloid_lj% init(epsilon, sigma_v, sigma_cut)

  call solvent% init(N, N_species, system_name='solvent') !there will be 2 species of solvent particles

  call hfile%create(args%output_file, 'RMPCDMD::single_janus_pbc', &
       'N/A', '<NAME>')
  call thermo_data%init(hfile, n_buffer=50, step=N_MD_steps, time=N_MD_steps*dt)

  data_filename = PTread_s(config, 'data_filename')
  link_treshold = PTread_d(config, 'link_treshold')

  call PTkill(config)

  call colloids%init_from_file(data_filename, 'janus', H5MD_FIXED)
  colloids%image = 0
  colloids%vel = 0
  colloids%force = 0
  colloids%mass = mass

  i_link = 0
  do i = 1, colloids%Nmax
     do j = i+1, colloids%Nmax
        x = rel_pos(colloids%pos(:,i), colloids%pos(:,j), solvent_cells%edges)
        dist = norm2(x)
        if (dist < link_treshold) then
           ! count link
           i_link = i_link + 1
        end if
     end do
  end do
  allocate(links(2,i_link), links_d(i_link))

  i_link = 0
  do i = 1, colloids%Nmax
     do j = i+1, colloids%Nmax
        x = rel_pos(colloids%pos(:,i), colloids%pos(:,j), solvent_cells%edges)
        dist = norm2(x)
        if (dist < link_treshold) then
           ! add link
           i_link = i_link + 1
           links(1, i_link) = i
           links(2, i_link) = j
           links_d(i_link) = dist
        end if
     end do
  end do

  write(*,*) 'number of links:', i_link

  janus_io%force_info%store = .false.
  janus_io%id_info%store = .false.
  janus_io%position_info%store = .true.
  janus_io%position_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  janus_io%position_info%step = N_MD_steps
  janus_io%position_info%time = N_MD_steps*dt
  janus_io%image_info%store = .true.
  janus_io%image_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  janus_io%image_info%step = N_MD_steps
  janus_io%image_info%time = N_MD_steps*dt
  janus_io%velocity_info%store = .true.
  janus_io%velocity_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  janus_io%velocity_info%step = N_MD_steps
  janus_io%velocity_info%time = N_MD_steps*dt
  janus_io%species_info%store = .true.
  janus_io%species_info%mode = H5MD_FIXED
  call janus_io%init(hfile, 'janus', colloids)

  solvent_io%force_info%store = .false.
  solvent_io%id_info%store = .false.
  solvent_io%position_info%store = .true.
  solvent_io%position_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  solvent_io%position_info%step = N_loop*N_MD_steps
  solvent_io%position_info%time = N_loop*N_MD_steps*dt
  solvent_io%image_info%store = .true.
  solvent_io%image_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  solvent_io%image_info%step = N_loop*N_MD_steps
  solvent_io%image_info%time = N_loop*N_MD_steps*dt
  solvent_io%velocity_info%store = .true.
  solvent_io%velocity_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  solvent_io%velocity_info%step = N_loop*N_MD_steps
  solvent_io%velocity_info%time = N_loop*N_MD_steps*dt
  solvent_io%species_info%store = .true.
  solvent_io%species_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  solvent_io%species_info%step = N_loop*N_MD_steps
  solvent_io%species_info%time = N_loop*N_MD_steps*dt
  call solvent_io%init(hfile, 'solvent', solvent)

  do k = 1, solvent%Nmax
     solvent% vel(1,k) = threefry_normal(state(1))*sqrt(T)
     solvent% vel(2,k) = threefry_normal(state(1))*sqrt(T)
     solvent% vel(3,k) = threefry_normal(state(1))*sqrt(T)
  end do
  solvent% vel = solvent% vel - spread(sum(solvent% vel, dim=2)/solvent% Nmax, 2, solvent% Nmax)
  solvent% force = 0
  solvent% species = 1

  call solvent_cells%init(L, 1.d0)

  call n_solvent_el%create_time(hfile%observables, 'n_solvent', &
       n_solvent, H5MD_LINEAR, step=N_MD_steps, &
       time=N_MD_steps*dt)

  call h5gcreate_f(janus_io%group, 'box', box_group, error)
  call h5md_write_attribute(box_group, 'dimension', 3)
  call dummy_element%create_fixed(box_group, 'edges', solvent_cells%edges)
  call h5gclose_f(box_group, error)

  call h5gcreate_f(solvent_io%group, 'box', box_group, error)
  call h5md_write_attribute(box_group, 'dimension', 3)
  call dummy_element%create_fixed(box_group, 'edges', solvent_cells%edges)
  call h5gclose_f(box_group, error)

  call h5gcreate_f(hfile%id, 'fields', fields_group, error)
  call radial_hist%init(0.d0, 5*sigma, 100, 2)
  call radial_hist_el%create_time(fields_group, 'radial_histogram', radial_hist%data, &
       H5MD_LINEAR, step=N_MD_steps, time=N_MD_steps*dt)
  call h5md_write_attribute(radial_hist_el%id, 'xmin', radial_hist%xmin)
  call h5md_write_attribute(radial_hist_el%id, 'dx', radial_hist%dx)
  call h5gclose_f(fields_group, error)

  call h5gcreate_f(hfile%id, 'connectivity', connectivity_group, error)
  call h5md_write_dataset(connectivity_group, 'janus_links', links-1)
  call h5dopen_f(connectivity_group, 'janus_links', tmp_id, error)
  call h5md_write_attribute(tmp_id, 'particles_group', 'janus')
  call h5dclose_f(tmp_id, error)
  call h5gclose_f(connectivity_group, error)

  call solvent% random_placement(solvent_cells% edges, colloids, solvent_colloid_lj, state(1))

  call solvent% sort(solvent_cells)

  call neigh% init(colloids% Nmax, int(300*sigma**3))

  skin = 1.5
  n_extra_sorting = 0

  call neigh% make_stencil(solvent_cells, max_cut+skin)

  call neigh% update_list(colloids, solvent, max_cut+skin, solvent_cells, solvent_colloid_lj)

  ! insert polar histogram

  e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
  e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)
  solvent% force_old = solvent% force
  colloids% force_old = colloids% force

  write(*,*) 'Running for', N_loop, 'loops'
  call main%tic()
  do i = 1, N_loop
     if (modulo(i,5) == 0) write(*,'(i05)',advance='no') i
     md_loop: do j = 1, N_MD_steps
        call md_pos(solvent, dt)

        call varia%tic()
        ! Extra copy for rattle
        colloids% pos_rattle = colloids% pos
        do k=1, colloids% Nmax
           colloids% pos(:,k) = colloids% pos(:,k) + dt * colloids% vel(:,k) + &
                dt**2 * colloids% force(:,k) / (2 * colloids% mass(colloids%species(k)))
        end do
        call varia%tac()

        call rattle_body_pos(colloids, links, links_d, dt, solvent_cells% edges, rattle_pos_tolerance)

        so_max = solvent% maximum_displacement()
        co_max = colloids% maximum_displacement()

        if ( (co_max >= skin/2) .or. (so_max >= skin/2) ) then
           call varia%tic()
           call apply_pbc(solvent, solvent_cells% edges)
           call apply_pbc(colloids, solvent_cells% edges)
           call varia%tac()
           call solvent% sort(solvent_cells)
           call neigh% update_list(colloids, solvent, max_cut + skin, solvent_cells, solvent_colloid_lj)
           call varia%tic()
           !$omp parallel do
           do k = 1, solvent%Nmax
              solvent% pos_old(:,k) = solvent% pos(:,k)
           end do
           colloids% pos_old = colloids% pos
           n_extra_sorting = n_extra_sorting + 1
           call varia%tac()
        end if

        call varia%tic()
        call switch(solvent% force, solvent% force_old)
        call switch(colloids% force, colloids% force_old)
{'fraction_non_alphanumeric': 0.1142337807606264, 'fraction_numerical': 0.014587994034302759, 'mean_word_length': 4.392560944961046, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4985455', 'n_tokens_mistral': 8797, 'n_tokens_neox': 8220, 'n_words': 1586}