  We propose to relax geometries throughout chemical compound space (CCS) using
alchemical perturbation density functional theory (APDFT). APDFT refers to
perturbation theory involving changes in nuclear charges within approximate
solutions to Schr\"odinger's equation. We give an analytical formula to
calculate the mixed second order energy derivatives with respect to both,
nuclear charges and nuclear positions (named "alchemical force"), within the
restricted Hartree-Fock case. We have implemented and studied the formula for
its use in geometry relaxation of various reference and target molecules. We
have also analysed the convergence of the alchemical force perturbation series,
as well as basis set effects. Interpolating alchemically predicted energies,
forces, and Hessian to a Morse potential yields more accurate geometries and
equilibrium energies than when performing a standard Newton Raphson step. Our
numerical predictions for small molecules including BF, CO, N2, CH$_4$, NH$_3$,
H$_2$O, and HF yield mean absolute errors of of equilibrium energies and bond
lengths smaller than 10 mHa and 0.01 Bohr for 4$^\text{th}$ order APDFT
predictions, respectively. Our alchemical geometry relaxation still preserves
the combinatorial efficiency of APDFT: Based on a single coupled perturbed
Hartree Fock derivative for benzene we provide numerical predictions of
equilibrium energies and relaxed structures of all the 17 iso-electronic
charge-netural BN-doped mutants with averaged absolute deviations of $\sim$27
mHa and $\sim$0.12 Bohr, respectively.
