I think we should implement a FreeAmplitude::fix() method so that
free_amplitude(*D, to(rho))->variableStatus() = VariableStatus::fixed;
*free_amplitude(*D, to(rho)) = std::complex<double>(1., 0.);
can be replaced by
free_amplitude(*D, to(rho))->fix(std::complex<double>(1., 0.));