The variational fitting of the Fock potential employing localized molecular orbitals requires either the inversion of the local two-center Coulomb matrices or alternatively the solution of corresponding linear equation systems with these matrices. In both cases, the method of choice is the Cholesky decomposition of the formally positive definite local two-center Coulomb matrices. However, due to finite-precision round-off errors, the local Coulomb matrices may be indefinite, and thus, the Cholesky decomposition is not applicable. To overcome this problem, we propose to make use of a modified Cholesky decomposition based on the indefinite factorization of local two-center Coulomb matrices. To this end, the working equations for the use of the modified Cholesky decomposition within the variational fitting of the Fock potential are presented. Benchmark calculations with global and range-separated hybrid functionals show that the proposed method can improve considerably the workload balance in parallel calculations.