In the Huckel method, you are just generating a very simplified version of the molecular Hamiltonian and determining it's eigenvalues. The molecular Hamiltonian will always be a Hermitian matrix (for the Huckel method, we can be more specific and say it's a real symmetric matrix by construction). Hermitian matrices are guaranteed to have all real eigenvalues ...
ASE has a FixSymmetry constraint that preserves spacegroup symmetries. It works with a variety of structure optimization algorithms. You could use LAMMPS as the engine or one of the many other calculators, including some DFT options.
Based on your input file, I tried to run the calculation in my machine. As you didn't specified how did you obtain the original structure (assembled by hand or downloaded from structure database), I first generated a SDF file of your structure with Avogadro, and did a conformer search with obconformer from Open Babel, just to be sure to pick a initial ...
If you look at your last five iterations, you will see that the density change (which is what matters most for SCF convergence) is getting smaller in each iteration, meaning that the calculation is still converging, and will converge if you give it more iterations:
GULP (General Utility Lattice Program) uses specific forcefields for different type of systems and it is designed to work with periodic conditions. Also, it has defined several type of potential models for two-body (Buckingham, Lennard-Jones, Morse, etc.), three-body (Three-body harmonic, Axilrod-Teller, Stillinger-Weber, etc.), four-body, six-body and ...