anchor.diff:    double precision :: g
     ! trap parameters
     double precision :: k, trap_center(3)
  +  logical :: move
   
     call PTparse(config,get_input_filename(),11)
   
  
positive.code: 
    call solvent% init(N,N_species)
  
    call colloids% init(1, N_species_colloids, mass)
    colloids% species(1) = 1
    colloids% vel = 0
  
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
  