The expectation value of force is defined in Cardosa-Guitierrez, Remacle, J Phys B 2024. The implementation has some problems with the term coming from the diple operator
$F_\mu = \sum_n \sum_m \langle\chi_n| \frac{\partial \vec{\mu}_{nm}}{\partial R} | \chi_m \rangle \cdot \vec{E}(t)$
The term coming from PES was tested to be correct, yet this term appears problematic. The derivation is not fully clear to me. It should be checked and implemented correctly.
The full force can be calculated directly which allows for a comparison.
This issue is only one function in qdyn_analyze.py and does not affect QDyn on its own.