We present a parallel implementation of a novel method to compute the NMR chemical shift in condensed matter. We designed the algorithm in a parallel fashion from the very beginning, and could reuse building blocks from our parallel eigensolver. By using a linear systems solver to avoid computing the full spectrum of a large matrix, our scheme becomes suitable for large systems. For the conjugate-gradient linear solver, we suggest a blocking strategy for better cache utilization. The most difficult part of the parallel implementation is a 3d-FFT, which we describe in detail. Fortran 90 and MPI are used to arrive at an efficient, portable and easy-to-use code. We present performance numbers for the computational kernels on our production platforms: Cray T3E, SGI PowerChallenge, and IBM SP2. The good return times, processing power, and memory available on parallel machines has allowed us to compute within a short time frame the chemical shifts of several interesting systems.