dean0x7d / pybinding

Scientific Python package for tight-binding calculations in solid state physics
https://pybinding.site
BSD 2-Clause "Simplified" License
188 stars 68 forks source link

What does `translational_symmetry` actually perform? #60

Open IlGufo opened 1 month ago

IlGufo commented 1 month ago

I don't understand properly the kind of modifications that translational_symmetry applies to the model. I tried to read the source code, but it is linked to some C++ code which is obscure to me. Does it perform some kind of Fourier transform of the Hamiltonian?

BertJorissen commented 1 month ago

The translation symmetry operator makes periodic boundary conditions at the edges of the model. This is by default given as the length of the unit cell, but it can change by giving the multiple of the unit vectors. If you omit this, the Hamilton is just one for a flake and this can not have the 'k' wave vector as it is a finite system. A finite system could only have real variables if there are none present in the onsite or hopping parameters, which can speed up the calculation.

The recent update on 'pip install pybinding-dev' has the e^(+ik.r) factor, but keep in mind that you have to give the hermitian conjugate of the required hopping when defining the hoppings.

In the older version (the one from 'pip install pybinding'), the changes in the Hamiltonian are that the hoppings that go over the 'boundary' (connect to the other side) get an additional phase of e^(+ir.r), with r the length of the length of the period. In the newer version (the one from 'pip install pybinding-dev'), there is an additional option that you can give to the pb.Model(), which is 'pb.force_phase()'. This will give every hopping the factor e^(ik.r), with r the length of the hopping. This last addition can be usefull when you compare the wave function for different k positions as otherwise the basis for this wave function will change depending on the k wave vector.

IlGufo commented 1 month ago

Thank you very much for your reply. Let me reformulate with some math what I have understood.

We can start from the Bloch sum

\Phi_{\mu i}(\vec{k}, \vec{r}) = \frac{1}{\sqrt{N}}\sum_{\vec{t}}e^{i\vec{k}\cdot(\vec{t}+\vec{\delta}_{\mu})}\phi_{\mu i}(\vec{r}-\vec{\delta}_{\mu}-\vec{t})

$\mu$ is the site index, $i$ is the orbital index, $\vec{\delta}{\mu}$ the position of site $\mu$, $N$ the number of periodic cells and $\vec{t}$ the Bravais vectors. $`\phi{\mu i}`$ are the orbital wave functions. In case the translation periodicity is absent, we simply have $N=1$ with $\vec{t}=\vec{0}$ and $\vec{k}$ can be put to $\vec{0}$ as well.

In the following I denote as $\langle \phi_{\mu i}(\vec{r}-\vec{\delta}_{\mu})|\hat{H}|\phi_{\nu i}(\vec{r}-\vec{\delta}_{\nu})\rangle|^{\text{hopping}}$ the hopping energies which are of course not computed using the average but are given i.e. with @hopping_energy_modifier and depend upon site positions $\mu$, $\nu$ and are matrices.

No periodicity

Hamiltonian is

H_{\mu i \nu j} = E^{\text{on-site}}_{ij}\delta_{\mu \nu} + \langle \phi_{\mu i}(\vec{r}-\vec{\delta}_{\mu})|\hat{H}|\phi_{\nu i}(\vec{r}-\vec{\delta}_{\nu})\rangle|^{\text{hopping}}

Old version

After adding translational_symmetry(), the Hamiltonian became

H_{\mu i \nu j} = E^{\text{on-site}}_{ij}\delta_{\mu \nu} +\sum_{\vec{t}}e^{i\vec{k}\cdot\vec{t}} \langle \phi_{\mu i}(\vec{r}-\vec{\delta}_{\mu})|\hat{H}|\phi_{\nu i}(\vec{r}-\vec{\delta}_{\nu}-\vec{t})\rangle|^{\text{hopping}}

where $\vec{k}$ is set with set_wave_vector() and the sum is over the cells near the primitive one. Solving the system we get the energy bands but the eigenstates are implicitly $\vec{k}$ dependent, as their basis are the Bloch sums and not the orbital states.

New version

After adding force_phase() and translational_symmetry(), the Hamiltonian became

H_{\mu i \nu j} = E^{\text{on-site}}_{ij}\delta_{\mu \nu} +\sum_{\vec{t}}e^{i\vec{k}\cdot\vec{t}}e^{i\vec{k}\cdot(\vec{\delta}_{\nu}-\vec{\delta}_{\mu})} \langle \phi_{\mu i}(\vec{r}-\vec{\delta}_{\mu})|\hat{H}|\phi_{\nu i}(\vec{r}-\vec{\delta}_{\nu}-\vec{t})\rangle|^{\text{hopping}}.

While this is the exact representation starting from the definition of Bloch sum, I can't figure out why the eigenstates should be independent upon $\vec{k}$.

IlGufo commented 1 month ago

The recent update on 'pip install pybinding-dev' has the e^(+ik.r) factor, but keep in mind that you have to give the hermitian conjugate of the required hopping when defining the hoppings.

Isn't the Hermitian conjugate added automatically?

BertJorissen commented 1 month ago

As you pointed out, tight-binding does not use the normal orbitals $\phi(\vec r)$ with a special dependence, but projects everything on these orbitals which forms the basis for the Hamiltonian, $(\phi_a, \phi_b)^T=|\phi>$. Your definition of the hopping term is thus certainly correct, but the spatial component $\vec r$ is not directly present, only via the location of the orbital. The only remark I want to point out is that your sum should run up to infinity, as the k-vector is typically associated with infinite systems. In a large system, the Brillouin zone gets folded until you end up with a gamma-point calculation. This infinite part in the tight-binding method is given by this periodic boundaries which thus give rise to the k-vector.

The k-vector is set in the pb.Model object. You can manually set the k-vector with model.set_wave_vector() or use the solver to calculate the eigenvalues (energies) or eigenvectors (wavefunction) along that path with solver.calc_bands(k_path) and solver.calc_wavefunction(k_path).

I meant with the Hermitian conjugate that the upper-right component of the matrix is filled in with the object that you pass, the lower left is the Hermitian conjugate automatically filled in with the Hermitian conjugate of thevcalue you passed.

You can also see the most recent documentation at bertjorissen.github.io/pybinding.

The translation symmetry is added in the normal way by adding the phase $e^{i\vec k \cdot \vec\delta}$ to the red hoppings like the hopping below. The term $\vec\delta$ is the length of the supercell. This basis however is dependent of the k-vector you set at the pb.Model object, thus you cannot compare two eigenvectors with each other. The pb.force_phase sets the guage of the Hamiltonian so that you can compare the eigenvectors from two different k-vectors. The factor $e^{i\vec k \cdot \vec\delta}$ is added for each hopping, with $\vec\delta$ the length of that hopping. image