An anisotropic auxiliary differential equation finite-difference time-domain formulation is presented in detail for the time-domain study of nematic liquid crystal devices in the terahertz spectrum. The termination of the computation domain is achieved by employing a properly designed convolution perfectly matched layer. The material dispersion and dichroism of the LC complex permittivities is modeled via a modified Lorentzian function that is demonstrated to provide an accurate description for a series of state-of-the-art materials used in LC-THz technology.