anchor.diff:  
     call colloid_lj% init(epsilon, sigma_v, sigma_cut)
   
  +  ! no wall in x and z
  +  wall_sigma = -1
  +  wall_shift = 0
  +
  +  do_ywall = PTread_l(config, 'do_ywall', loc=params_group)
  +  if (do_ywall) then
  +     wall_sigma(2,:) = PTread_d(config, 'wall_sigma', loc=params_group)
  +     wall_shift(2) = PTread_d(config, 'wall_shift', loc=params_group)
  +     wall_epsilon = PTread_d(config, 'wall_epsilon', loc=params_group)
  +     call walls_colloid_lj% init(wall_epsilon, &
  +          wall_sigma, 3.d0**(1.d0/6.d0)*wall_sigma, wall_shift)
  +  end if
  +
     ! solvent index first, colloid index second, in solvent_colloid_lj
     sigma = PTread_d(config, 'sigma', loc=params_group)
     sigma_v = sigma
  
positive.code:   call axial_cf%init(block_length, N_loop, N_loop*N_MD_steps)
  
    call colloids%init_from_file(data_filename, data_group, H5MD_FIXED)
  
    call input_data_file%open(data_filename, H5F_ACC_RDONLY_F)
  
  
negative.code:         e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
          if (do_lennard_jones) e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)
          if (do_elastic) e3 = elastic_network(colloids, links, links_d, elastic_k, solvent_cells%edges)
          if (do_ywall) e_wall = lj93_zwall(colloids, solvent_cells% edges, walls_colloid_lj)
  
          call cell_md_vel(solvent_cells, solvent, dt, md_flag=.true.)
  
  