Comments on implementation

[zerothi]

I have some general remarks on the implementation of your SCF model

I know it is a development code, so please see these as comments and not as criticism ;)

I'll refer to HH as the HubbardHamiltonian.

• Generally I think having separate IO and models would be better, i.e. you should spilt HH into HH and IO. The IO class should then inherit from the SileCDF class which enables all methods to easily create variables, groups etc.

  See e.g. SileCDF._crt_*

   grp = SileCDF._crt_group(self, 'GROUP')
   grp2 = SileCDF._crt_group(grp, 'GROUP') # inside above group
  

  If the group already exists, the handle for the group is returned.

  The same functionality applies to _crt_dim/_crt_var, see for examples on how to use these sisl/io/siesta/siesta_nc.py

  Once you have a class that enables writing a <>.HU.nc file (HubbardU-file) then you simply do:

  get_sile('name.HU.nc', 'w').write(HubbardHamiltonian) # however see further down
  

• From my opinion it seems that the HH is very similar to the regular Hamiltonian, in this regard it may be more ideal to only implement a class which does Hubbard U modelling using an external Hamiltonian.

Thus I would make a new class, HubbardSCF which implements the required functionality

  # U could be a single number (same U for all orbitals)
  #   or an array of numbers (different U for different orbitals)
  hscf = HubbardSCF(Hamiltonian, U, tol=<SCF-tolerance>)
  hscf.SCF(...) # instead of converge

Instead of hm.iterate you could perhaps look into using Python's build-in iterator

  hscf.__iter__(self)

this would let you do:

  for i, dq in enumerate(hscf):
  	  # dq is the current change in charge
      print('SCF step {} dq={}'.format(i, dq))

Then the HubbardSCF can be extended more generically which doesn't put restraints on the Hamiltonian.

Doing this I think is more generic because the Hubbard-U is not specific to the Hamiltonian, but rather to the way the resulting Hamiltonian looks like. It is the SCF that determines H.

To conclude I think you shouldn't implement a separate Hamiltonian, but rather a separate HubbardSCF which does everything, hold U, define the SCF etc.

From my experience, the more general one class becomes, the easier it is to adapt to different environments. And it will be possible to much more easily extend the code with new functionality.