mnt: enabled tiling an eigenstate #355

zerothi · 2021-07-07T13:02:34Z

Using the Bloch expansion we can get the same state for
a new lattice by expanding it.

Closes Can one also get the state vector(s) in the *extended* cell in a similar way? #354
Tests added
Documentation for functionality
User visible changes (including notable bug fixes) are documented in CHANGELOG

@tfrederiksen could you have a look at this

tfrederiksen · 2021-07-07T15:04:06Z

How does one use tile to construct a given state in the supercell?

I was thinking something like:

es = H.eigenstate(k)
es_ext = np.concatenate([np.exp(2j * np.pi * k.dot(isc)) * es.state for _, isc in H.geom.sc], axis=1)

where one could then compute matrix elements (in this case band energies) in "real space" for the operator/matrix like

ev1 = np.diag(es.state.conj().dot(H.tocsr().dot(es_ext.T)))
ev2 = es.inner(matrix=H.Hk(k))
assert np.allclose(ev1, ev2)

zerothi · 2021-07-07T18:47:53Z

How does one use tile to construct a given state in the supercell?
I think you can just do:

es = H.eigenstate(...)
es_nsc = es
for axis, nsc in enumerate(H.geometry.sc.nsc):
   es_nsc = es_nsc.tile(nsc, axis)

# possibly one should re-arrange data in some other way to make them *diagonal*
es.inner(es_nsc.state.T, H.Hk(...)

and then you have the full eigenstate vectors in a geometry as large as the supercell, and the phase-factors will run from [0 ; nsc-1] instead of [-nsc//2; nsc//2]. But I think this is just a wrap-around and would be equivalent (at the top of my head, could be wrong)?

I was thinking something like:
es = H.eigenstate(k)
es_ext = np.concatenate([np.exp(2j * np.pi * k.dot(isc)) * es.state for _, isc in H.geom.sc], axis=1)
where one could then compute matrix elements (in this case band energies) in "real space" for the operator/matrix like
ev1 = np.diag(es.state.conj().dot(H.tocsr().dot(es_ext.T)))
ev2 = es.inner(matrix=H.Hk(k))
assert np.allclose(ev1, ev2)

Also, these big calculations will quickly take up lots of memory, for a system with 1000 orbitals and nsc = [3, 3, 1] you will find that doing this for a single ev1 takes up 1000*1000 * 3 * 3 reals (68 MB).

Wouldn't the above method suffice?

Using the Bloch expansion we can get the same state for a new lattice by expanding it.

zerothi · 2021-07-08T06:55:25Z

I have now amended the code to have an offset parameter, this will enable you to do es.tile(3, 1, offset=1) which then uses the phases [-1,0,1].

tfrederiksen · 2021-07-08T07:08:59Z

Nice! Could we have a property to an eigenstate object, like es.state_sc, that automatically calls tile with the correct offsets to have the state vectors in the supercell?

Reordering seems also necessary before one can do the dot product H.tocsr().dot(es.state_sc)

zerothi · 2021-07-08T08:09:15Z

Nice! Could we have a property to an eigenstate object, like es.state_sc, that automatically calls tile with the correct offsets to have the state vectors in the supercell?

Hmm. I am a little bit hesitant since you can't really use it as you imagine it, i.e. the problem of re-ordering. And if one has a hard time figuring out what it actually means, then probably not so useful.. ?

Reordering seems also necessary before one can do the dot product H.tocsr().dot(es.state_sc)

Yeah, this makes it rather problematic. I have now added a translate method so that one can do:

es_exts = map(lambda x: es.translate(x[1]), H.geometry.sc)
np.concatenate((es.inner(ext) for ext in es_exts), axis=1)

would this suffice?

It seems this is a bit more natural?

zerothi · 2021-07-08T09:59:41Z

I have now amended the code to have an offset parameter, this will enable you to do es.tile(3, 1, offset=1) which then uses the phases [-1,0,1].

Actually, the offset should not be subtractive, I'll fix this, the above should be es.tile(3, 1, offset=-1)

tfrederiksen · 2021-07-08T10:32:19Z

es_exts = map(lambda x: es.translate(x[1]), H.geometry.sc)
np.concatenate((es.inner(ext) for ext in es_exts), axis=1)

would this suffice?

It seems this is a bit more natural?

I would really still like to see this as something the user can invoke directly from the object. I think H.tocsr(-1).dot(es.state_sc) is what we need in ìnner to handle the non-ortho case and bra/ket at different k-points.

Said more generally, I think there are important use cases for inner with matrix being a sparse matrix defined for the supercell.

Couldn't you make it a property?

    @property
    def state_sc(self):
        r""" Eigenstate(s) in the extended supercell """
        return np.concatenate([np.exp(2j * np.pi * self.k.dot(isc)) * self.state for _, isc in self.geom.sc], axis=1)

zerothi · 2021-07-08T10:42:45Z

I would really still like to see this as something the user can invoke directly from the object. I think H.tocsr(-1).dot(es.state_sc) is what we need in ìnner to handle the non-ortho case and bra/ket at different k-points.

Couldn't you make it a property?
    @property
    def state_sc(self):
        r""" Eigenstate(s) in the extended supercell """
        es_exts = map(lambda x: es.translate(x[1]), H.geometry.sc)
        return np.concatenate((es.inner(ext) for ext in es_exts), axis=1)

I guess you mean:

@property
def state_sc(self):
    return np.concatenate([self.translate(isc).state for _, isc in H.geometry.sc], axis=1)

?

However, it still isn't natural to use, I think. It requires bookkeeing to figure out where what belongs, and since you can directly use the above in an es.inner(state_sc) then its use is rather limited.

Also, perhaps another way to use the state would be:

@property
def state_sc(self):
    return np.stack([self.translate(isc).state for _, isc in H.geometry.sc])

which would return (n_s, len(es), es.shape[1]), this would have the correct order by doing, state_sc.reshape(-1, es.shape[1]), I think

But... Is this intuitive? What does sc mean in a State object? Due to its periodicity it would be indefinite.

tfrederiksen · 2021-07-08T10:48:47Z

Sorry, yes, I copied the wrong piece of code.

The book-keeping is already taken care of by adhering to the supercell ordering, no?

From my perspective, being able to evaluate A.tocsr().dot(es.state_sc) seems important (and, from the user perspective, should be as simple as that).

Regarding the wording, my logic was that as state provides the state vectors in the primitive cell, state_sc would denote the corresponding state vectors with Bloch phases and ordering corresponding to the supercell.

This allows one to retrieve the state-vectors in the auxiliary supercell. This is important for further post-processing.

zerothi · 2021-07-08T11:11:15Z

Ok, I have added what should work.

Before merging, one final question:

Currently the ket argument of inner is transposed internally, does this make sense? Wouldn't one expect one to pass the correct order?

tfrederiksen · 2021-07-08T11:31:26Z

sisl/physics/state.py

+        else:
+            def translate(isc):
+                return self.state * exp(2j*_pi * k @ isc)
+        return np.concatenate([translate(isc) for _, isc in sc], axis=0)


Shouldn't it be axis=1 here? Such that self.state_sc.shape == (len(state), len(state) * state.parent.n_s)?

The first index should naturally be the band index as for self.state. The number of bands dnoes not change by extending to the supercell. The second index the orbital index which grows.

No, that would leave the ordering wrong. You would have to increase the number of bands to nbands * n_s to get this correct. Alternatively it should return (n_s, nbands, ncoeffs)?

or (nbands, n_s, ncoeffs)

actually, I am here more in favour of (n_s, nbands, ncoeffs)

I mean you want to do es.inner(es.state_sc), right, then you would have to increase the number of bands, and not number of coefficients, no?

Yes ok, so what you want is to use a non-square matrix as the operator. However, I don't really follow what you then get, since [i, J] . [m, J].T -> [i, m] and becomes a square matrix again, and you loose the supercell information? No?

Note: generally you shouldn't use tocsr for this, rather use H.Hk(format='sc:csr'), this will ensure to be working for NC and SOC as well.

It seems to me that this is a greater problem than what this PR proposes to do, I think I should delete the state_sc, and then another PR can take care of that, then the discussion can continue there, and more ftitingly so I think. :)

What I want is basically the matrix elements between primitive cell band indices n, m but where the operator A is given in "real space". I do not want any phases applied to A.

As far as I see this is also the correct approach to do the inner product between states at different k-points with non-ortho basis, i.e.,

<k1(n) | k2(m)> = k1.state[n, i].S[i, J].((k2.state_sc[m, J]).T)

becomes a square matrix again, and you loose the supercell information? No?

That's the point, I do not want to retain information about the supercell. But I need to consider the supercell to perform the inner product in the most general case.

lets continue discussion in #356

zerothi · 2021-07-08T12:27:25Z

This is now merged (the last commit here isn't and will be followed by #356)

zerothi · 2021-07-08T12:35:11Z

For completeness sake, the code left out of this PR is:

    @property
    def state_sc(self):
        r""" Calculate the state coefficients in the auxiliary supercell
        The format of this state vector is usable for `inner` as the ``ket`` argument.
        The order of the state vectors is equivalent to the auxiliary supercell order.
        See Also
        --------
        inner : these state vectors may be used in the inner product
        tile : equivalent form along a single given axis
        translate : a single translation of the state vectors
        Examples
        --------
        The returned state vectors can be formulated like this:
        >>> state = State(...)
        >>> state_sc = np.stack([state.translate(isc)
        ...                      for _, isc in state.parent.sc])
        >>> assert state_sc.shape == (len(state) * state.parent.n_s, state.shape[1])
        """
        k = _a.asarrayd(self.info.get("k", [0]*3))
        sc = self.parent.sc
        if np.allclose(k, 0):
            dtype = self.state.dtype
            def translate(isc):
                return self.state
        else:
            dtype = np.complex128
            def translate(isc):
                return self.state * exp(2j*_pi * k @ isc)
        state = np.empty((sc.n_s, ) + self.shape, dtype=dtype)
        for i, (_, isc) in enumerate(sc):
            state[i] = translate(isc)
        return state.reshape(-1, self.shape[1])

and a test

    def test_eigenstate_state_sc(self, setup):
        # Test of eigenvalues
        R, param = [0.1, 1.5], [0., 2.7]
        H = setup.H.copy()
        H.construct((R, param))

        k = [0] * 3
        # we must select a k that does not fold on
        # itself which then creates degenerate states
        for k1 in [0.5, 1/3]:
            k[1] = k1
            es = H.eigenstate(k)
            state_sc = es.state_sc.reshape(-1, *es.shape)
            for i, (_, isc) in enumerate(H.geometry.sc):
                est = es.translate(isc)
                assert np.allclose(state_sc[i], est.state)

zerothi added 3 commits July 8, 2021 08:47

mnt: enabled tiling an eigenstate

ed39a94

Using the Bloch expansion we can get the same state for a new lattice by expanding it.

mnt: rebased

9221a9d

enh: added offset argument to displace the phases

a039bbc

zerothi force-pushed the 354-tile-state branch from b4b66f6 to a039bbc Compare July 8, 2021 06:54

enh: added translate to the state objects

6f4dc81

enh: clarified offset and added test for offset

ad7ad06

enh: added state_sc for states

50a98f6

This allows one to retrieve the state-vectors in the auxiliary supercell. This is important for further post-processing.

tfrederiksen reviewed Jul 8, 2021

View reviewed changes

zerothi mentioned this pull request Jul 8, 2021

How to do inner products between different k #356

Closed

zerothi closed this Jul 8, 2021

zerothi deleted the 354-tile-state branch July 8, 2021 12:36

This pull request was closed.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

mnt: enabled tiling an eigenstate #355

mnt: enabled tiling an eigenstate #355

zerothi commented Jul 7, 2021 •

edited

Loading

tfrederiksen commented Jul 7, 2021

zerothi commented Jul 7, 2021

zerothi commented Jul 8, 2021

tfrederiksen commented Jul 8, 2021 •

edited

Loading

zerothi commented Jul 8, 2021

zerothi commented Jul 8, 2021

tfrederiksen commented Jul 8, 2021 •

edited

Loading

zerothi commented Jul 8, 2021 •

edited

Loading

tfrederiksen commented Jul 8, 2021

zerothi commented Jul 8, 2021

tfrederiksen Jul 8, 2021

zerothi Jul 8, 2021

zerothi Jul 8, 2021

zerothi Jul 8, 2021

zerothi Jul 8, 2021

zerothi Jul 8, 2021

zerothi Jul 8, 2021

tfrederiksen Jul 8, 2021

tfrederiksen Jul 8, 2021 •

edited

Loading

zerothi Jul 8, 2021

zerothi commented Jul 8, 2021

zerothi commented Jul 8, 2021 •

edited

Loading

mnt: enabled tiling an eigenstate #355

mnt: enabled tiling an eigenstate #355

Conversation

zerothi commented Jul 7, 2021 • edited Loading

tfrederiksen commented Jul 7, 2021

zerothi commented Jul 7, 2021

zerothi commented Jul 8, 2021

tfrederiksen commented Jul 8, 2021 • edited Loading

zerothi commented Jul 8, 2021

zerothi commented Jul 8, 2021

tfrederiksen commented Jul 8, 2021 • edited Loading

zerothi commented Jul 8, 2021 • edited Loading

tfrederiksen commented Jul 8, 2021

zerothi commented Jul 8, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tfrederiksen Jul 8, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

zerothi commented Jul 8, 2021

zerothi commented Jul 8, 2021 • edited Loading

zerothi commented Jul 7, 2021 •

edited

Loading

tfrederiksen commented Jul 8, 2021 •

edited

Loading

tfrederiksen commented Jul 8, 2021 •

edited

Loading

zerothi commented Jul 8, 2021 •

edited

Loading

tfrederiksen Jul 8, 2021 •

edited

Loading

zerothi commented Jul 8, 2021 •

edited

Loading