8. Charge Transport

8.1. Obtain Mobility by Green-Kubo Formula

class renormalizer.transport.kubo.TransportKubo(model: Model, temperature: Quantity, distance_matrix: Optional[ndarray] = None, insteps: int = 1, ievolve_config=None, compress_config=None, evolve_config=None, dump_dir: Optional[str] = None, job_name: Optional[str] = None, thermal_dump_path: Optional[str] = None, properties: Optional[Property] = None)

Calculate mobility via Green-Kubo formula:

\[\mu = \frac{1}{k_B T} \int_0^\infty dt \langle \hat j (t) \hat j(0) \rangle = \frac{1}{k_B T} \int_0^\infty dt C(t)\]

where

\[\hat j = -\frac{i}{\hbar}[\hat P, \hat H]\]

and \(\hat P = e_0 \sum_m R_m a^\dagger_m a_m\) is the polarization operator.

Note

In principle \(\hat H\) can take any form, however only Holstein-Peierls Hamiltonian is well tested.

More explicitly, \(C(t)\) has the form:

\[C(t) = \textrm{Tr}\{\rho(T) e^{i \hat H t} \hat j(0) e^{- i \hat H t} \hat j (0)\}\]

where we have assumed \(\rho(T)\) is normalized (i.e. it is divided by the partition function \(\textrm{Tr}\{\rho(T)\}\)).

In terms of practical implementation, it is ideal if \(\rho(T)\) is split into two parts to (hopefully) speed up calculation and minimize time evolution error:

\[\begin{split}\begin{align} C(t) & = \textrm{Tr}\{\rho(T) e^{i \hat H t} \hat j(0) e^{- i \hat H t} \hat j(0)\} \\ & = \textrm{Tr}\{e^{-\beta \hat H} e^{i \hat H t} \hat j(0) e^{- i \hat H t} \hat j(0)\} \\ & = \textrm{Tr}\{e^{-\beta \hat H / 2} e^{i \hat H t} \hat j(0) e^{- i \hat H t} \hat j(0) e^{-\beta \hat H / 2}\} \end{align}\end{split}\]

In this class, imaginary time propagation from infinite temperature to \(\beta/2\) is firstly carried out to obtain \(e^{-\beta \hat H / 2}\), and then real time propagation is carried out for \(e^{-\beta \hat H / 2}\) and \(\hat J(0) e^{-\beta \hat H / 2}\) respectively to obtain \(e^{-\beta \hat H / 2} e^{i \hat H t}\) and \(e^{- i \hat H t} \hat j(0) e^{-\beta \hat H / 2}\). The correlation function at time \(t\) can thus be calculated via expectation calculation.

Note

Although the class is able to carry out imaginary time propagation, in practice for large scale calculation it is usually preferable to carry out imaginary time propagation in another job and load the dumped initial state directly in this class.

Parameters:

• model (Model) – system information.

• temperature (Quantity) – simulation temperature. Zero temperature is not supported.

• distance_matrix (np.ndarray) –

  two-dimensional array \(D_{ij} = P_i - P_j\) representing distance between the \(i\) th electronic degree of freedom and the \(j\) th degree of freedom. The index is the same with model.e_dofs. The parameter takes the role of \(\hat P\) and can better handle periodic boundary condition. The default value is None in which case the distance matrix is constructed assuming the system is a one-dimensional chain.

  Note

  The construction of the matrix should be taken with great care if periodic boundary condition is applied. Take a one-dimensional chain as an example, the distance between the leftmost site and the rightmost site is \(\pm R\) where \(R\) is the intersite distance, rather than \(\pm (N-1)R\) where \(N\) is the total number of electronic degrees of freedom.

• insteps (int) – steps for imaginary time propagation.

• ievolve_config (EvolveConfig) – config when carrying out imaginary time propagation.

• compress_config (CompressConfig) – config when compressing MPS. Note that if TDVP based methods are chosen for time evolution, compression is necessary when \(\hat j\) is applied on \(e^{-\beta \hat H / 2}\) but no compression is carried out for \(e^{-\beta \hat H / 2} e^{i \hat H t}\).

• evolve_config (EvolveConfig) – config when carrying out real time propagation.

• dump_dir (str) – the directory for logging and numerical result output.

• job_name (str) – the name of the calculation job which determines the file name of the logging and numerical result output.

• thermal_dump_path (str) – the path to load previously thermal propagated initial state if it exists or the path to dump the thermal state if it does not exist. The default value is determined by appending "_impdm.npz" to the path joined by dump_dir and job_name.

• properties (Property) – other properties to calculate during real time evolution. Currently only supports Holstein model.

property auto_corr: ndarray

Correlation function \(C(t)\).

Returns:

1-d numpy array containing the correlation function evaluated at each time step.

property auto_corr_decomposition: ndarray

Correlation function \(C(t)\) decomposed into contributions from different parts of the current operator. Generally, the current operator can be split into two parts: current without phonon assistance and current with phonon assistance. For example, if the Holstein-Peierls model is considered:

\[\hat H = \sum_{mn} [\epsilon_{mn} + \sum_\lambda \hbar g_{mn\lambda} \omega_\lambda (b^\dagger_\lambda + b_\lambda) ] a^\dagger_m a_n + \sum_\lambda \hbar \omega_\lambda b^\dagger_\lambda b_\lambda\]

Then current operator without phonon assistance is defined as:

\[\hat j_1 = \frac{e_0}{i\hbar} \sum_{mn} (R_m - R_n) \epsilon_{mn} a^\dagger_m a_n\]

and the current operator with phonon assistance is defined as:

\[\hat j_2 = \frac{e_0}{i\hbar} \sum_{mn} (R_m - R_n) \hbar g_{mn\lambda} \omega_\lambda (b^\dagger_\lambda + b_\lambda) a^\dagger_m a_n\]

With \(\hat j = \hat j_1 + \hat j_2\), the correlation function can be decomposed into four parts:

\[\begin{split}\begin{align} C(t) & = \langle \hat j(t) \hat j(0) \rangle \\ & = \langle ( \hat j_1(t) + \hat j_2(t) ) (\hat j_1(0) + \hat j_2(0) ) \rangle \\ & = \langle \hat j_1(t) \hat j_1(0) \rangle + \langle \hat j_1(t) \hat j_2(0) \rangle + \langle \hat j_2(t) \hat j_1(0) \rangle + \langle \hat j_2(t) \hat j_2(0) \rangle \end{align}\end{split}\]

Returns:

\(n \times 4\) array for the decomposed correlation function defined as above where \(n\) is the number of time steps.

calc_mobility()

evolve_single_step(evolve_dt)

Evolve the mps for a single step with step size evolve_dt.

Returns:

new mps after the evolution

get_dump_dict()

Obtain calculated properties to dump in dict type.

Returns:

return a (ordered) dict to dump as npz

init_mps()

Returns:

initial mps of the system

process_mps(mps)

Process the newly evolved the mps. Primarily for the calculation of properties. Note that currently self.latest_mps has not been updated.

Parameters:

mps – The evolved new mps of the time step.

stop_evolve_criteria()

8.2. Charge Diffusion Dynamics Simulation

class renormalizer.transport.dynamics.ChargeDiffusionDynamics(model: ~renormalizer.model.model.HolsteinModel, temperature: ~renormalizer.utils.quantity.Quantity = <renormalizer.utils.quantity.Quantity object>, compress_config: ~typing.Optional[~renormalizer.utils.configs.CompressConfig] = None, evolve_config: ~typing.Optional[~renormalizer.utils.configs.EvolveConfig] = None, stop_at_edge: bool = True, init_electron=InitElectron.relaxed, rdm: bool = False, dump_dir: ~typing.Optional[str] = None, job_name: ~typing.Optional[str] = None)

Simulate charge diffusion dynamics by TD-DMRG. It is possible to obtain mobility from the simulation, but care must be taken to ensure that mean square displacement grows linearly with time.

Parameters:

• model (HolsteinModel) – system information. Currently only support HolsteinModel.

• temperature (Quantity) – simulation temperature. Default is zero temperature. For finite temperature charge dynamics, it is recommended to use thermal field dynamics and transform the Hamiltonian. See the documentation of SpectralFunctionZT for more details.

• compress_config (CompressConfig) – config when compressing MPS.

• evolve_config (EvolveConfig) – config when evolving MPS.

• stop_at_edge (bool) – whether stop when charge has diffused to the boundary of the system. Default is True.

• init_electron (InitElectron) – the method to prepare the initial state.

• rdm (bool) – whether calculate reduced density matrix and k-space representation for the electron. Default is False because usually the calculation is time consuming. Using scheme 4 might partly solve the problem.

• dump_dir (str) – the directory for logging and numerical result output. Also the directory from which to load previous thermal propagated initial state (if exists).

• job_name (str) – the name of the calculation job which determines the file name of the logging and numerical result output. For thermal propagated initial state input/output the file name is appended with "_impdm.npz".

energies

calculated energy of the states during the time evolution. Without dissipation or TD-Hartree the value should remain unchanged and to some extent can be used to measure the error during time evolution.

Type:

np.ndarray

r_square_array

calculated mean square displacement \(\langle \psi | \hat r^2 | \psi \rangle - \langle \psi | \hat r | \psi \rangle^2\) at each evolution time step.

Type:

np.ndarray

e_occupations_array

calculated electron occupations in real space on each site for each evolution time step.

Type:

np.ndarray

ph_occupations_array

calculated phonon occupations on each site for each evolution time step.

Type:

np.ndarray

reduced_density_matrices

calculated reduced density matrices of the electron for each evolution time step. Only available when rdm is set to True.

Type:

list

k_occupations_array

calculated electron occupations in momentum (k) space on each site for each evolution time step. Only available when rdm is set to True. The basis transformation is based on:

\[| k \rangle = \sum_j e^{-ijk} | j \rangle\]

where \(k\) starts from \(-\pi\) to \(\pi\) with interval \(2\pi/N\). \(N\) represents total number of electronic sites.

Type:

np.ndarray

coherent_length_array

coherent length \(L\) calculated for each evolution time step.

\[L = \sum_{ij, i \neq j} | \rho_{ij} |\]

where \(\rho\) is the density matrix of the electron. Naturally this is only available when rdm is set to True.

Type:

np.ndarray

create_electron(gs_mp)

create_electron_fc(gs_mp)

create_electron_relaxed(gs_mp)

evolve_single_step(evolve_dt)

Evolve the mps for a single step with step size evolve_dt.

Returns:

new mps after the evolution

get_dump_dict()

Obtain calculated properties to dump in dict type.

Returns:

return a (ordered) dict to dump as npz

init_mps()

Returns:

initial mps of the system

is_similar(other: ChargeDiffusionDynamics, rtol=0.001)

property mol_num

process_mps(mps)

Process the newly evolved the mps. Primarily for the calculation of properties. Note that currently self.latest_mps has not been updated.

Parameters:

mps – The evolved new mps of the time step.

stop_evolve_criteria()

class renormalizer.transport.dynamics.InitElectron(value)

Available methods to prepare initial state of charge diffusion

fc = 'franck-condon excitation'

relaxed = 'analytically relaxed phonon(s)'

renormalizer.transport.dynamics.calc_r_square(e_occupations)

8.3. One-Particle Spectral Function

class renormalizer.transport.spectral_function.SpectralFunctionZT(model: TI1DModel, compress_config: Optional[CompressConfig] = None, evolve_config: Optional[EvolveConfig] = None, dump_dir: Optional[str] = None, job_name: Optional[str] = None)

Calculate one-particle retarded Green’s function at zero temperature for translational invariant one-dimensional model:

\[iG_{ij}(t) = \langle 0 |c_i(t) c^\dagger_j | 0 \rangle\]

The value of the array is stored with key "G array" in the dumped output file. Because of the translational invariance, there are only two dimension in this array. The first dimension is \(t\) and the second dimension is \(|i-j|\). In k-space, \(iG_{ij}\) is transformed to:

\[iG_k(t) = \langle 0 |c_k(t) c^\dagger_k | 0 \rangle\]

and \(G_k(t)\) is saved with key "Gk array" in the dumped output file. Fourier transformation of \(G_k(t)\) results in the spectral function:

\[A(k, \omega) = -\frac{1}{\pi} \rm{Im} \int_0^\infty dt e^{i \omega t} G_k(t)\]

However, this value is not calculated directly in this class, since usually an broadening to \(G_k(t)\) is required and an optimal range of \(\omega\) can not be determined without additional knowledge.

For finite temperature spectral function, it is recommended to use thermal field dynamics and transform the Hamiltonian. An example is included in the test case. See J. Chem.Phys.145, 224101 (2016) and Phys. Rev. Lett.123, 090402 (2019) for details.

Parameters:

• model (TI1DModel) – system information. In principle should use TI1DModel. Using custom Model is possible however in this case the user should be responsible to ensure the translational invariance. It is recommended to set the number of the electronic DoFs to be an even number so that the \(k=\pi\) point is well defined.

• compress_config (CompressConfig) – config when compressing MPS.

• evolve_config (EvolveConfig) – config when carrying out real time propagation.

• dump_dir (str) – the directory for logging and numerical result output.

• job_name (str) – the name of the calculation job which determines the file name of the logging and numerical result output.

property G_array

\(G_{ij}(t)\) as a two dimensional array. The first dimension is \(t\) and the second dimension is \(|i-j|\).

Returns:

G_array – The Green’s function array

Return type:

np.ndarray

evolve_single_step(evolve_dt)

Evolve the mps for a single step with step size evolve_dt.

Returns:

new mps after the evolution

get_dump_dict()

Obtain calculated properties to dump in dict type.

Returns:

return a (ordered) dict to dump as npz

init_mps()

Returns:

initial mps of the system

process_mps(mps)

Process the newly evolved the mps. Primarily for the calculation of properties. Note that currently self.latest_mps has not been updated.

Parameters:

mps – The evolved new mps of the time step.