Nonlinear processes in hot, magnetized plasma are notoriously difficult to understand without the use of numerical simulations. In recent decades, first principles, kinetic simulations have been widely and successfully used to study plasma turbulence and reconnection in weakly collisional systems. In this thesis, extensions of well-known, Lagrangian, particle-in-cell (PIC) simulation algorithms for problems such as these are derived and implemented. The algorithms are tested for multiple species (electrons and ions, with the physical mass ratio) in non-trivial magnetic geometry (cylindrical/toroidal). The advances presented here address two major shortcomings of conventional gyrokinetic PIC algorithms, with demonstrated excellent performance on large, parallel supercomputers. Although the gyrokinetic formalism rigorously describes the evolution of fluctuations which are small compared to a typical Larmor radius, most existing algorithms use low-order approximations of the gyroaveraging operator, and cannot accurately describe small scale fluctuations. The gyroaveraging algorithm presented here accurately and uniquely treats a wide range of fluctuation scales, above and below the thermal gyroradius. The second shortcoming of traditional algorithms relates to the slow loss of accuracy that is associated with the build-up of noise. In this thesis, a PIC pitch-angle scattering collision operator is developed. This collision operator is physically motivated and controls the growth of noise without introducing non-physical dissipation. Basic tests of the new algorithms are presented in linear and nonlinear regimes, using one to thousands of processors simultaneously.