A similar expression can be obtained when using the adiabatic-Fock basis. The nuclear motion can be solved using a velocity-verlet algorithm. The nuclear DOFs can be initialized by sampling its thermal distribution on the ground state potential energy surface around the Frank-Condon region at a given temperature either by use of BO molecular dynamics (BOMD) using randomly sampled positions and velocities over long-time dynamics or via sampling the classical Wigner distribution. Both methods can be performed at arbitrary temperatures up to the point where the normal mode analysis breaks down, at which point the system needs to be sampled via BOMD to obtain a meaningful distribution in a highly anharmonic ground state potential.