Numerical problems with small unit cells

[godotalgorithm]

This isn't really a bug, but more a fuzzy limitation to MOPAC's current functionality. When simulating periodic systems, MOPAC is primarily intended for large unit cells since it only uses the Gamma point of the Brillouin zone. Specifically, its implementation of periodic boundary conditions sometimes assumes that an atom is not coupled to its periodic image. Because MOPAC's models are based on Hartree-Fock calculations, the length scale of coupling through Fock exchange matrix elements depends on the decay of the density matrix. For large-gap insulators, this decay is rapid and supercells with lattice vectors of ~10 Angstroms should work fine. For metals and semiconductors, this decay is going to be a lot slower and it may not be possible to calculate a unit cell large enough to suppress this problem completely. It is not clear what precisely is afflicted by this problem, but energy derivatives and functionality that depends on them (e.g. geometry relaxation and vibrations) seem to be.

Although atoms being coupled to their periodic images is an unphysical finite-size effect, it is very helpful to be able to run reliable simulations in this regime to estimate the magnitude of finite-size effects and possibly extrapolate away their effects. It should be possible to repair MOPAC's functionality in this regime without major code revisions, but it will require careful testing and review of MOPAC's infrastructure for periodic systems.