anchor.diff:    double precision :: g
     ! trap parameters
     double precision :: k, trap_center(3)
  +  logical :: move
   
     call PTparse(config,get_input_filename(),11)
   
  
positive.code:         e2 = compute_force_harmonic_trap(colloids, k, trap_center)
  
          call md_vel(solvent, dt)
          colloids% vel = colloids% vel + &
               dt * ( colloids% force + colloids% force_old ) / (2 * colloids% mass(1))
  
       end do md_loop
  
negative.code:      md_loop: do j = 1, N_MD_steps
          call mpcd_stream_nogravity_zwall(solvent, solvent_cells, dt)
  
  
          if (move) colloids% pos = colloids% pos + dt * colloids% vel + &
               dt**2 * colloids% force / (2*colloids% mass(1))
  
          so_max = solvent% maximum_displacement()
  