Any combination of methods for geometry optimization and wavefunction optimization is allowed. Possible options for geometry optimization are GDIIS, RFO, BFGS, and steepest descent. If you choose steepest descent for both, geometry variables and the wavefunction, a combined method is used. For all other combinations a full wavefunction optimization is performed between changes of the ionic coordinates. The default options are GDIIS and ODIIS. The quasi-Newton methods (GDIIS, RFO, and BFGS) are using the BFGS method to update an approximate Hessian. At the beginning of a run the Hessian can either be initialized as a unit matrix HESSIAN UNIT or with an empirical force field. Two force fields are implemented: The DISCO and the SCHLEGEL force field. The algorithm for the empirical force fileds has to identify bonds in the system. For unusual geometries this may fail and the Hessian becomes singular. To prevent this you can add or delete bonds with the keyword CHANGE BONDS.