It probably won't happen, but it's a possible enhancement.

It would require a lot of changes, the computing power requirement will increase a lot...

Instead of a single variational parameter, a lot of them would exist. Better optimization methods could be used, for some of them, look into the DFT Quantum Dot project.

The current stats are pretty bad, because of autocorrelation. Need to be improved, using binning.

As it is, the code won't find the optimum beta for the more complex cases and also the error estimation is very bad.

As the title says. I committed the code before having it finished.

A lot of steps are currently required for 'thermalization' since the initial guess is very bad. Use Thomas-Fermi to improve that https://en.wikipedia.org/wiki/Thomas%E2%80%93Fermi_model

Right now the valence orbitals that are combined for the molecule are not chosen very well. In fact, with certain combinations of atoms the program will crash because not enough orbitals are combined. This should be improved.

This could probably be done more generally, not only for diatomic molecules.

I added some comments in the code about it. Basically some code from the Hartree-Fock could be reused to compute the overlap matrix for the atomic contracted gaussian orbitals from the valence shells of the atoms. The overlap value could be used to select the pair of orbitals to use to generate norm * (O1 + O2) and norm * (O1 - O2) 'molecular' orbitals to be used in the Slater determinant.