  Efficient representations of the Hamiltonian such as double factorization
drastically reduce circuit depth or number of repetitions in error corrected
and noisy intermediate scale quantum (NISQ) algorithms for chemistry. We report
a Lagrangian-based approach for evaluating relaxed one- and two-particle
reduced density matrices from double factorized Hamiltonians, unlocking
efficiency improvements in computing the nuclear gradient and related
derivative properties. We demonstrate the accuracy and feasibility of our
Lagrangian-based approach to recover all off-diagonal density matrix elements
in classically-simulated examples with up to 327 quantum and 18470 total atoms
in QM/MM simulations, with modest-sized quantum active spaces. We show this in
the context of the variational quantum eigensolver (VQE) in case studies such
as transition state optimization, ab initio molecular dynamics simulation and
energy minimization of large molecular systems.
