  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 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 based on Brownian bridges
connecting permuted electron coordinates. The computational work for $n$
electrons is $\mathcal O(n^3)$, which reduces the computational complexity
associated with the fermion sign problem. We analyze a bias resulting from this
approximation and provide a computational error indicator. It remains to
rigorously explain the surprisingly high accuracy.
