2.8. Advanced manybody settings¶

We show some advanced settings for solving the manybody Schrödinger equation which is based on Solving the many-body problem. In the following, syst is a finalized kwant system with leads and we import the manybody module from tkwant:

from tkwant import manybody

2.8.1. Lead occupation¶

The occupation of the manybody state depends on the temperature \(T\), the chemical potential \(\mu\), as well as the distribution function \(f\) of the leads. The function manybody.lead_occupation() offers an easy way to set the lead occupation. This is possible both in the “high-level” approach with manybody.State() and also in the “low-level” approach with manybody.WaveFunction(), see Solving the many-body problem. The default occupation

occupation = manybody.lead_occupation()

corresponds to \(T = \mu = 0\) and noninteracting Fermi-Dirac distribution function. In the case that several leads are attached to the system, the occupation can be a sequence of occupation objects, one for each lead. If only one occupation object is used, as in this example, the occupation is assumed to be identical for each lead.

In the following we show how to change the default behavior of manybody.lead_occupation().

Chemical potential¶

The chemical potential (\(\mu\)) is zero by default. In the following example, we set the chemical potential of the lead to the finite value \(\mu = 1\):

occupation = manybody.lead_occupation(chemical_potential=1.)

For zero temperature and Fermi Dirac distribution, the chemical potential is identical to the Fermi energy.

Temperature¶

The temperature (\(T\)) is zero by default. In the following example, we set the temperature of the lead to the finite value \(T = 0.5\):

occupation = manybody.lead_occupation(temperature=0.5)

Distribution function¶

The distribution function \(f_\alpha(E)\) can be changed with the keyword distribution. As an example, we set up a lead with Bose distribution at finite temperature:

def bose_function(mu, T, energy):
    return 1 / (exp((energy - mu) / T) - 1)

occupation = manybody.lead_occupation(temperature=0.5, distribution=bose_function)

Note that the distribution function \(f\) must have the calling signature (chemical_potential, temperature, energy).

Energy range¶

Upper and lower energy values can be set with the keyword energy_range, which is a sequence of intervals in the form (emin, emax). If present, only modes with energies emin \(\leq E_n \leq\) emax are considered. As an example, we select only the modes with energies \(0.4 \leq E_n \leq 0.6\):

occupation = manybody.lead_occupation(energy_range=[(0.4, 0.6)])

Include / exclude bands¶

The statistical average can be performed only over a subset of energy bands \(E_n\). Specifying the keyword bands with a list of band indices, only the bands with band index \(n\) specified in the list are included. As an example, we perform the statistical average only over the two lowest energy \(E_n\) bands with band index \(n = 0\) and 1:

occupation = manybody.lead_occupation(bands=[0, 1])

Note that the provided band indicees must be in the physically valid range, that is, they must not exceed the maximum number of bands.

Include / exclude leads¶

For system with several leads, the occupation passed to manybody.State() can be a sequence with one element per leads. If a lead should not contribute to the statistial average, the corresponding element of the occupation sequence must be set to None. In the following example, our system is expected to has two leads, but only the contribution of the lead with index 0 is taken into account:

occup = manybody.lead_occupation()
occupations = [occup, None]

Occupation data format¶

The function manybody.lead_occupation() returns a manybody.Occupation instance with the physical information on how to perform the statistical average. One can directly inspect the information stored:

occupation = manybody.lead_occupation()
print(occupation)

lead occupation: distribution=<function lead_occupation.<locals>._zero_temperature_fermi_dirac_distribution at 0x7ff9de031158>, energy_range=[(None, 0)] , bands=None

The occupation object is used to extract energy cutoffs in order to calculate quadrature intervals and boundary conditions. It also stores the distribution function \(f(E)\).

2.8.2. Numerical integration¶

Quadrature intervals for the manybody integral are calculated from the lead occupation and the lead spectrum:

spectra = kwantspectrum.spectra(syst.leads)
occupation = manybody.lead_occupation()
intervals = manybody.calc_intervals(spectra, occupation)

This sequence is part of the pre-processing in the “low-level” approach, see Solving the many-body problem. In the “high-level” approach, one can precalculate the intervals and pass then to the manybody state, in order to bypass the default interval calculation in manybody.State():

tmax = 10
spectra = kwantspectrum.spectra(syst.leads)
occupation = manybody.lead_occupation()
intervals = manybody.calc_intervals(spectra, occupation)
state = manybody.State(syst, tmax, occupation, intervals=intervals)

This mechanism is especially useful to directly manipulate intervals in the intervals list. A second way to manipulate the interval calculation of manybody.State() is to change the default interval type. The interval type can be passed also by the interval argument:

state = manybody.State(syst, tmax, occupation, intervals=manybody.Interval)

manybody.Interval is the default data class to construct the quadrature interval. One can change the default to alter the behavior, as will be shown in the following.

Quadrature order¶

The quadrature order can be changed via:

import functools
interval_type = functools.partial(manybody.Interval, order=10)
intervals = manybody.calc_intervals(spectra, occupation, interval_type=interval_type)

The list of intervals can be used further in the “low-level” approach or passed to manybody.State() in the “high-level” approach:

state = manybody.State(syst, tmax, occupation, intervals=intervals)

Alternatively, for the second way of the “high-level”, the modified manybody.Interval data class can passed directly to the manybody state:

import functools
interval_type = functools.partial(manybody.Interval, order=10)
state = manybody.State(syst, tmax, occupation, intervals=interval_type)

The order is usually taken between 10 and 20. If the integration is not accurate enough one should rather divide each interval into subintervals with the keyword number_subintervals.

Interval subdivision¶

Momentum split Each momentum quadrature interval can be equidistantly divided into \(n\) subintervals to increase the numerical accuracy. Here we divide into \(n=10\) subintervals with the keyword number_subintervals:

# (kmin, kmax) -> [(kmin, k_0), (k_0, k_1).. (k_{n-1}, kmax)]
intervals = manybody.split_intervals(intervals, number_subintervals=10)

Energy split Splitting intervals in energy space is possible by using energy_range. The snippet below is to reproduce older calculations performed with energy integrations on an energy interval splitted equidistantly into subintervals (here 10).

from itertools import tee

def pairwise(iterable):
    "s -> (s0,s1), (s1,s2), (s2, s3), ..."
    a, b = tee(iterable)
    next(b, None)
    return zip(a, b)

chemical_potential = 1

# energy interval split (0, 1) -> [(0, 0.1), (0.1, 0.2).. (0.9, 1)]
energy_range = [i for i in pairwise(np.linspace(0, chemical_potential, 11))]
occupation = manybody.lead_occupation(chemical_potential,
                                      energy_range=energy_range)

Interval = functools.partial(manybody.Interval, integration_variable='energy')
intervals = manybody.calc_intervals(spectra, occupation, interval_type=Interval)

intervals is now already a sequence of 10 intervals, no additional split with manybody.split_intervals is needed. If integration_variable="momentum", the split is still performed in the energy range, but the quadrature will be integrate over momentum.

Quadrature rule¶

The quadrature rule can be changed via:

import functools
interval_type = functools.partial(manybody.Interval, quadrature = 'gausslegendre')
intervals = manybody.calc_intervals(spectra, occupation, interval_type=interval_type)

Note

Changing the quadrature rule is only useful in the low-level approach. In manybody.State() the error estimate and the refinement is based on Gauss-Kronrod quadrature, which cannot be changed. If intervals with another rule than Gauss-Kronrod are passed to manybody.State(), these intervals do not take place for refinement and error estimate.

Quadrature error¶

The solver manybody.State() has a method estimate_error() to estimate the quadrature error. Here we show an alternative error estimate in the low-level approach. We define the integration error (\(\delta\)) as:

\[\delta = \text{max}(|\rho_{n} - \rho_{2n+1}|)\]

where \(\rho_{n}\) is the density calculated with the lower (\(n\)) order rule and \(\rho_{2 n + 1}\) is the density calculated with the higher (\(2 n +1\)) quadrature rule of the Gauss-Kronrod method.

def maximal_absolute_error(result):
    low_order, high_order = result
    return np.max(np.abs(low_order - high_order))

intervals = manybody.calc_intervals(spectra, occupation)
tasks = manybody.calc_tasks(intervals, spectra, occupation)
emin, emax = manybody.calc_energy_cutoffs(occupation)
boundaries = tkwant.leads.automatic_boundary(spectra, tmax, emin=emin, emax=emax)
psi_init = manybody.calc_initial_state(syst, tasks, boundaries)
state = manybody.WaveFunction(psi_init, tasks)
density = state.evaluate(density_operator)
print('integration error delta= {}'.format(maximal_absolute_error(density)))

integration error delta= 3.0679847551340345e-11

The quadrature rule applied to each interval can be accessed using the keyword quadrature. Note that the adaptive state manybody.State needs a Gauss-Kronrod like rule for the error estimate.

Energy vs. momentum integration¶

The expectation value of an observable \(\hat{\mathbf{A}}\) in the manybody state is defined as

\[\langle \hat{\mathbf{A}} \rangle (t) = \sum_{\alpha} \int \frac{dE}{2 \pi} f_\alpha(E) \psi_{\alpha E}^*(t,i) \mathbf{A}_{ij} \psi_{\alpha E}(t,j)\]

The integrand of the above energy integral diverges weakly as each scattering wave function scales as \(\sim 1/\sqrt{v_\alpha(E)}\), with velocities \(v_\alpha(E) = \frac{d E_{\alpha}}{d k}\) in the vicinity of the band openings. We can rewrite the integral in the form

\[\langle \hat{\mathbf{A}} \rangle (t) = \sum_{\alpha} \int_{k_{{\rm min}, \alpha}}^{k_{{\rm max}, \alpha}} \frac{dk}{2 \pi} f_\alpha(E_\alpha(k)) v_\alpha(k) \psi_{\alpha k}^*(t,i) \mathbf{A}_{ij} \psi_{\alpha k}(t,j)\]

which is analytically equivalent but performs better numerically, as the diverging points are eliminated.

One can switch between the two ways to perform the manybody average with the keyword integration_variable. By default, integration_variable = momentum, and the integration is performed over the momentum, which corresponds the second equation. With integration_variable = energy one switches to the first equation and integrates explicitly over

Intervals can be pre-calculated with

import functools
interval_type = functools.partial(manybody.Interval, integration_variable='energy')
intervals = manybody.calc_intervals(spectra, occupation, interval_type=interval_type)

Note that manybody.Intervals returned from manybody.intervals() always store momentum intervals, independent of integration_variable.

Alternatively, the integral type used by the manybody state is changed via

interval_type = functools.partial(manybody.Interval, integration_variable='energy')
state = manybody.State(syst, tmax, occupation, intervals=interval_type)

Interval data format¶

One can directly inspect the stored information of a quadrature interval:

intervals = manybody.calc_intervals(spectra, occupation)
for interval in intervals:
    print(interval)

quadrature interval: lead=0, band=0, kmin=1.047198, kmax=1.104031, order=10, quadrature=kronrod, integration_variable=momentum
quadrature interval: lead=0, band=0, kmin=1.104031, kmax=1.159279, order=10, quadrature=kronrod, integration_variable=momentum
quadrature interval: lead=0, band=0, kmin=1.159279, kmax=1.213225, order=10, quadrature=kronrod, integration_variable=momentum
quadrature interval: lead=0, band=0, kmin=1.213225, kmax=1.266104, order=10, quadrature=kronrod, integration_variable=momentum
quadrature interval: lead=0, band=0, kmin=1.266104, kmax=1.318116, order=10, quadrature=kronrod, integration_variable=momentum
quadrature interval: lead=0, band=0, kmin=1.318116, kmax=1.369438, order=10, quadrature=kronrod, integration_variable=momentum
quadrature interval: lead=0, band=0, kmin=1.369438, kmax=1.420228, order=10, quadrature=kronrod, integration_variable=momentum
quadrature interval: lead=0, band=0, kmin=1.420228, kmax=1.470629, order=10, quadrature=kronrod, integration_variable=momentum
quadrature interval: lead=0, band=0, kmin=1.470629, kmax=1.520775, order=10, quadrature=kronrod, integration_variable=momentum
quadrature interval: lead=0, band=0, kmin=1.520775, kmax=1.570796, order=10, quadrature=kronrod, integration_variable=momentum

lead and band corresonds to the lead respectively band index, and kmin and kmax are the lower respectively upper momentum values of the interval. The actual values change depending on the specific problem. The other parameters have default values: order is the quadrature order applied to the interval (generally the number of point at which the interval is discretized) and quadrature is the concrete quadrature rule that will be applied. integration_variable is an additional information to switch between energy and momentum integration.

Tasks data format¶

The routine calc_tasks() calculates the modes and weights for all intervals:

spectra = kwantspectrum.spectra(syst.leads)
occupation = manybody.lead_occupation()
intervals = manybody.calc_intervals(spectra, occupation)
tasks = manybody.calc_tasks(intervals, spectra, occupation)

To be more precise, calc_tasks() returns a dictionary with information about all one-body states and the corresponding weighting factor. Each element in the dictionary corresponds to a onebody state and looks like this:

for key, task in tasks.items():
    if key < 3:  # print only the first three tasks
        print(task)

onebody task: lead=0, mode=0, energy=-0.9999948293377526, momentum=0.002273904121928405, weight=[0.00000000e+00 4.43209212e-06], math_weight=[0.        0.0061233], phys_weight=0.0007238079361937524
onebody task: lead=0, mode=0, energy=-0.9998133386878931, momentum=0.013662509717224935, weight=[1.51811637e-04 7.41354171e-05], math_weight=[0.03490903 0.01704741], phys_weight=0.004348777912409993
onebody task: lead=0, mode=0, energy=-0.9986628241791404, momentum=0.03656945200472678, weight=[0.         0.00033366], math_weight=[0.         0.02867012], phys_weight=0.011637821525310944

For each set of (lead, mode, energy) attributes which are stored in each task element, the function manybody.calc_initial_state() calculates the corresponding one-body scattering state \(\psi_{\alpha}(t=0, x)\):

tasks = manybody.calc_tasks(intervals, spectra, occupation)
boundaries = tkwant.leads.automatic_boundary(spectra, tmax)
psi_init = manybody.calc_initial_state(syst, tasks, boundaries)

that form the manybody state at the initial time t=0. The initial state can be used as an initial state of manybody wave function manybody.WaveFunction():

psi_init = manybody.calc_initial_state(syst, tasks, boundaries)
state = manybody.WaveFunction(psi_init, tasks)

Evaluating an observable over the manybody state, the state manybody.WaveFunction() will sum over all one-body states present in the tasks list and weight each term in the sum with the respective weight weight.

Note

Tasks can not passed to the “high-level” solver manybody.State. This is because the quadrature error (via estimate_error()) can be estimated only on intervals, not on individual evaluation points. Moreover, psi_init is a dictionary, but which is distributed over all MPI ranks. We refer to the source code for technical details.

Time integration¶

The time integration is performed at the level of onebody states. One can change the default settings by prebinding values with the module functool.partial to the onebody state. In the current example, we change the relative tolerance rtol of the time-stepping algorithm. The first binding to solver_type changes the tolerance rtol of the time integrator. The second binding defines a new onebody solver onebody_solver that uses this time integrator. The onebody solver can be passed via keyword to one of the three solvers. Here we show the example for manybody.State():

import functools
solver_type = functools.partial(tkwant.onebody.solvers.default, rtol=1E-5)
onebody_wavefunction_type = functools.partial(tkwant.onebody.WaveFunction.from_kwant,
                                              solver_type=solver_type)
scattering_state_type = functools.partial(tkwant.onebody.ScatteringStates,
                                          wavefunction_type=onebody_wavefunction_type)
state = manybody.State(syst, tmax=10, scattering_state_type=scattering_state_type)

A similar strategy is possible to change the onebody kernels onebody.kernels that evaluate the right-hand-side of the one-body Schrödinger equations.

2.8.3. Retrieve one-body states¶

The manybody state manybody.WaveFunction provides several methods to retrieve the underlying one-body states. They can be useful for realizing specific problems beyond the examples mentioned above.

State identifier¶

The attribute state.tasks is a dictionary of all one-body states stored in the state instance. The state.tasks dictionary is identical on all MPI ranks. The key (\(\equiv \alpha\)) uniquely labels each one-body state \(\psi_{\alpha}(t)\) :

state = manybody.WaveFunction(psi_init, tasks)
keys = state.get_keys()
print('number of one-body states inside manybody = {}'.format(len(keys)))
print('keys={}'.format(keys))

number of one-body states inside manybody = 21
keys=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]

State information¶

Further information about the one-body state \(\psi_{\alpha}(t)\) can be found in the state.tasks dictionary. It contains the lead and the mode index, the mode energy, as well as the weighting factor \(w_{\alpha, l}\) in the statistical average. As an example, we print the information for one-body state corresponding to the key zero:

key = 0
print(state.tasks[key])

onebody task: lead=0, mode=0, energy=-0.9999948293377526, momentum=0.002273904121928405, weight=[0.00000000e+00 4.43209212e-06], math_weight=[0.        0.0061233], phys_weight=0.0007238079361937524

One-body states¶

One-body states \(\psi_{\alpha}(t)\) can be retrieved using the get_state() method:

psi = state.get_onebody_state(key=0)

Note that the attribute state.psi is a dictionary of one-body states that basically provides the same information. In parallel MPI environments, however, some caution is required. When multiple ranks are available manybody.State() distributes the one-body states to all MPI ranks. Since each one-body state is unique, a particluar one-body state is located only on one of all the available MPI ranks. Although not currently implemented, manybody.State() could additionally use dynamic load balancing in the future and redistribute the one-body states at runtime. The get_state() method guarantees that the requested state is always returned without worring about the MPI rank distribution.

In addition, we would like to point out that the result of evaluating an observable is independent from the distribution of the one-body states on the MPI ranks (although the result is not bitwise identical due to rounding errors), because the statistical average is taken over all states listed in state.tasks. Only the performance and memory consumption might become bad, if the distribution to the MPI ranks is not balanced.

Evolution time¶

The current time of the state can be obtained via the attribute state.time:

print('time= {}'.format(state.time))

time= 0

The initial time of the manybody wavefunction manybody.WaveFunction is zero, so that state.time=0 after initialization. It is not checked if the time attributes in the individual one-body states are similar or consistent, only backward-propagation of the states in time is not allowd. After calling state.evolve(time) all one-body states that compose the manybody wavefunction have evolved to time time and also state.time equals time.

2.8.4. Miscellaneous¶

Band structure¶

The band spectra for all leads are calculated on demand inside manybody.State(). An own spectrum can be passed to the state with the keyword spectra:

spectra = kwantspectrum.spectra(syst.leads)
state = manybody.State(syst, tmax=10, spectra=spectra)

Boundary conditions¶

Boundary conditions are calculated internally for all leads of the system. In order to provide own boundary conditions, they can be precalculate and transferred to the state with keyword boundaries:

boundaries = [tkwant.leads.SimpleBoundary(tmax=10)] * len(syst.leads)
state = manybody.State(syst, boundaries=boundaries)

Initial state¶

Arbitrary initial states can be provided to the manybody state with the keyword psi_init. The state is simply a dictionary of one-body states. In addition, a dictionary with the weighting factors for each state must be provided with the keyword tasks:

new_state = manybody.WaveFunction(psi_init=psi_init, tasks=tasks)

Note that the boundary condition is attached to each (initial) one-body state that forms the manybody state. Therefore, one cannot redefine boundary conditions by initializing manybody.State() from psi_init. Passing the keywords boundary and/or tmax together with psi_init leads to a ValueError.

The time attribute new_state.time=0 after initialization. As stated above, is not checked if the time attributes in the individual one-body states are similar or consistent.

Boundstates¶

The following code illustrates how to include boundstates in the manybody dynamics. First, the boundstates are calculated in some user-written function, here called calc_my_boundstates():

boundstate_psi, boundstate_tasks = calc_my_boundstates()

boundstate_psi is a dictionary of all boundstate wavefunctions. It has basically the same form as an initial manybody state calculated from manybody.calc_initial_state(). boundstate_tasks contain the weighting factor and the energy of each boundstate and is similar to the output of manybody.calc_tasks().

In the case that state is a manybody state instance created with manybody.State, boundstates are easily added by the method add_boundstates():

state = manybody.State(syst, tmax)
state.add_boundstates(boundstate_psi, boundstate_tasks)

In the case that the manybody state instance has been created with manybody.WaveFunction, the boundstates must be added with the function manybody.add_boundstates(). If the boundstates are not yet time dependent, they must first transformed with make_boundstate_time_dependent. In both cases, the keys for the new boundstates must not be present in the manybody state already.