We develop generalization of the fixed-phase diffusion Monte Carlo method for Hamiltonians which explicitly depend on particle spins such as for spin-orbit interactions. The method is formulated in zero variance manner and is similar to treatment of nonlocal operators in commonly used staticspin calculations. Tests on atomic and molecular systems show that it is very accurate, on par with the fixed-node method. This opens electronic structure quantum Monte Carlo methods to a vast research area of quantum phenomena in which spin-related interactions play an important role.PACS numbers: 02.70. Ss, 31.30.jc Quantum Monte Carlo (QMC) methods are making significant contributions to our understanding of manybody effects in quantum systems. Although hampered by the infamous fermion sign problem, a number of approaches have been explored for dealing with inefficiencies whenever sampled distributions possess varying signs or complex values. One of the commonly used strategies is the fixed-node approximation that replaces the fermionic antisymmetry with boundaries given by trial wave function nodes. For broken time-reversal Hamiltonians or for twisted boundary conditions [1] with inherently complex eigenstates, the fixed-node condition has been generalized to the fixed-phase approximation [2]. Benchmark quality results for both models and real materials have been obtained in many settings such as molecules, solids, non-covalently bonded complexes, ultracold condensates and other systems [3,4].Electronic structure QMC calculations are usually done with particle spins being assigned fixed labels, up or down. Since spins commute with Hamiltonians without explicit spin terms, the problem simplifies to the spatial solution of the stationary Schrödinger equation. Treating the spins as quantum variables for more complicated Hamiltonians was explored very early [5] in variational Monte Carlo (VMC) of nuclear systems. However, extending this to projection methods such as the diffusion Monte Carlo (DMC) in position space [3,4] is much more involved. Building upon results for nuclear systems [6][7][8], a DMC method has been proposed and applied to a 2D electron gas with Rashba spin-orbit interaction [9]. In this approach the spinors are stochastically updated by the action of the spin-orbit operator that is absorbed into the path sampling part of the propagator. It effectively samples the space of spinor states rather than (spin) coordinates and a similar VMC approach has been implemented for spin-orbit in atoms [10] very recently.Here we propose a new development that is formulated as the DMC method in coordinate space with spinors in the trial state kept intact during the imaginary time evolution. This implies the zero variance property, ie, the bias in the obtained energy is proportional to the square of the trial function error.The method builds upon our previous work [11] on nonlocal pseudopotentials (PP) since the spin-orbit operator is just another case of inherent nonlocality. It is also well suited for calculations of re...