anchor.diff:  
     call solvent% init(N,N_species)
   
  +  move = PTread_l(config, 'move')
     call colloids% init(1, N_species_colloids, mass)
     colloids% species(1) = 1
     colloids% vel = 0
  
positive.code:   double precision :: g
    ! trap parameters
    double precision :: k, trap_center(3)
  
    call PTparse(config,get_input_filename(),11)
  
  
negative.code:         e2 = compute_force_harmonic_trap(colloids, k, trap_center)
  
          call md_vel(solvent, dt)
          if (move) colloids% vel = colloids% vel + &
               dt * ( colloids% force + colloids% force_old ) / (2 * colloids% mass(1))
  
       end do md_loop
  