dipc-cc / hubbard Goto Github PK

Revise code, in particular with respect to the print statement.

Currently, over the course of iterations, we only mix mix in a fraction of the instantaneous solution for the electron density. It would be nice to have some more advanced schemes, such as Pulay/DIIS mixing, for systems that are more tricky to converge.

https://en.wikipedia.org/wiki/DIIS

Restructure code to separate the core functionality and post-processing tools (plotting etc)

Currently only 1D periodicity along the x-axis is correctly handled. This should be generalized to arbitrary directions and dimensionalities.

The random_density could be improved slightly

Here are some general thoughts that may be helpful for you when reconsidering what to do:

• It is a little annoying that the random_density prints to std-out. I would simply remove it, generally remove as many print-statements and provide a logger functionality if required (I plan on doing something with the mixer, but I haven't looked at it yet).
• It would be really nice if one could (optionally) add a random_density instead of overwriting it, imagine a converged system, just to give it a bump one could add a random density (with some weight), for instance HH.dm += np.random.rand((2, no)) * 0.1 - 0.05. The thing is that sometimes one is not fully sure the SCF convergence finds the global minimum, for some of your systems I found that switching between Pulay and Linear mixing could result in two different minima. I didn't check their configurations.

As mentioned in 1359760 I find that the Hamiltonian is set up incorrectly (hopping terms in the second spin index not written) if one removes the explicit spin indexing with :. I would like to understand this and find a solution.

A simple test of sislis not revealing any such issues:

import sisl
import numpy as np

b = 10
a = 1000
chain = np.zeros((a,3))
for i in range(a):
    chain[i, 0] += i

g = sisl.Geometry(chain)
H = sisl.Hamiltonian(g, dim=2)

for ia in g:
    idx = g.close(ia, R=[0.1, 1.1])
    H[ia, idx[0]] = ia
    H[ia, idx[1]] = -ia

for ispin in range(2):
    print H.Hk(format='dense', spin=ispin)[:b, :b]

In the above, both spin channels are correctly initialized (and identical in this case).

When the user reads a geometry (.XV or .xyz formats) it seems natural to determine at this stage center and extension for the various plotting functionalities. Perhaps the user wants origo to be the center-of-mass or some specific atom.

Automatic documentation with Sphinx was initiated in de99c53. Try:

source install.sh
cd doc
source run.sh
firefox latest/index.html

However, for this to be useful we need a lot of in-code documentation strings.

Hi, I was wondering if there is a way as we do with .set_polarization() for maximizing up and down polarization to impose a zero polarization on a specific atomic site.

Thanks in advance!

Hi, I'm trying to plot the DOS using the plot.DOS for different values of U in the same plot using a loop. But the usual plt.plot and plt.show of matplotlib doesn't seem to work.

I can see in your code you are using:

def update_hamiltonian(...)

    atom.Z

to calculate the initial charge.
However, I think you will be better of by letting the Atom carry the charge (in that way you have a much better control).

I have a small test implementation:

1. Code change

diff --git a/Hubbard/hamiltonian.py b/Hubbard/hamiltonian.py
index 057c484..b3a8cb8 100644
--- a/Hubbard/hamiltonian.py
+++ b/Hubbard/hamiltonian.py
@@ -107,21 +107,12 @@ class HubbardHamiltonian(object):
 
     def update_hamiltonian(self):
         # Update spin Hamiltonian
-        g = self.geom
-
-        # Faster to loop individual species
-        E = np.empty([len(g), 2])
-        E[:, 0] = self.U * self.ndn
-        E[:, 1] = self.U * self.nup
-        for atom, ias in g.atoms.iter(True):
-            # charge on neutral atom
-            n0 = atom.Z - 5
-
-            # Faster to do it for more than one element at a time
-            E[ias, :] -= self.U * n0
-
-            for ia in ias:
-                self.H[ia, ia, [0, 1]] = E[ia] + self.e0[ia]
+        q0 = self.geom.atoms.q0
+        E = self.e0.copy()
+        E[:, 0] += self.U * (self.ndn - q0)
+        E[:, 1] += self.U * (self.nup - q0)
+        a = np.arange(len(self.H))
+        self.H[a, a, [0, 1]] = E

And then your kondo-paper example:

diff --git a/examples/molecules/kondo-paper/fig_S11-S13.py b/examples/molecules/kondo-paper/fig_S11-S13.py
index 2fc0d77..9d996ee 100644
--- a/examples/molecules/kondo-paper/fig_S11-S13.py
+++ b/examples/molecules/kondo-paper/fig_S11-S13.py
@@ -8,6 +8,10 @@ import sisl
 
 # Build sisl Geometry object
 mol = sisl.get_sile('junction-2-2.XV').read_geometry()
+for atom in mol.atoms.iter():
+    # Initialize charges
+    if atom.Z in [5, 6, 7]:
+        atom.orbital[0].q0 = atom.Z - 5
 mol.sc.set_nsc([1,1,1])
 mol = mol.move(-mol.center(what='xyz')).rotate(220, [0,0,1])
 Hsp2 = sp2(mol, t1=2.7, t2=0.2, t3=.18)

In this way you can turn on-and off selected sites without having to rely on atomic numbers.

Well, just a thought. :)
Probably it should also be somewhat faster.
This will let your geometry carry information about the charge and thus you can now do:
geometry.q0 to get an array of initial charges per atom.

This will allow you to make your update_hamiltonian faster.

We need a general method to split a given geometry into (perhaps just approximate) sublattices.

At the moment the classification is based on the assumption that alternating atom indices is about the right division, but this may not generally be the atom ordering supplied by the user.

This would boil down to forcing nup=ndn. In other words, it would be sufficient to work with a Hamiltonian of dim=1 and iterate over one spin component only.

One could have the code to choose this method if the initial single-particle Hamiltonian is only of dim=1 or if the user explicitly sets spin_degenerate=True.

Using sisl.MonkhorstPack should be enough

With test-hubbard-open.py, and many examples, I get AttributeError: module 'sisl' has no attribute 'WideBandSE'. Looking at the documentation of sisl, the correct class to use is sisl.physics.WideBandSE instead of sisl.WideBandSE.
My packages:
sisl 0.10.0 py38h1b0416a_0 conda-forge
python 3.8.8 hdb3f193_4
scipy 1.6.2 py38h91f5cce_0
numpy 1.19.2 py38h54aff64_0

Full error message:

Traceback (most recent call last):

File "/home/leo/Documents/hubbard/tests/test-hubbard-open.py", line 39, in
negf = NEGF(MFH_HC, [(MFH_elec, '-A'), (MFH_elec, '+A')], elec_indx)

File "/opt/anaconda3/lib/python3.8/site-packages/hubbard/negf.py", line 158, in init
if all(map(lambda obj: isinstance(obj, sisl.WideBandSE), self.elec_SE)):

File "/opt/anaconda3/lib/python3.8/site-packages/hubbard/negf.py", line 158, in
if all(map(lambda obj: isinstance(obj, sisl.WideBandSE), self.elec_SE)):

AttributeError: module 'sisl' has no attribute 'WideBandSE'

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.

I am having troubles when I try to install the Hubbard package on Windows, I got this error message:

(env1) PS C:\Users\jhieulle\Anaconda3\Lib\site-packages\Hubbard> python setup.py install -fcompiler=gfortran -compiler=mingw32
Warning: Assuming default configuration (.\plot/{setup_plot,setup}.py was not found)
Appending Hubbard.plot configuration to Hubbard
Ignoring attempt to set 'name' (from 'Hubbard' to 'Hubbard.plot')
Warning: Assuming default configuration (.\geometry/{setup_geometry,setup}.py was not found)
Appending Hubbard.geometry configuration to Hubbard
Ignoring attempt to set 'name' (from 'Hubbard' to 'Hubbard.geometry')
usage: setup.py [global_opts] cmd1 [cmd1_opts] [cmd2 [cmd2_opts] ...]
or: setup.py --help [cmd1 cmd2 ...]
or: setup.py --help-commands
or: setup.py cmd --help

error: option -o not recognized
(env1) PS C:\Users\jhieulle\Anaconda3\Lib\site-packages\Hubbard>

Does someone already tried to install it on windows?
Thanks a lot for the help.
Best Regards,
Jeremy.

The idea is to have a number of scripts (and possibly a master script calling all individual ones) in the Tests folder to quickly check/test the various developments with little focus on the physics:

• Quick tests for Hubbard and its dependencies
• Detailed tests comparing obtained densities, energies, etc with well-converged, reference results
• Checks for both molecules and periodic systems

The Examples folder could on the other hand showcase some real and interesting physical systems.

Implement option to use a non-orthogonal basis for p_z orbitals representation

The idea is to have the ability to run a calculation of a molecule on a substrate using the negf class, i.e., by including the self-energy of the substrate in the calculation of the molecule, and/or the self-energy of the STM tip.

Tests/test-hubbard.py is buggy

To initialize the spin configuration in SIESTA calculations it could be very useful to use the MFH solution. We would need to implement the writing of a DM.InitSpin block in .fdf format.

Build a common iterative method where only the calculation of the occupations changes for the different calculations, instead of having several iterative methods.

It would be nice to plot the real-space electron and/or spin density using tools from sisl.

The various densities (spin polarization and wavefunction densities) seems to be constructed by squaring first the coefficients and then multiplying onto individual p_z(r) orbitals on the 3D grid. Is this correct? Shouldn't it rather be combined with p_z^2(r)?

When tiling a structure you force an int conversion (which always takes the floor of the floating point).

I would suggest you to do another way since some times 1 electron will disappear.

I.e. in the test-hubbard-open.py if you tile the electrode 2 times, sometimes Nup or Ndn becomes 3 instead of the wanted 4.

There are several independent operations that could be parallelized.
For example, for periodic systems the diagonalization must be done per k-point. However these are independent processes, since the diagonalization for a particular k-point does not depend on other k-points.

We should be able to recover this functionality from our Jupyter notebooks.

I'm really sorry for writing another issue. Feel free to tell me when this is getting to annoying for you. So I have two new problems:

1. I used the function plot.SpinPolarization and it worked totally fine. Now I wanted to to store the maximum value of the spin-polarization. For that I would need an array where all the polarization values are stored. Is there an easy way to get these values?
2. I want to calculate the gap energy for the following band-structure:
   

So in a similar fashion: Is there a way to get the values which are plotted in the band-structure plot?

Sorry again for asking that many questions

So far we have a single-determinant solving method to approach the problem.
We could develop an additional method that builds the wavefunction as a linear combination of Slater-determinants by taking all possible excitations within the active space.

Hi both,

I have now created a utility which reads in TS.ChemPots and TS.Contour.* and then creates the required contours.

E.g. it could be useful for you when running open systems since you can change stuff in an fdf file, run tscontour CONTOUR.fdf and you will get the TSCC* contour outputs.

Perhaps this may be useful?

You can find it here:
https://gitlab.com/nickpapior/siesta/-/archive/tbtrans-store/siesta-tbtrans-store.zip

Util/TS/tscontour

Implement tile and repeat for HubbardHamiltonian?

Currently the filling scheme assumes an electronic gap between valence (fillled) and conduction (empty) bands. To generalize the methods to metallic systems one would need to introduce fractional occupations and determine the Fermi energy.

Hello,
I wanted to calculate the bandgap of an doped 7-AGNR. Therefore I wanted to change the parameters in the sp2-function. I looked in the paper which is referenced in the periodic-example and saw that for the armchair-configuration it might be the best thing to chose the parameters as following:

t1 = 2.7
t2 = 0.09
t3 = 0.27
s1 = 0.11
s2 = 0.045
s3 = 0.026

But as soon as I had an input of 3 non zero-values for s1, s2, s3 I received an error:

ValueError: shape mismatch: value array of shape (3,) could not be broadcast to indexing result of shape (2,). Let me know if I should sent more of my code :)

Hello,
I want to replace a C-atom with a N-atom in the backbone of a AGNR. I only found a way to remove an atom but not replace it.

My code looks like this:
agnr =agnr.remove(0)

When I looked in the xyz-file it consisted only of 3 numbers per row.

The xyz-file looks like this:
x y z

and not like this:
Atom x y z right?

Does this mean that there is no way to distinguish between atomic-species?

So my question basically comes down to these two points:

• Is it even possible to include atoms besides carbon?
• If so, how do I insert impurities?

Note: I constructed everything using the sisl.geom and the .tile-functions.

I just started working with hubbard and sisl so I hope that there is an easy solution to this.

It could be nice to implement in the NEGF class a method to invert the Green's function using the block tridiagonal inversion algorithm to reduce the time that this operation consumes.

Allow iterative3 method for either 1) convergence towards a total number of electrons, and 2) towards a user-defined spin-configuration in open-systems.

Currently this repo is private. Before the workshop I suspect it to be published, or do you want me to ship a zip file?

For more generality, we could use an instance of sisl.DensityMatrix instead of the current nup and ndn arrays.

Introduce self-energies from electrodes etc. NEGF formalism.

The averaging of electron density over the first Brillouin zone does currently not make use of any symmetries. Time-reversal symmetry, i.e., that H(k) = H(-k) could rather easily be implemented and speed up calculations for periodic systems by roughly a factor 2.

Here you should do:

edges = H.edges(ia, exclude=ia)
for ib in edges:

easier (and faster) ;)

Currently, with the recent restructuring of the Hubbard code, some scripts in the examples folder need an overhaul, e.g., examples/molecules/kondo-paper/fig_S11-S13.py that is not producing the relevant plots. General updates are needed to bring all examples up to date with the current version of the code.

For curved carbon structures we would generally need to include (s, px, py, pz) as basis functions. This would in turn imply a list of U-values, one for each kind of orbital, say [U(C_s), U(C_p)] for the carbon orbitals.

Further, we could think of different elements (C, H, N, B, etc), each type with different onsite repulsion parameters per orbital.

Currently we have two independent instances of sisl.Hamiltonian for the two two spin channels. We should rather try to find out how to take advantage of the spin-index in sisl (and thus work with just one Hamiltonian instance).

Function to plot the projected DOS

Tersoff-Hamann theory with both s and p-wave tips?

Instead of using first- or third-nearest neighbour models (1NN/3NN) it would be nice to construct the Hamiltonian from Slater-Koster type parametrizations.

Add these molecular structures in Examples/molecules/triangulene and Examples/molecules/clar-goblet, respectively.