Mean-field molecular dynamics based on path integrals is used to approximate canonical quantum observables for particle systems consisting of nuclei and electrons. A computational bottleneck is the Monte Carlo sampling from the Gibbs density of the electron operator, which due to the fermion sign problem has a computational complexity that scales exponentially with the number of electrons. In this work, we construct an algorithm that approximates the mean-field Hamiltonian by path integrals for fermions. The algorithm is based on the determinant of a matrix with components built on Brownian bridges connecting permuted electron coordinates. The computational work for n electrons is O(n3), which reduces the computational complexity associated with the fermion sign problem. We analyze a bias resulting from this approximation and provide a rough computational error indicator. It remains to rigorously explain the surprisingly high accuracy for high temperatures. The method becomes infeasible at low temperatures due to a large sample variance.