This work presents a variationally fitted long-range exact exchange algorithm that can be used for the computation of range-separated hybrid density functionals in the linear combination of Gaussian type orbital (LCGTO) approximation. The obtained LCGTO energy and gradient expressions are free of four-center integrals and employ modified three-center integral recurrence relations to obtain optimal computational performance. The accuracy and performance of selected range-separated hybrid functionals with variational fitted long-range exact exchange are analyzed and discussed. A parallel implementation is presented and benchmarked.