iDEA package¶

Submodules¶

iDEA.EXT1 module¶

Calculates the exact ground-state electron density and energy for a system of one electron through solving the Schroedinger equation. If the system is perturbed, the time-dependent electron density and current density are calculated. Excited states of the unperturbed system can also be calculated.

iDEA.EXT1.calculate_current_density(pm, density)[source]¶

Calculates the current density from the time-dependent electron density by solving the continuity equation.

\[\frac{\partial n}{\partial t} + \frac{\partial j}{\partial x} = 0\]

Parameters:	pm : object Parameters object density : array_like 2D array of the time-dependent density, indexed as density[time_index,space_index] returns array_like 2D array of the current density, indexed as current_density[time_index,space_index]

iDEA.EXT1.construct_A(pm, H)[source]¶

Constructs the sparse matrix A, used when solving Ax=b in the Crank-Nicholson propagation.

\[A = I + i \frac{\delta t}{2} H\]

Parameters:	pm : object Parameters object H : array_like 2D array containing the band elements of the Hamiltonian matrix, indexed as H[band,space_index] returns sparse_matrix The sparse matrix A

iDEA.EXT1.construct_K(pm)[source]¶

Stores the band elements of the kinetic energy matrix in lower form. The kinetic energy matrix is constructed using a three-point, five-point or seven-point stencil. This yields an NxN band matrix (where N is the number of grid points). For example with N=6 and a three-point stencil:

\[\begin{split}K = -\frac{1}{2} \frac{d^2}{dx^2}= -\frac{1}{2} \begin{pmatrix} -2 & 1 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 1 & -2 \end{pmatrix} \frac{1}{\delta x^2} = [\frac{1}{\delta x^2},-\frac{1}{2 \delta x^2}] \end{split}\]

Parameters:	pm : object Parameters object returns array_like 2D array containing the band elements of the kinetic energy matrix, indexed as K[band,space_index]

iDEA.EXT1.main(parameters)[source]¶

Calculates the ground-state of the system. If the system is perturbed, the time evolution of the perturbed system is then calculated.

Parameters:	parameters : object Parameters object returns object Results object

iDEA.EXT2 module¶

Calculates the exact ground-state electron density and energy for a system of two interacting electrons through solving the many-electron Schroedinger equation. If the system is perturbed, the time-dependent electron density and current density are calculated.

iDEA.EXT2.calculate_current_density(pm, density)[source]¶

Calculates the current density from the time-dependent electron density by solving the continuity equation.

\[\frac{\partial n}{\partial t} + \frac{\partial j}{\partial x} = 0\]

Parameters:	pm : object Parameters object density : array_like 2D array of the time-dependent density, indexed as density[time_index,space_index] returns array_like 2D array of the current density, indexed as current_density[time_index,space_index]

iDEA.EXT2.calculate_density(pm, wavefunction_2D)[source]¶

Calculates the electron density from the two-electron wavefunction.

\[n(x) = 2 \int_{-x_{\mathrm{max}}}^{x_{\mathrm{max}}} |\Psi(x,x_{2})|^{2} dx_{2}\]

Parameters:	pm : object Parameters object wavefunction : array_like 2D array of the wavefunction, indexed as wavefunction_2D[space_index_1,space_index_2] returns array_like 1D array of the density, indexed as density[space_index]

iDEA.EXT2.calculate_energy(pm, wavefunction_reduced, wavefunction_reduced_old)[source]¶

Calculates the energy of the system.

\[E = - \ln\bigg(\frac{|\Psi(x_{1},x_{2},\tau)|}{|\Psi(x_{1},x_{2},\tau - \delta \tau)|}\bigg) \frac{1}{\delta \tau}\]

Parameters:	pm : object Parameters object wavefunction_reduced : array_like 1D array of the reduced wavefunction at t, indexed as wavefunction_reduced[space_index_1_2] wavefunction_reduced_old : array_like 1D array of the reduced wavefunction at t-dt, indexed as wavefunction_reduced_old[space_index_1_2] returns float Energy of the system

iDEA.EXT2.construct_A_reduced(pm, reduction_matrix, expansion_matrix, td)[source]¶

Constructs the reduced form of the sparse matrix A.

\[\begin{split}\text{Imaginary time}: \ &A = I + \frac{\delta \tau}{2}H \\ \text{Real time}: \ &A = I + i\frac{\delta t}{2}H \\ \\ &A_{\mathrm{red}} = RAE \end{split}\]

where $R =$ reduction matrix and $E =$ expansion matrix

Parameters:

pm : object: Parameters object
reduction_matrix : sparse_matrix: Sparse matrix used to reduce the wavefunction (remove indistinct elements) by exploiting the exchange antisymmetry
expansion_matrix : sparse_matrix: Sparse matrix used to expand the reduced wavefunction (insert indistinct elements) to get back the full wavefunction
td : integer: 0 for imaginary time, 1 for real time
returns sparse_matrix: Reduced form of the sparse matrix A, used when solving the equation Ax=b

iDEA.EXT2.construct_antisymmetry_matrices(pm)[source]¶

Constructs the reduction and expansion matrices that are used to exploit the exchange antisymmetry of the wavefunction.

\[\Psi(x_{1},x_{2}) = -\Psi(x_{2},x_{1})\]

Parameters:	pm : object Parameters object returns sparse_matrix and sparse_matrix Reduction matrix used to reduce the wavefunction (remove indistinct elements). Expansion matrix used to expand the reduced wavefunction (insert indistinct elements) to get back the full wavefunction.

iDEA.EXT2.initial_wavefunction(pm)[source]¶

Generates the initial condition for the Crank-Nicholson imaginary time propagation.

\[\Psi(x_{1},x_{2}) = \frac{1}{\sqrt{2}}\big(\phi_{1}(x_{1})\phi_{2}(x_{2}) - \phi_{2}(x_{1})\phi_{1}(x_{2})\big)\]

Parameters:	pm : object Parameters object returns array_like 1D array of the reduced wavefunction, indexed as wavefunction_reduced[space_index_1_2]

iDEA.EXT2.main(parameters)[source]¶

Calculates the ground-state of the system. If the system is perturbed, the time evolution of the perturbed system is then calculated.

Parameters:	parameters : object Parameters object returns object Results object

iDEA.EXT2.non_approx(pm)[source]¶

Calculates the two lowest non-interacting eigenstates of the system. These can then be expressed in Slater determinant form as an approximation to the exact many-electron wavefunction.

\[\bigg(-\frac{1}{2} \frac{d^{2}}{dx^{2}} + V_{\mathrm{ext}}(x) \bigg) \phi_{j}(x) = \varepsilon_{j} \phi_{j}(x)\]

Parameters:	pm : object Parameters object returns array_like and array_like 1D array of the 1st non-interacting eigenstate, indexed as eigenstate_1[space_index]. 1D array of the 2nd non-interacting eigenstate, indexed as eigenstate_2[space_index].

iDEA.EXT2.qho_approx(pm, n)[source]¶

Calculates the nth eigenstate of the quantum harmonic oscillator, and shifts to ensure it is neither an odd nor an even function (necessary for the Gram-Schmidt algorithm).

\[ \begin{align}\begin{aligned}\bigg(-\frac{1}{2} \frac{d^{2}}{dx^{2}} + \frac{1}{2} \omega^{2} x^{2} \bigg) \phi_{n}(x) = \varepsilon_{n} \phi_{n}(x)\\\phi_{n}(x) = \frac{1}{\sqrt{2^{n}n!}} \bigg(\frac{\omega}{\pi}\bigg)^{1/4} e^{-\frac{\omega x^{2}}{2}} H_{n}\bigg(\sqrt{\omega}x \bigg)\end{aligned}\end{align} \]

Parameters:	pm : object Parameters object n : integer Principle quantum number returns array_like 1D array of the nth eigenstate, indexed as eigenstate[space_index]

iDEA.EXT2.solve_imaginary_time(pm, A_reduced, C_reduced, wavefunction_reduced, expansion_matrix)[source]¶

Propagates the initial wavefunction through imaginary time using the Crank-Nicholson method to find the ground-state of the system.

\[\begin{split}&\Big(I + \frac{\delta \tau}{2}H\Big) \Psi(x_{1},x_{2},\tau+\delta \tau) = \Big(I - \frac{\delta \tau}{2}H\Big) \Psi(x_{1},x_{2},\tau) \\ &\Psi(x_{1},x_{2},\tau) = \sum_{m}c_{m}e^{-\varepsilon_{m} \tau}\phi_{m} \implies \lim_{\tau \to \infty} \Psi(x_{1},x_{2},\tau) = \phi_{0}\end{split}\]

Parameters:

pm : object: Parameters object
A_reduced : sparse_matrix: Reduced form of the sparse matrix A, used when solving the equation Ax=b
C_reduced : sparse_matrix: Reduced form of the sparse matrix C, defined as C=-A+2I
wavefunction_reduced : array_like: 1D array of the reduced wavefunction, indexed as wavefunction_reduced[space_index_1_2]
expansion_matrix : sparse_matrix: Sparse matrix used to expand the reduced wavefunction (insert indistinct elements) to get back the full wavefunction
returns float and array_like: Energy of the ground-state system. 1D array of the ground-state wavefunction, indexed as wavefunction[space_index_1_2].

iDEA.EXT2.solve_real_time(pm, A_reduced, C_reduced, wavefunction, reduction_matrix, expansion_matrix)[source]¶

Propagates the ground-state wavefunction through real time using the Crank-Nicholson method to find the time-evolution of the perturbed system.

\[\Big(I + i\frac{\delta t}{2}H\Big) \Psi(x_{1},x_{2},t+\delta t) = \Big(I - i\frac{\delta t}{2}H\Big) \Psi(x_{1},x_{2},t) \]

Parameters:

pm : object: Parameters object
A_reduced : sparse_matrix: Reduced form of the sparse matrix A, used when solving the equation Ax=b
C_reduced : sparse_matrix: Reduced form of the sparse matrix C, defined as C=-A+2I
wavefunction : array_like: 1D array of the ground-state wavefunction, indexed as wavefunction[space_index_1_2]
reduction_matrix : sparse_matrix: Sparse matrix used to reduce the wavefunction (remove indistinct elements) by exploiting the exchange antisymmetry
expansion_matrix : sparse_matrix: Sparse matrix used to expand the reduced wavefunction (insert indistinct elements) to get back the full wavefunction
returns array_like and array_like: 2D array of the time-dependent density, indexed as density[time_index,space_index]. 2D array of the current density, indexed as current_density[time_index,space_index].

iDEA.EXT3 module¶

Calculates the exact ground-state electron density and energy for a system of three interacting electrons through solving the many-electron Schroedinger equation. If the system is perturbed, the time-dependent electron density and current density are calculated.

iDEA.EXT3.calculate_current_density(pm, density)[source]¶

Calculates the current density from the time-dependent electron density by solving the continuity equation.

\[\frac{\partial n}{\partial t} + \frac{\partial j}{\partial x} = 0\]

Parameters:	pm : object Parameters object density : array_like 2D array of the time-dependent density, indexed as density[time_index,space_index] returns array_like 2D array of the current density, indexed as current_density[time_index,space_index]

iDEA.EXT3.calculate_density(pm, wavefunction_3D)[source]¶

Calculates the electron density from the three-electron wavefunction.

\[n(x) = 3 \int_{-x_{\mathrm{max}}}^{x_{\mathrm{max}}} |\Psi(x,x_{2},x_{3})|^{2} dx_{2} \ dx_{3}\]

Parameters:	pm : object Parameters object wavefunction : array_like 3D array of the wavefunction, indexed as wavefunction_3D[space_index_1,space_index_2,space_index_3] returns array_like 1D array of the density, indexed as density[space_index]

iDEA.EXT3.calculate_energy(pm, wavefunction_reduced, wavefunction_reduced_old)[source]¶

Calculates the energy of the system.

\[E = - \ln\bigg(\frac{|\Psi(x_{1},x_{2},x_{3},\tau)|}{|\Psi(x_{1},x_{2},x_{3},\tau - \delta \tau)|}\bigg) \frac{1}{\delta \tau}\]

Parameters:	pm : object Parameters object wavefunction_reduced : array_like 1D array of the reduced wavefunction at t, indexed as wavefunction_reduced[space_index_1_2_3] wavefunction_reduced_old : array_like 1D array of the reduced wavefunction at t-dt, indexed as wavefunction_reduced_old[space_index_1_2_3] returns float Energy of the system

iDEA.EXT3.construct_A_reduced(pm, reduction_matrix, expansion_matrix, td)[source]¶

Constructs the reduced form of the sparse matrix A.

\[\begin{split}\text{Imaginary time}: \ &A = I + \frac{\delta \tau}{2}H \\ \text{Real time}: \ &A = I + i\frac{\delta t}{2}H \\ \\ &A_{\mathrm{red}} = RAE \end{split}\]

where $R =$ reduction matrix and $E =$ expansion matrix

Parameters:

pm : object: Parameters object
reduction_matrix : sparse_matrix: Sparse matrix used to reduce the wavefunction (remove indistinct elements) by exploiting the exchange antisymmetry
expansion_matrix : sparse_matrix: Sparse matrix used to expand the reduced wavefunction (insert indistinct elements) to get back the full wavefunction
td : integer: 0 for imaginary time, 1 for real time
returns sparse_matrix: Reduced form of the sparse matrix A, used when solving the equation Ax=b

iDEA.EXT3.construct_antisymmetry_matrices(pm)[source]¶

Constructs the reduction and expansion matrices that are used to exploit the exchange antisymmetry of the wavefunction.

\[\Psi(x_{1},x_{2},x_{3}) = -\Psi(x_{2},x_{1},x_{3}) = \Psi(x_{2},x_{3},x_{1})\]

Parameters:	pm : object Parameters object returns sparse_matrix and sparse_matrix Reduction matrix used to reduce the wavefunction (remove indistinct elements). Expansion matrix used to expand the reduced wavefunction (insert indistinct elements) to get back the full wavefunction.

iDEA.EXT3.initial_wavefunction(pm)[source]¶

Generates the initial wavefunction for the Crank-Nicholson imaginary time propagation.

\[\Psi(x_{1},x_{2},x_{3}) = \frac{1}{\sqrt{6}}\big(\phi_{1}(x_{1}) \phi_{2}(x_{2})\phi_{3}(x_{3}) - \phi_{1}(x_{1})\phi_{3}(x_{2}) \phi_{2}(x_{3}) + \phi_{3}(x_{1})\phi_{1}(x_{2})\phi_{2}(x_{3}) - \phi_{3}(x_{1})\phi_{2}(x_{2})\phi_{1}(x_{3}) + \phi_{2}(x_{1}) \phi_{3}(x_{2})\phi_{1}(x_{3}) - \phi_{2}(x_{1})\phi_{1}(x_{2}) \phi_{3}(x_{3}) \big)\]

Parameters:	pm : object Parameters object returns array_like 1D array of the reduced wavefunction, indexed as wavefunction_reduced[space_index_1_2_3]

iDEA.EXT3.main(parameters)[source]¶

Calculates the ground-state of the system. If the system is perturbed, the time evolution of the perturbed system is then calculated.

Parameters:	parameters : object Parameters object returns object Results object

iDEA.EXT3.non_approx(pm)[source]¶

Calculates the three lowest non-interacting eigenstates of the system. These can then be expressed in Slater determinant form as an approximation to the exact many-electron wavefunction.

\[\bigg(-\frac{1}{2} \frac{d^{2}}{dx^{2}} + V_{\mathrm{ext}}(x) \bigg) \phi_{j}(x) = \varepsilon_{j} \phi_{j}(x)\]

Parameters:	pm : object Parameters object returns array_like and array_like and array_like 1D array of the 1st non-interacting eigenstate, indexed as eigenstate_1[space_index]. 1D array of the 2nd non-interacting eigenstate, indexed as eigenstate_2[space_index]. 1D array of the 3rd non-interacting eigenstate, indexed as eigenstate_3[space_index].

iDEA.EXT3.qho_approx(pm, n)[source]¶

Calculates the nth eigenstate of the quantum harmonic oscillator, and shifts to ensure it is neither an odd nor an even function (necessary for the Gram-Schmidt algorithm).

\[ \begin{align}\begin{aligned}\bigg(-\frac{1}{2} \frac{d^{2}}{dx^{2}} + \frac{1}{2} \omega^{2} x^{2} \bigg) \phi_{n}(x) = \varepsilon_{n} \phi_{n}(x)\\\phi_{n}(x) = \frac{1}{\sqrt{2^{n}n!}} \bigg(\frac{\omega}{\pi}\bigg)^{1/4} e^{-\frac{\omega x^{2}}{2}} H_{n}\bigg(\sqrt{\omega}x \bigg)\end{aligned}\end{align} \]

Parameters:	pm : object Parameters object n : integer Principle quantum number returns array_like 1D array of the nth eigenstate, indexed as eigenstate[space_index]

iDEA.EXT3.solve_imaginary_time(pm, A_reduced, C_reduced, wavefunction_reduced, expansion_matrix)[source]¶

Propagates the initial wavefunction through imaginary time using the Crank-Nicholson method to find the ground-state of the system.

\[\begin{split}&\Big(I + \frac{\delta \tau}{2}H\Big) \Psi(x_{1},x_{2},x_{3},\tau+\delta \tau) = \Big(I - \frac{\delta \tau}{2}H\Big) \Psi(x_{1},x_{2},x_{3},\tau) \\ &\Psi(x_{1},x_{2},x_{3},\tau) = \sum_{m}c_{m}e^{-\varepsilon_{m} \tau}\phi_{m} \implies \lim_{\tau \to \infty} \Psi(x_{1},x_{2},x_{3},\tau) = \phi_{0}\end{split}\]

Parameters:

pm : object: Parameters object
A_reduced : sparse_matrix: Reduced form of the sparse matrix A, used when solving the equation Ax=b
C_reduced : sparse_matrix: Reduced form of the sparse matrix C, defined as C=-A+2I
wavefunction_reduced : array_like: 1D array of the reduced wavefunction, indexed as wavefunction_reduced[space_index_1_2_3]
expansion_matrix : sparse_matrix: Sparse matrix used to expand the reduced wavefunction (insert indistinct elements) to get back the full wavefunction
returns float and array_like: Energy of the ground-state system. 1D array of the ground-state wavefunction, indexed as wavefunction[space_index_1_2_3].

iDEA.EXT3.solve_real_time(pm, A_reduced, C_reduced, wavefunction, reduction_matrix, expansion_matrix)[source]¶

Propagates the ground-state wavefunction through real time using the Crank-Nicholson method to find the time-evolution of the perturbed system.

\[\Big(I + i\frac{\delta t}{2}H\Big) \Psi(x_{1},x_{2},x_{3},t+\delta t) = \Big(I - i\frac{\delta t}{2}H\Big) \Psi(x_{1},x_{2},x_{3},t) \]

Parameters:

pm : object: Parameters object
A_reduced : sparse_matrix: Reduced form of the sparse matrix A, used when solving the equation Ax=b
C_reduced : sparse_matrix: Reduced form of the sparse matrix C, defined as C=-A+2I
wavefunction : array_like: 1D array of the ground-state wavefunction, indexed as wavefunction[space_index_1_2_3]
reduction_matrix : sparse_matrix: Sparse matrix used to reduce the wavefunction (remove indistinct elements) by exploiting the exchange antisymmetry
expansion_matrix : sparse_matrix: Sparse matrix used to expand the reduced wavefunction (insert indistinct elements) to get back the full wavefunction
returns array_like and array_like: 2D array of the time-dependent density, indexed as density[time_index,space_index]. 2D array of the current density, indexed as current_density[time_index,space_index].

iDEA.EXT_cython module¶

Contains the cython modules that are called within EXT2 and EXT3. Cython is used for operations that are very expensive to do in Python, and performance speeds are close to C.

iDEA.EXT_cython.continuity_eqn()¶

Calculates the electron current density of the system at a particular time step by solving the continuity equation.

\[\frac{\partial n}{\partial t} + \frac{\partial j}{\partial x} = 0\]

Parameters:	pm : object Parameters object density_new : array_like 1D array of the electron density at time t density_old : array_like 1D array of the electron density at time t-dt returns array_like 1D array of the electron current density at time t

iDEA.EXT_cython.expansion_three()¶

Calculates the coordinates and data of the non-zero elements of the three-electron expansion matrix that is used to exploit the exchange antisymmetry of the three-electron wavefunction.

Parameters:

pm : object: Parameters object
coo_1 : array_like: 1D COOrdinate holding array for the non-zero elements of the expansion matrix
coo_2 : array_like: 1D COOrdinate holding array for the non-zero elements of the expansion matrix
coo_data : array_like: 1D array of the non-zero elements of the expansion matrix
returns array_like and array_like and array_like: Populated 1D COOrdinate holding arrays and 1D data holding array for the non-zero elements of the expansion matrix

iDEA.EXT_cython.expansion_two()¶

Calculates the coordinates and data of the non-zero elements of the two-electron expansion matrix that is used to exploit the exchange antisymmetry of the two-electron wavefunction.

Parameters:

pm : object: Parameters object
coo_1 : array_like: 1D COOrdinate holding array for the non-zero elements of the expansion matrix
coo_2 : array_like: 1D COOrdinate holding array for the non-zero elements of the expansion matrix
coo_data : array_like: 1D array of the non-zero elements of the expansion matrix
returns array_like and array_like and array_like: Populated 1D COOrdinate holding arrays and 1D data holding array for the non-zero elements of the expansion matrix

iDEA.EXT_cython.hamiltonian_three()¶

Calculates the coordinates and data of the non-zero elements of the three-electron Hamiltonian matrix.

\[\hat{H}=\sum_{i=1}^{3} \hat{K}_{i} + \sum_{i=1}^{3} \hat{V}_{\mathrm{ext}}(x_{i}) + \sum_{i=1}^{3}\sum_{j > i}^{3}\hat{V}_{\mathrm{int}}(x_{i}, x_{j})\]

Parameters:

pm : object: Parameters object
coo_1 : array_like: 1D COOrdinate holding array for the non-zero elements of the Hamiltonian matrix
coo_2 : array_like: 1D COOrdinate holding array for the non-zero elements of the Hamiltonian matrix
coo_data : array_like: 1D array of the non-zero elements of the Hamiltonian matrix
td : integer: 0 for imaginary time, 1 for real time
returns array_like and array_like and array_like: Populated 1D COOrdinate holding arrays and 1D data holding array for the non-zero elements of the Hamiltonian matrix

iDEA.EXT_cython.hamiltonian_two()¶

Calculates the coordinates and data of the non-zero elements of the two-electron Hamiltonian matrix.

\[\hat{H}=\sum_{i=1}^{2} \hat{K}_{i} + \sum_{i=1}^{2} \hat{V}_{\mathrm{ext}}(x_{i}) + \sum_{i=1}^{2}\sum_{j > i}^{2}\hat{V}_{\mathrm{int}}(x_{i}, x_{j})\]

Parameters:

pm : object: Parameters object
coo_1 : array_like: 1D COOrdinate holding array for the non-zero elements of the Hamiltonian matrix
coo_2 : array_like: 1D COOrdinate holding array for the non-zero elements of the Hamiltonian matrix
coo_data : array_like: 1D array of the non-zero elements of the Hamiltonian matrix
td : integer: 0 for imaginary time, 1 for real time
returns array_like and array_like and array_like: Populated 1D COOrdinate holding arrays and 1D data holding array for the non-zero elements of the Hamiltonian matrix

iDEA.EXT_cython.imag_pot_three()¶

Calculates the imaginary component of the perturbing potential to be added to the main diagonal of the three-electron Hamiltonian matrix if imaginary potentials are used.

Parameters:	pm : object Parameters object returns array_like 1D array of the perturbing potential

iDEA.EXT_cython.imag_pot_two()¶

Calculates the imaginary component of the perturbing potential to be added to the main diagonal of the two-electron Hamiltonian matrix if imaginary potentials are used.

Parameters:	pm : object Parameters object returns array_like 1D array of the perturbing potential

iDEA.EXT_cython.reduction_three()¶

Calculates the coordinates and data of the non-zero elements of the three-electron reduction matrix that is used to exploit the exchange antisymmetry of the three-electron wavefunction.

Parameters:	pm : object Parameters object coo_1 : array_like 1D COOrdinate holding array for the non-zero elements of the reduction matrix coo_2 : array_like 1D COOrdinate holding array for the non-zero elements of the reduction matrix returns array_like and array_like Populated 1D COOrdinate holding arrays for the non-zero elements of the reduction matrix

iDEA.EXT_cython.reduction_two()¶

Calculates the coordinates and data of the non-zero elements of the two-electron reduction matrix that is used to exploit the exchange antisymmetry of the two-electron wavefunction.

Parameters:	pm : object Parameters object coo_1 : array_like 1D COOrdinate holding array for the non-zero elements of the reduction matrix coo_2 : array_like 1D COOrdinate holding array for the non-zero elements of the reduction matrix returns array_like and array_like Populated 1D COOrdinate holding arrays for the non-zero elements of the reduction matrix

iDEA.EXT_cython.single_index_three()¶

Takes every permutation of the three electron indices and creates a single unique index.

Parameters:	j : integer First electron index k : integer Second electron index l : integer Third electron index grid : integer Number of spatial grid points in the system returns integer Single unique index

iDEA.EXT_cython.single_index_two()¶

Takes every permutation of the two electron indices and creates a single unique index.

Parameters:	j : integer First electron index k : integer Second electron index grid : integer Number of spatial grid points in the system returns integer Single unique index

iDEA.EXT_cython.wavefunction_three()¶

Constructs the initial three-electron wavefunction in reduced form from three single-electron eigenstates.

\[\Psi(x_{1},x_{2},x_{3}) = \frac{1}{\sqrt{6}}\big(\phi_{1}(x_{1}) \phi_{2}(x_{2})\phi_{3}(x_{3}) - \phi_{1}(x_{1})\phi_{3}(x_{2}) \phi_{2}(x_{3}) + \phi_{3}(x_{1})\phi_{1}(x_{2})\phi_{2}(x_{3}) - \phi_{3}(x_{1})\phi_{2}(x_{2})\phi_{1}(x_{3}) + \phi_{2}(x_{1}) \phi_{3}(x_{2})\phi_{1}(x_{3}) - \phi_{2}(x_{1})\phi_{1}(x_{2}) \phi_{3}(x_{3}) \big)\]

Parameters:

pm : object: Parameters object
eigenstate_1 : array_like: 1D array of the 1st single-electron eigenstate, indexed as eigenstate_1[space_index]
eigenstate_2 : array_like: 1D array of the 2nd single-electron eigenstate, indexed as eigenstate_2[space_index]
eigenstate_3 : array_like: 1D array of the 3rd single-electron eigenstate, indexed as eigenstate_3[space_index]
returns array_like: 1D array of the reduced wavefunction, indexed as wavefunction_reduced[space_index_1_2_3]

iDEA.EXT_cython.wavefunction_two()¶

Constructs the two-electron initial wavefunction in reduced form from two single-electron eigenstates.

\[\Psi(x_{1},x_{2}) = \frac{1}{\sqrt{2}}\big(\phi_{1}(x_{1})\phi_{2}(x_{2}) - \phi_{2}(x_{1})\phi_{1}(x_{2})\big)\]

Parameters:	pm : object Parameters object eigenstate_1 : array_like 1D array of the 1st single-electron eigenstate, indexed as eigenstate_1[space_index] eigenstate_2 : array_like 1D array of the 2nd single-electron eigenstate, indexed as eigenstate_2[space_index] returns array_like 1D array of the reduced wavefunction, indexed as wavefunction_reduced[space_index_1_2]

iDEA.HF module¶

Computes ground-state and time-dependent charge density in the Hartree-Fock approximation. The code outputs the ground-state charge density, the energy of the system and the Hartree-Fock orbitals.

iDEA.HF.calculate_current_density(pm, density)[source]¶

Calculates the current density from the time-dependent (and ground-state) electron density by solving the continuity equation.

\[\frac{\partial n}{\partial t} + \nabla \cdot j = 0\]

Parameters:	pm : object Parameters object density : array_like 2D array of the time-dependent density, indexed as density[time_index,space_index] returns array_like 2D array of the current density, indexed as current_density[time_index,space_index]

iDEA.HF.crank_nicolson_step(pm, waves, H)[source]¶

Solves Crank Nicolson Equation

\[\left(\mathbb{1} + i\frac{dt}{2} H\right) \Psi(x,t+dt) = \left(\mathbb{1} - i \frac{dt}{2} H\right) \Psi(x,t)\]

for $\Psi(x,t+dt)$.

Parameters:	pm : object Parameters object total_td_density : array_like Time dependent density of the system indexed as total_td_density[time_index][space_index] returns array_like Time dependent current density indexed as current_density[time_index][space_index]

iDEA.HF.electron_density(pm, orbitals)[source]¶

Compute density for given orbitals

Parameters:	pm : object Parameters object orbitals: array_like Array of properly normalised orbitals
Returns:	n: array_like electron density

iDEA.HF.fock(pm, eigf)[source]¶

Constructs Fock operator from a set of orbitals

\[\Sigma_{\text{x}}(x,y) = - \sum_{j=1}^{N} \varphi_{j}^{*}(y) \varphi_{j}(x) u(x,y)\]

where u(x,y) denotes the appropriate Coulomb interaction.

Parameters:	pm : object Parameters object eigf : array_like Eigenfunction orbitals
Returns:	F: array_like Fock matrix

iDEA.HF.groundstate(pm, H)[source]¶

Diagonalises Hamiltonian H

\[\hat{H}\varphi_{j} = \varepsilon_{j}\varphi_{j}\]

Parameters:	pm: object Parameters object H: array_like Hamiltonian matrix (band form)
Returns:	n: array_like density eigf: array_like normalised orbitals, index as eigf[space_index,orbital_number] eigv: array_like orbital energies

iDEA.HF.hamiltonian(pm, wfs, perturb=False)[source]¶

Compute HF Hamiltonian

Computes HF Hamiltonian from a given set of single-particle states

Parameters:	pm : object Parameters object wfs array_like single-particle states perturb: bool If True, add perturbation to external potential (for time-dep. runs) returns array_like Hamiltonian matrix

iDEA.HF.hartree(pm, density)[source]¶

Computes Hartree potential for a given density

\[v_{\text{H}}(x) = \int \! n(y)u(x,y) \, \mathrm{d}y\]

Parameters:	pm : object Parameters object density : array_like given density returns array_like Hartree potential

iDEA.HF.main(parameters)[source]¶

Performs Hartree-fock calculation

Parameters:	parameters : object Parameters object returns object Results object

iDEA.HF.total_energy(pm, eigf, eigv)[source]¶

Calculates total energy of Hartree-Fock wave function

Parameters:	pm : array_like external potential eigf : array_like eigenfunctions eigv : array_like eigenvalues returns float

iDEA.HYB module¶

Computes ground-state and time-dependent charge density in the Hybrid Hartree-Fock-LDA approximation.

The code outputs the ground-state charge density, the energy of the system and the single-quasiparticle orbitals.

iDEA.HYB.calc_with_alpha(pm, alpha, occupations)[source]¶

Calculate with given alpha.

Perform hybrid calculation with given alpha.

Parameters:	alpha float HF-LDA mixing parameter (1 = all HF, 0 = all LDA) occupations: array_like orbital occupations returns density, eigf, eigv, E Hybrid Density, orbitals and total energy

iDEA.HYB.fractional_run(pm, results, occupations, fractions)[source]¶

iDEA.HYB.hamiltonian(pm, eigf, density, alpha, occupations, perturb=False)[source]¶

Compute HF Hamiltonian.

Computes HYB Hamiltonian from a given set of single-particle states.

\[H_{\alpha}(x,y) = \delta(x-y)\hat{T} + \delta(x-y)v_{\text{ext}}(y) + \delta(x-y)v_{\text{H}}(y) + \alpha\Sigma_{\text{x}}(x,y) + (1-\alpha)\delta(x-y)v_{\text{xc}}^{\text{LDA}}(y)\]

Parameters:	eigf array_like single-particle states density array_like electron density alpha float HF-LDA mixing parameter (1 = all HF, 0 = all LDA) occupations: array_like orbital occupations perturb: bool If True, add perturbation to external potential (for time-dep. runs) returns array_like Hamiltonian matrix

iDEA.HYB.main(parameters)[source]¶

Performs Hybrid calculation.

Parameters:	parameters : object Parameters object returns object Results object

iDEA.HYB.n_minus_one_run(pm, results, alphas, occupations)[source]¶

Calculate for $N-1$ electron run.

Calculate total energy and LUMO eigenvalue of $N-1$ electron system.

Parameters:	results Results Object object to add results to alphas: array_like range of alphas to use occupations: array_like orbital occupations

iDEA.HYB.n_run(pm, results, alphas, occupations)[source]¶

Calculate for $N$ electron run.

Calculate total energy and HOMO eigenvalue of N electron system.

Parameters:	results Results Object object to add results to alphas: array_like range of alphas to use occupations: array_like orbital occupations

iDEA.HYB.optimal_alpha(pm, results, alphas, occupations)[source]¶

Calculate optimal alpha.

Calculate over range of alphas to determine optimal alpha.

Parameters:	results Results Object object to add results to alphas: array_like range of alphas to use occupations: array_like orbital occupations

iDEA.HYB.save_results(pm, results, density, E, eigf, eigv, alpha)[source]¶: Saves hybrid results to outputs directory

iDEA.LDA module¶

Uses the [adiabatic] local density approximation ([A]LDA) to calculate the [time-dependent] electron density [and current] for a system of N electrons.

Computes approximations to V_KS, V_H, V_xc using the LDA self-consistently. For ground state calculations the code outputs the LDA orbitals and energies of the system, the ground-state charge density and Kohn-Sham potential. For time dependent calculations the code also outputs the time-dependent charge and current densities and the time-dependent Kohn-Sham potential.

Note: Uses the LDAs developed in [Entwistle2018] from finite slab systems and the HEG, in one dimension.

iDEA.LDA.DXC(pm, n)[source]¶

Calculates the derivative of the exchange-correlation potential, necessary for the RPA preconditioner.

Parameters:	pm : object Parameters object n : array_like 1D array of the electron density, indexed as n[space_index] returns array_like 1D array of the derivative of the exchange-correlation potential, indexed as D_xc[space_index]

iDEA.LDA.banded_to_full(pm, H)[source]¶

Converts the Hamiltonian matrix in band form to the full matrix.

Parameters:	pm : object Parameters object H : array_like 2D array of the Hamiltonian matrix in band form, indexed as H[band,space_index]
Returns:	H_full : array_like 2D array of the Hamiltonian matrix in full form, indexed as H_full[space_index,space_index]

iDEA.LDA.calculate_current_density(pm, density)[source]¶

Calculates the current density from the time-dependent electron density by solving the continuity equation.

\[\frac{\partial n}{\partial t} + \frac{\partial j}{\partial x} = 0\]

Parameters:	pm : object Parameters object density : array_like 2D array of the time-dependent density, indexed as density[time_index,space_index] returns array_like 2D array of the current density, indexed as current_density[time_index,space_index]

iDEA.LDA.crank_nicolson_step(pm, orbitals, H_full)[source]¶

Solves Crank Nicolson Equation

\[\left(I + i\frac{dt}{2} H\right) \Psi(x,t+dt) = \left(I - i \frac{dt}{2} H\right) \Psi(x,t)\]

Parameters:	pm : object Parameters object orbitals : array_like 2D array of the Kohn-Sham orbitals, index as orbitals[space_index,orbital_number] H_full : array_like 2D array of the Hamiltonian matrix in full form, indexed as H_full[space_index,space_index] returns

iDEA.LDA.electron_density(pm, orbitals)[source]¶

Calculates the electron density from the set of orbitals.

\[n(x) = \sum_{j=1}^{N}|\phi_{j}(x)|^{2}\]

Parameters:	pm : object Parameters object orbitals : array_like 2D array of the Kohn-Sham orbitals, index as orbitals[space_index,orbital_number]
Returns:	density : array_like 1D array of the electron density, indexed as density[space_index]

iDEA.LDA.groundstate(pm, H)[source]¶

Calculates the ground-state of the system for a given potential.

\[\hat{H} \phi_{j} = \varepsilon_{j} \phi_{j}\]

Parameters:	pm : object Parameters object H : array_like 2D array of the Hamiltonian matrix in band form, indexed as H[band,space_index]
Returns:	density : array_like 1D array of the electron density, indexed as density[space_index] orbitals : array_like 2D array of the Kohn-Sham orbitals, index as orbitals[space_index,orbital_number] eigenvalues : array_like 1D array of the Kohn-Sham eigenvalues, indexed as eigenvalues[orbital_number]

iDEA.LDA.hamiltonian(pm, v_ks=None, orbitals=None, perturbation=False)[source]¶

Constructs the Hamiltonian matrix in band form for a given Kohn-Sham potential.

\[\hat{H} = \hat{K} + \hat{V}_{\mathrm{KS}}\]

Parameters:	pm : object Parameters object v_ks : array_like 1D array of the Kohn-Sham potential, indexed as v_ks[space_index] orbitals : array_like 2D array of the Kohn-Sham orbitals, index as orbitals[space_index,orbital_number] perturbation: bool True: Perturbed external potential False: Unperturbed external potential
Returns:	H : array_like 2D array of the Hamiltonian matrix in band form, indexed as H[band,space_index]

iDEA.LDA.hartree_energy(pm, v_h, density)[source]¶

Calculates the Hartree energy of the ground-state system.

\[E_{\mathrm{H}}[n] = \frac{1}{2} \int \int U(x,x') n(x) n(x') dx dx' = \frac{1}{2} \int V_{\mathrm{H}}(x) n(x) dx\]

Parameters:	pm : object Parameters object v_h : array_like 1D array of the ground-state Hartree potential, indexed as v_h[space_index] density : array_like 1D array of the ground-state electron density, indexed as density[space_index] returns float The Hartree energy of the ground-state system

iDEA.LDA.hartree_potential(pm, density)[source]¶

Calculates the Hartree potential for a given electron density.

\[V_{\mathrm{H}}(x) = \int U(x,x') n(x') dx'\]

Parameters:	pm : object Parameters object density : array_like 1D array of the electron density, indexed as density[space_index] returns array_like 1D array of the Hartree potential, indexed as v_h[space_index]

iDEA.LDA.kinetic(pm)[source]¶

Stores the band elements of the kinetic energy matrix in lower form. The kinetic energy matrix is constructed using a three-point, five-point or seven-point stencil. This yields an NxN band matrix (where N is the number of grid points). For example with N=6 and a three-point stencil:

\[\begin{split}K = -\frac{1}{2} \frac{d^2}{dx^2}= -\frac{1}{2} \begin{pmatrix} -2 & 1 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 1 & -2 \end{pmatrix} \frac{1}{\delta x^2} = [\frac{1}{\delta x^2},-\frac{1}{2 \delta x^2}] \end{split}\]

Parameters:	pm : object Parameters object returns array_like 2D array containing the band elements of the kinetic energy matrix, indexed as K[band,space_index]

iDEA.LDA.kinetic_energy(pm, orbitals)[source]¶

Calculates the kinetic energy from the Kohn-Sham orbitals.

\[T_{s}[n] = \sum_{j=1}^{N} \langle \phi_{j} | K | \phi_{j} \rangle\]

Parameters:	pm : object Parameters object orbitals : array_like 2D array of the Kohn-Sham orbitals, index as orbitals[space_index,orbital_number]

iDEA.LDA.ks_potential(pm, density, perturbation=False)[source]¶

Calculates the Kohn-Sham potential from the electron density.

\[V_{\mathrm{KS}} = V_{\mathrm{ext}} + V_{\mathrm{H}} + V_{\mathrm{xc}}\]

Parameters:	pm : object Parameters object density : array_like 1D array of the electron density, indexed as density[space_index] perturbation: bool True: Perturbed external potential False: Unperturbed external potential
Returns:	v_ks : array_like 1D array of the Kohn-Sham potential, indexed as v_ks[space_index]

iDEA.LDA.main(parameters)[source]¶

Performs LDA calculation

Parameters:	parameters : object Parameters object returns object Results object

iDEA.LDA.total_energy_eigf(pm, orbitals, density=None, v_h=None)[source]¶

Calculates the total energy from the Kohn-Sham orbitals.

\[E[n] = \sum_{j=1}^{N} \langle \phi_{j} | K | \phi_{j} \rangle + E_H[n] + E_{xc}[n] + \int n(x) V_{\mathrm{ext}}(x)dx\]

Parameters:	pm : object Parameters object orbitals : array_like 2D array of the Kohn-Sham orbitals, index as orbitals[space_index,orbital_number] density : array_like 1D array of the electron density, indexed as density[space_index] v_h : array_like 1D array of the Hartree potential, indexed as v_h[space_index] returns float Total energy

iDEA.LDA.total_energy_eigv(pm, eigenvalues, orbitals=None, density=None, v_h=None, v_xc=None)[source]¶

Calculates the total energy from the Kohn-Sham eigenvalues.

\[E[n] = \sum_{j=1}^{N} \varepsilon_j + E_{xc}[n] - E_H[n] - \int n(x) V_{xc}(x)dx\]

Parameters:

pm : object: Parameters object
eigenvalues : array_like: 1D array of the Kohn-Sham eigenvalues, indexed as eigenvalues[orbital_number]
orbitals : array_like: 2D array of the Kohn-Sham orbitals, index as orbitals[space_index,orbital_number]
density : array_like: 1D array of the electron density, indexed as density[space_index]
v_h : array_like: 1D array of the Hartree potential, indexed as v_h[space_index]
v_xc : array_like: 1D array of the exchange-correlation potential, indexed as v_xc[space_index]
returns float: Total energy

iDEA.LDA.xc_energy(pm, n, separate=False)[source]¶

LDA approximation for the exchange-correlation energy. Uses the LDAs developed in [Entwistle et al. 2018] from finite slab systems and the HEG.

\[E_{\mathrm{xc}}^{\mathrm{LDA}}[n] = \int \varepsilon_{\mathrm{xc}}(n) n(x) dx\]

Parameters:	pm : object Parameters object n : array_like 1D array of the electron density, indexed as n[space_index] separate: bool True: Split the HEG exchange-correlation energy into separate exchange and correlation terms False: Just return the exchange-correlation energy returns float Exchange-correlation energy

iDEA.LDA.xc_potential(pm, n, separate=False)[source]¶

LDA approximation for the exchange-correlation potential. Uses the LDAs developed in [Entwistle et al. 2018] from finite slab systems and the HEG.

\[V_{\mathrm{xc}}^{\mathrm{LDA}}(x) = \frac{\delta E_{\mathrm{xc}}^{\mathrm{LDA}}[n]}{\delta n(x)} = \varepsilon_{\mathrm{xc}}(n(x)) + n(x)\frac{d\varepsilon_{\mathrm{xc}}}{dn} \bigg|_{n(x)}\]

Parameters:	pm : object Parameters object n : array_like 1D array of the electron density, indexed as n[space_index] separate: bool True: Split the HEG exchange-correlation potential into separate exchange and correlation terms False: Just return the exchange-correlation potential returns array_like 1D array of the exchange-correlation potential, indexed as v_xc[space_index]

iDEA.LDA_parameters module¶

These are the parameters for the LDAs developed in [Entwistle2018] from finite slab systems and the HEG.

iDEA.NON module¶

Calculates the ground-state electron density and energy for a system of N non-interacting electrons through solving the Schroedinger equation. If the system is perturbed, the time-dependent electron density and current density are calculated.

iDEA.NON.calculate_current_density(pm, density)[source]¶

Calculates the current density from the time-dependent electron density by solving the continuity equation.

\[\frac{\partial n}{\partial t} + \frac{\partial j}{\partial x} = 0\]

Parameters:	pm : object Parameters object density : array_like 2D array of the time-dependent density, indexed as density[time_index,space_index] returns array_like 2D array of the current density, indexed as current_density[time_index,space_index]

iDEA.NON.construct_A(pm, H)[source]¶

Constructs the sparse matrix A, used when solving Ax=b in the Crank-Nicholson propagation.

\[A = I + i \frac{\delta t}{2} H\]

Parameters:	pm : object Parameters object H : array_like 2D array containing the band elements of the Hamiltonian matrix, indexed as H[band,space_index] returns sparse_matrix The sparse matrix A

iDEA.NON.construct_K(pm)[source]¶

Stores the band elements of the kinetic energy matrix in lower form. The kinetic energy matrix is constructed using a three-point, five-point or seven-point stencil. This yields an NxN band matrix (where N is the number of grid points). For example with N=6 and a three-point stencil:

\[\begin{split}K = -\frac{1}{2} \frac{d^2}{dx^2}= -\frac{1}{2} \begin{pmatrix} -2 & 1 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 1 & -2 \end{pmatrix} \frac{1}{\delta x^2} = [\frac{1}{\delta x^2},-\frac{1}{2 \delta x^2}] \end{split}\]

Parameters:	pm : object Parameters object returns array_like 2D array containing the band elements of the kinetic energy matrix, indexed as K[band,space_index]

iDEA.NON.main(parameters)[source]¶

Calculates the ground-state of the system. If the system is perturbed, the time evolution of the perturbed system is then calculated.

Parameters:	parameters : object Parameters object returns object Results object

iDEA.RE module¶

Calculates the exact Kohn-Sham potential and exchange-correlation potential for a given electron density using the reverse-engineering algorithm. This works for both a ground-state and time-dependent density.

iDEA.RE.calculate_current_density(pm, density_ks)[source]¶

Calculates the Kohn-Sham electron current density, at time t+dt, from the time-dependent Kohn-Sham electron density by solving the continuity equation.

\[\frac{\partial n_{\mathrm{KS}}}{\partial t} + \nabla \cdot j_{\mathrm{KS}} = 0\]

Parameters:	pm : object Parameters object density_ks : array_like 2D array of the time-dependent Kohn-Sham electron density, at time t and t+dt, indexed as density_ks[time_index, space_index] returns array_like 1D array of the Kohn-Sham electron current density, at time t+dt, indexed as current_density_ks[space_index]

iDEA.RE.calculate_ground_state(pm, approx, density_approx, v_ext, K)[source]¶

Calculates the exact ground-state Kohn-Sham potential by solving the ground-state Kohn-Sham equations and iteratively correcting v_ks. The exact ground-state Kohn-Sham eigenfunctions, eigenenergies and electron density are then calculated.

\[V_{\mathrm{KS}}(x) \rightarrow V_{\mathrm{KS}}(x) + \mu[n_{\mathrm{KS}}^{p}(x)-n_{\mathrm{approx}}^{p}(x)]\]

Parameters:

pm : object: Parameters object
approx : string: The approximation used to calculate the electron density
density_approx : array_like: 1D array of the ground-state electron density from the approximation, indexed as density_approx[space_index]
v_ext : array_like: 1D array of the unperturbed external potential, indexed as v_ext[space_index]
K : array_like: 2D array of the kinetic energy matrix, index as K[band,space_index]
returns array_like and array_like and array_like and array_like and Boolean: 1D array of the ground-state Kohn-Sham potential, indexed as v_ks[space_index]. 1D array of the ground-state Kohn-Sham electron density, indexed as density_ks[space_index]. 2D array of the ground-state Kohn-Sham eigenfunctions, indexed as wavefunctions_ks[space_index,eigenfunction]. 1D array containing the ground-state Kohn-Sham eigenenergies, indexed as energies_ks[eigenenergies]. Boolean - True if file containing exact Kohn-Sham potential is found, False if file is not found.

iDEA.RE.calculate_hartree_energy(pm, density_ks, v_h)[source]¶

Calculates the Hartree energy of the ground-state system.

\[E_{\mathrm{H}} = \frac{1}{2} \int \int U(x,x') n(x) n(x') dx dx' = \frac{1}{2} \int V_{\mathrm{H}}(x) n(x) dx\]

Parameters:	pm : object Parameters object density_ks : array_like 1D array of the ground-state Kohn-Sham electron density, indexed as density_ks[space_index] v_h : array_like 1D array of the ground-state Hartree potential, indexed as v_h[space_index] returns float The Hartree energy of the ground-state system

iDEA.RE.calculate_hartree_potential(pm, density_ks)[source]¶

Calculates the Hartree potential for a given electron density.

\[V_{\mathrm{H}}(x) = \int U(x,x') n(x') dx'\]

Parameters:	pm : object Parameters object density_ks : array_like 1D array of the Kohn-Sham electron density, either the ground-state or at a particular time step, indexed as density_ks[space_index] returns array_like 1D array of the Hartree potential for a given electron density, indexed as v_h[space_index]

iDEA.RE.calculate_time_dependence(pm, A_initial, momentum, A_ks, damping, wavefunctions_ks, density_ks, density_approx, current_density_approx)[source]¶

Calculates the exact time-dependent Kohn-Sham vector potential, at time t+dt, by solving the time-dependent Kohn-Sham equations and iteratively correcting A_ks. The exact time-dependent Kohn-Sham eigenfunctions, electron density and electron current density are then calculated.

\[A_{\mathrm{KS}}(x,t) \rightarrow A_{\mathrm{KS}}(x,t) + \nu\bigg[\frac{j_{\mathrm{KS}}(x,t)-j_{\mathrm{approx}}(x,t)} {n_{\mathrm{approx}}(x,t) + a}\bigg]\]

Parameters:

pm : object: Parameters object
A_initial : sparse_matrix: The sparse matrix A at t=0, A_initial, used when solving the equation Ax=b
momentum : array_like: 1D array of the lower band elements of the momentum matrix, index as momentum[band]
A_ks : array_like: 1D array of the time-dependent Kohn-Sham vector potential, at time t+dt, indexed as A_ks[space_index]
damping : array_like: 1D array of the damping function used to filter out noise, indexed as damping[frequency_index]
wavefunctions_ks : array_like: 2D array of the time-dependent Kohn-Sham eigenfunctions, at time t, indexed as wavefunctions_ks[space_index,eigenfunction]
density_ks : array_like: 2D array of the time-dependent Kohn-Sham electron density, at time t and t+dt, indexed as density_ks[time_index, space_index]
density_approx : array_like: 1D array of the time-dependent electron density from the approximation, at time t+dt, indexed as density_approx[space_index]
current_density_approx : array_like: 1D array of the electron current density from the approximation, at time t+dt, indexed as current_density_approx[space_index]
returns array_like and array_like and array_like and array_like and float
and float: 1D array of the time-dependent Kohn-Sham vector potential, at time t+dt, indexed as A_ks[space_index]. 1D array of the time-dependent Kohn-Sham electron density, at time t+dt, indexed as density_ks[space_index]. 1D array of the Kohn-Sham electron current density, at time t+dt, indexed as current_density_ks[space_index]. 2D array of the time-dependent Kohn-Sham eigenfunctions, at time t+dt, indexed as wavefunctions_ks[space_index,eigenfunction]. The error between the Kohn-Sham electron density and electron density from the approximation. The error between the Kohn-Sham electron current density and electron current density from the approximation.

iDEA.RE.calculate_xc_energy(pm, approx, density_ks, v_h, v_xc, energies_ks)[source]¶

Calculates the exchange-correlation energy of the ground-state system.

\[E_{\mathrm{xc}} = E_{\mathrm{total}} - \sum_{j=1}^{N}\varepsilon_{j} + \int \bigg[\frac{1}{2}V_{\mathrm{H}}(x) + V_{\mathrm{xc}}(x)\bigg]n_{\mathrm{KS}}(x)dx\]

Parameters:

pm : object: Parameters object
approx : string: The approximation used to calculate the electron density
density_ks : array_like: 1D array of the ground-state Kohn-Sham electron density, indexed as density_ks[space_index]
v_h : array_like: 1D array of the ground-state Hartree potential, indexed as v_h[space_index]
v_xc : array_like: 1D array of the ground-state exchange-correlation potential, indexed as v_xc[space_index]
energies_ks : array_like: 1D array containing the ground-state Kohn-Sham eigenenergies, indexed as energies_ks[eigenenergies]
returns float: The exchange-correlation energy of the ground-state system

iDEA.RE.construct_A(pm, A_initial, A_ks, momentum)[source]¶

Constructs the sparse matrix A at time t.

\[\begin{split}A = I + i\dfrac{\delta t}{2}H \\ H = \frac{1}{2}(p-A_{\mathrm{KS}})^{2} + V_{\mathrm{KS}}(t=0) \\ = K + V_{\mathrm{KS}}(t=0) + \frac{A_{\mathrm{KS}}^{2}}{2} - \frac{pA_{\mathrm{KS}}}{2} - \frac{A_{\mathrm{KS}}p}{2} \\ = H(t=0) + \frac{A_{\mathrm{KS}}^{2}}{2} - \frac{pA_{\mathrm{KS}}}{2} - \frac{A_{\mathrm{KS}}p}{2}\end{split}\]

Parameters:	pm : object Parameters object A_initial : sparse_matrix The sparse matrix A at t=0 A_ks : array_like 1D array of the time-dependent Kohn-Sham vector potential, at time t, indexed as A_ks[space_index] momentum : array_like 1D array of the lower band elements of the momentum matrix, index as momentum[band,space_index] returns sparse_matrix The sparse matrix A at time t

iDEA.RE.construct_A_initial(pm, K, v_ks)[source]¶

Constructs the sparse matrix A at t=0, once the external perturbation has been applied.

\[\begin{split}A_{\mathrm{initial}} = I + i\dfrac{\delta t}{2}H(t=0) \\ H(t=0) = K + V_{\mathrm{KS}}(t=0)\end{split}\]

Parameters:	pm : object Parameters object K : array_like 2D array of the kinetic energy matrix, index as K[band,space_index] v_ks : array_like 1D array of the ground-state Kohn-Sham potential + the external perturbation, indexed as v_ks[space_index] returns sparse_matrix The sparse matrix A at t=0, A_initial, used when solving the equation Ax=b

iDEA.RE.construct_K(pm)[source]¶

Stores the band elements of the kinetic energy matrix in lower form. The kinetic energy matrix is constructed using a three-point, five-point or seven-point stencil. This yields an NxN band matrix (where N is the number of grid points). For example with N=6 and a three-point stencil:

\[\begin{split}K = -\frac{1}{2} \frac{d^2}{dx^2}= -\frac{1}{2} \begin{pmatrix} -2 & 1 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 1 & -2 \end{pmatrix} \frac{1}{\delta x^2} = [\frac{1}{\delta x^2},-\frac{1}{2 \delta x^2}]\end{split}\]

Parameters:	pm : object Parameters object returns array_like 2D array containing the band elements of the kinetic energy matrix, indexed as K[band,space_index]

iDEA.RE.construct_damping(pm)[source]¶

Stores the damping function which is used to filter out noise in the Kohn-Sham vector potential.

\[f_{\mathrm{damping}}(x) = e^{-10^{-12}(\beta x)^{\sigma}}\]

Parameters:	pm : object Parameters object returns array_like 1D array of the damping function used to filter out noise, indexed as damping[frequency_index]

iDEA.RE.construct_momentum(pm)[source]¶

Stores the band elements of the momentum matrix in lower form. The momentum matrix is constructed using a five-point or seven-point stencil. This yields an NxN band matrix (where N is the number of grid points). For example with N=6 and a five-point stencil:

\[\begin{split}p = -i \frac{d}{dx}= -\frac{i}{12} \begin{pmatrix} 0 & 8 & -1 & 0 & 0 & 0 \\ -8 & 0 & 8 & -1 & 0 & 0 \\ 1 & -8 & 0 & 8 & -1 & 0 \\ 0 & 1 & -8 & 0 & 8 & -1 \\ 0 & 0 & 1 & -8 & 0 & 8 \\ 0 & 0 & 0 & 1 & -8 & 0 \end{pmatrix} \frac{1}{\delta x} = \frac{1}{\delta x} \bigg[0,\frac{2}{3},-\frac{1}{12}\bigg]\end{split}\]

Parameters:	pm : object Parameters object returns array_like 1D array of the lower band elements of the momentum matrix, index as momentum[band]

iDEA.RE.filter_noise(pm, A_ks, damping)[source]¶

Filters out noise in the Kohn-Sham vector potential by suppressing high-frequency terms in the Fourier transform.

Parameters:

pm : object: Parameters object
A_ks : array_like: 1D array of the time-dependent Kohn-Sham vector potential, at time t+dt, indexed as A_ks[space_index]
damping : array_like: 1D array of the damping function used to filter out noise, indexed as damping[frequency_index]
returns array_like: 1D array of the time-dependent Kohn-Sham vector potential, at time t+dt, indexed as A_ks[space_index]

iDEA.RE.main(parameters, approx)[source]¶

Calculates the exact Kohn-Sham potential and exchange-correlation potential for a given electron density using the reverse-engineering algorithm. This works for both a ground-state and time-dependent system.

Parameters:	parameters : object Parameters object approx : string The approximation used to calculate the electron density returns object Results object

iDEA.RE.read_current_density(pm, approx)[source]¶

Reads in the electron current density that was calculated using the selected approximation.

Parameters:	pm : object Parameters object approx : string The approximation used to calculate the electron current density returns array_like 2D array of the electron current density from the approximation, indexed as current_density_approx[time_index,space_index]

iDEA.RE.read_density(pm, approx)[source]¶

Reads in the electron density that was calculated using the selected approximation.

Parameters:	pm : objectmath Parameters object approx : string The approximation used to calculate the electron density returns array_like 2D array of the ground-state/time-dependent electron density from the approximation, indexed as density_approx[time_index,space_index]

iDEA.RE.remove_gauge(pm, A_ks, v_ks, v_ks_gs)[source]¶

Removes the gauge transformation that was applied to the Kohn-Sham potential, so that it becomes a fully scalar quantity.

\[V_{\mathrm{KS}}(x,t) \rightarrow V_{\mathrm{KS}}(x,t) + \int_{-x_{\mathrm{max}}}^{x} \frac{\partial A_{\mathrm{KS}}(x',t)}{\partial t} dx'\]

Parameters:

pm : object: Parameters object
A_ks : array_like: 2D array of the time-dependent Kohn-Sham vector potential, at time t and t+dt, indexed as A_ks[time_index, space_index]
v_ks : array_like: 1D array of the time-dependent Kohn-Sham potential, at time t+dt, indexed as v_ks[space_index]
v_ks_gs : array_like: 1D array of the ground-state Kohn-Sham potential, indexed as v_ks[space_index]
returns array_like: 1D array of the time-dependent Kohn-Sham potential, at time t +dt, indexed as v_ks[space_index]

iDEA.RE.solve_gsks_equations(pm, hamiltonian)[source]¶

Solves the ground-state Kohn-Sham equations to find the ground-state Kohn-Sham eigenfunctions, energies and electron density.

\[\begin{split}\hat{H}\phi_{j}(x) = \varepsilon_{j}\phi_{j}(x) \\ n_{\mathrm{KS}}(x) = \sum_{j=1}^{N}|\phi_{j}(x)|^{2}\end{split}\]

Parameters:

pm : object: Parameters object
hamiltonian : array_like: 2D array of the Hamiltonian matrix, index as K[band,space_index]
returns array_like and array_like and array_like: 2D array of the ground-state Kohn-Sham eigenfunctions, indexed as wavefunctions_ks[space_index,eigenfunction]. 1D array containing the ground-state Kohn-Sham eigenenergies, indexed as energies_ks[eigenenergies]. 1D array of the ground-state Kohn-Sham electron density, indexed as density_ks[space_index].

iDEA.RE.solve_tdks_equations(pm, A, wavefunctions_ks)[source]¶

Solves the time-dependent Kohn-Sham equations to find the time-dependent Kohn-Sham eigenfunctions and electron density.

\[\begin{split}\hat{H}\phi_{j}(x,t) = i\frac{\partial\phi_{j}(x,t)}{\partial t} \\ n_{\mathrm{KS}}(x,t) = \sum_{j=1}^{N}|\phi_{j}(x,t)|^{2}\end{split}\]

Parameters:

pm : object: Parameters object
A : sparse_matrix: The sparse matrix A, used when solving the equation Ax=b
wavefunctions_ks : array_like: 2D array of the time-dependent Kohn-Sham eigenfunctions, at time t+dt, indexed as wavefunctions_ks[space_index,eigenfunction]
returns array_like and array_like: 1D array of the time-dependent Kohn-Sham electron density, at time t+dt, indexed as density_ks[space_index]. 2D array of the time-dependent Kohn-Sham eigenfunctions, at time t+dt, indexed as wavefunctions_ks[space_index,eigenfunction]

iDEA.RE.xc_correction(pm, v_xc)[source]¶

Calculates an approximation to the constant that needs to be added to the exchange-correlation potential so that it asymptotically approaches zero at large $|x|$. The approximate error (standard deviation) on the constant is also calculated.

\[V_{\mathrm{xc}}(x) \rightarrow V_{\mathrm{xc}}(x) + a \ , \ \ \text{s.t.} \ \lim_{|x| \to \infty} V_{\mathrm{xc}}(x) = 0\]

Parameters:	pm : object Parameters object v_xc : array_like 1D array of the ground-state exchange-correlation potential, indexed as v_xc[space_index] returns float and float An approximation to the constant that needs to be added to the exchange-correlation potential. The approximate error (standard deviation) on the constant.

iDEA.RE.xc_fit(grid, correction)[source]¶

Applies a fit to the exchange-correlation potential over a specified range near the edge of the system’s grid to determine the correction that needs to be applied to give the correct asymptotic behaviour at large $|x|$

\[V_{\mathrm{xc}}(x) \approx \frac{1}{x} + a\]

Parameters:	grid : array_like 1D array of the spatial grid over a specified range near the edge of the system correction : float An approximation to the constant that needs to be added to the exchange-correlation potential returns array_like A fit to the exchange-correlation potential over the specified range

iDEA.RE_cython module¶

Contains the cython modules that are called within RE. Cython is used for operations that are very expensive to do in Python, and performance speeds are close to C.

iDEA.RE_cython.continuity_eqn()¶

Calculates the electron current density of the system for a particular time step by solving the continuity equation.

Parameters:	pm : object Parameters object density_new : array_like 1D array of the electron density at time t density_old : array_like 1D array of the electron density at time t-dt returns array_like 1D array of the electron current density at time t

iDEA.cli module¶

Functions for the command line interface.

iDEA.cli.examples_cli()[source]¶

Start jupyter notebook server in example directory.

If this fails, print the path to the example directory.

iDEA.cli.optimize()[source]¶

Calculate optimal alpha for hybrid Hamiltonian.

Return a tuple of three optimal values of alpha, corresponding to each condition: (LUMO-A, LUMO-HOMO, LUMO-I). If no alpha satisfies a condition, return None.

iDEA.cli.optimize_cli()[source]¶

iDEA.cli.run_cli(fname='parameters.py')[source]¶: Function to run iDEA using parameters.py file in current working directory.

iDEA.cli.video_cli()[source]¶: Produce plots, animations and data files from pickle files generated by iDEA

iDEA.info module¶

Contains information on version, authors, etc.

iDEA.info.get_sha1()[source]¶

Returns sha1 hash of last commit from git

Works only, if the code resides inside a git repository. Returns None otherwise.

iDEA.input module¶

Stores input parameters for iDEA calculations.

class iDEA.input.Input[source]¶

Bases: object

Stores variables of input parameters file

Includes automatic generation of dependent variables, checking of input parameters, printing back to file and more.

__dict__ = mappingproxy({'__module__': 'iDEA.input', '__doc__': 'Stores variables of input parameters file\n\n Includes automatic generation of dependent variables,\n checking of input parameters, printing back to file and more.\n ', 'priority_dict': {'low': 2, 'default': 1, 'high': 0}, '__init__': <function Input.__init__>, 'check': <function Input.check>, '__str__': <function Input.__str__>, 'sprint': <function Input.sprint>, 'from_python_file': <classmethod object>, 'read_from_python_file': <function Input.read_from_python_file>, 'output_dir': <property object>, 'make_dirs': <function Input.make_dirs>, 'setup_space': <function Input.setup_space>, 'execute': <function Input.execute>, '__dict__': <attribute '__dict__' of 'Input' objects>, '__weakref__': <attribute '__weakref__' of 'Input' objects>})¶

__init__(self)[source]¶: Sets default values of some properties.

__module__ = 'iDEA.input'¶

__str__(self)[source]¶: Prints different sections in input file

__weakref__¶: list of weak references to the object (if defined)

check(self)[source]¶: Checks validity of input parameters.

execute(self)[source]¶: Run this job

classmethod from_python_file(filename)[source]¶: Create Input from Python script.

make_dirs(self)[source]¶: Set up ouput directory structure

output_dir¶: Returns full path to output directory

priority_dict = {'default': 1, 'high': 0, 'low': 2}¶

read_from_python_file(self, filename)[source]¶: Update Input from Python script.

setup_space(self)[source]¶

Prepares for performing calculations

precomputes quantities on grids, etc.

sprint(self, string='', priority=1, newline=True, refresh=5e-06, savelog=True)[source]¶

Customized print function

Prints to screen and appends to log.

If newline == False, overwrites last line, but refreshes only every refresh seconds.

Parameters:	string : string string to be printed priority: int priority of message, possible values are 0: debug 1: normal 2: important newline : bool If False, overwrite the last line refresh : float If newline == False, print only every “refresh” seconds savelog : bool If True, save string to log file

class iDEA.input.InputSection[source]¶

Bases: object

Generic section of input file

__dict__ = mappingproxy({'__module__': 'iDEA.input', '__doc__': 'Generic section of input file', '__str__': <function InputSection.__str__>, '__dict__': <attribute '__dict__' of 'InputSection' objects>, '__weakref__': <attribute '__weakref__' of 'InputSection' objects>})¶

__module__ = 'iDEA.input'¶

__str__(self)[source]¶: Print variables of section and their values

__weakref__¶: list of weak references to the object (if defined)

class iDEA.input.SpaceGrid(pm)[source]¶

Bases: object

Stores basic real space arrays

These arrays should be helpful in many types of iDEA calculations. Storing them in the Input object avoids having to recompute them and reduces code duplication.

__dict__ = mappingproxy({'__module__': 'iDEA.input', '__doc__': 'Stores basic real space arrays\n\n These arrays should be helpful in many types of iDEA calculations.\n Storing them in the Input object avoids having to recompute them\n and reduces code duplication.\n ', '__init__': <function SpaceGrid.__init__>, '__str__': <function SpaceGrid.__str__>, '__dict__': <attribute '__dict__' of 'SpaceGrid' objects>, '__weakref__': <attribute '__weakref__' of 'SpaceGrid' objects>})¶

__init__(self, pm)[source]¶: Initialize self. See help(type(self)) for accurate signature.

__module__ = 'iDEA.input'¶

__str__(self)[source]¶: Print variables of section and their values

__weakref__¶: list of weak references to the object (if defined)

class iDEA.input.SystemSection[source]¶

Bases: iDEA.input.InputSection

System section of input file

Includes some derived quantities.

__module__ = 'iDEA.input'¶

deltat¶: Spacing of temporal grid

deltax¶: Spacing of real space grid

grid_points¶: Real space grid

iDEA.input.input_string(key, value)[source]¶: Prints a line of the input file

iDEA.minimize module¶

Direct minimisation of the Hamiltonian

class iDEA.minimize.CGMinimizer(pm, total_energy=None, nstates=None, cg_restart=3, line_fit='quadratic')[source]¶

Bases: object

Performs conjugate gradient minimization

Performs Pulay mixing with Kerker preconditioner, as described on pp 1071 of [Payne1992]

__dict__ = mappingproxy({'__module__': 'iDEA.minimize', '__doc__': 'Performs conjugate gradient minimization\n\n Performs Pulay mixing with Kerker preconditioner,\n as described on pp 1071 of [Payne1992]_\n ', '__init__': <function CGMinimizer.__init__>, 'exact_dirs': <function CGMinimizer.exact_dirs>, 'step': <function CGMinimizer.step>, 'line_search': <function CGMinimizer.line_search>, 'braket': <function CGMinimizer.braket>, 'steepest_dirs': <function CGMinimizer.steepest_dirs>, 'conjugate_directions': <function CGMinimizer.conjugate_directions>, 'total_energy': <function CGMinimizer.total_energy>, 'subspace_diagonalization': <function CGMinimizer.subspace_diagonalization>, '__dict__': <attribute '__dict__' of 'CGMinimizer' objects>, '__weakref__': <attribute '__weakref__' of 'CGMinimizer' objects>})¶

__init__(self, pm, total_energy=None, nstates=None, cg_restart=3, line_fit='quadratic')[source]¶

Initializes variables

Parameters:

pm: object: input parameters
total_energy: callable: function f(pm, waves) that returns total energy
nstates: int: how many states to retain. currently, this must equal the number of occupied states (for unoccupied states need to re-diagonalize)
cg_restart: int: cg history is restarted every cg_restart steps
line_fit: str: select method for line-search ‘quadratic’ (default), ‘quadratic-der’ or ‘trigonometric’

__module__ = 'iDEA.minimize'¶

__weakref__¶: list of weak references to the object (if defined)

braket(self, bra=None, O=None, ket=None)[source]¶

Compute braket with operator O

bra and ket may hold multiple vectors or may be empty. Variants:

Parameters:	bra: array_like (grid, nwf) lhs of braket O: array_like (grid, grid) operator. defaults to identity matrix ket: array_like (grid, nwf) rhs of braket

conjugate_directions(self, steepest_dirs, wfs)[source]¶

Compute conjugate gradient descent for one state

Updates internal arrays accordingly

See eqns (5.8-9) in [Payne1992]

Parameters:	steepest_dirs: array_like steepest-descent directions (grid, nwf) wfs: array_like wave functions (grid, nwf)
Returns:	cg_dirs: array_like conjugate directions (grid, nwf)

exact_dirs(self, wfs, H)[source]¶

Search direction from exact diagonalization

Just for testing purposes (you can easily end up with multiple minima along the exact search directions)

line_search(self, wfs, dirs, H, mode)[source]¶

Performs a line search along cg direction

Trying to minimize total energy

\[E(s) = \langle \psi(s) | H[\psi(s)] | \psi(s) \rangle\]

steepest_dirs(self, wfs, H)[source]¶

Compute steepest descent directions

Compute steepest descent directions and project out components pointing along other orbitals (equivalent to steepest descent with the proper Lagrange multipliers).

See eqns (5.10), (5.12) in [Payne1992]

Parameters:	H: array_like Hamiltonian matrix (grid,grid) wavefunctions: array_like wave function array (grid, nwf)
Returns:	steepest_orth: array_like steepest descent directions (grid, nwf)

step(self, wfs, H)[source]¶

Performs one cg step

After each step, the Hamiltonian should be recomputed using the updated wave functions.

Note that we currently don’t enforce the wave functions to remain eigenfunctions of the Hamiltonian. This should not matter for the total energy but means we need to performa a diagonalisation at the very end.

Parameters:	wfs: array_like (grid, nwf) input wave functions H: array_like input Hamiltonian
Returns:	wfs: array_like (grid, nwf) updated wave functions

subspace_diagonalization(self, v, H)[source]¶

Diagonalise suspace of wfs

Parameters:	v: array_like (grid, nwf) array of orthonormal vectors H: array_like (grid,grid) Hamiltonian matrix
Returns:	v_rot: array_like (grid, nwf) array of orthonormal eigenvectors of H (or at least close to eigenvectors)

total_energy(self, wfs)[source]¶

Compute total energy for given wave function

This method must be provided by the calling module and is initialized in the constructor.

class iDEA.minimize.DiagMinimizer(pm, total_energy=None)[source]¶

Bases: object

Performs minimization using exact diagonalisation

This would be too slow for ab initio codes but is something we can afford in the iDEA world.

Not yet fully implemented though (still missing analytic derivatives and a proper line search algorithm).

__dict__ = mappingproxy({'__module__': 'iDEA.minimize', '__doc__': 'Performs minimization using exact diagonalisation\n\n This would be too slow for ab initio codes but is something we can afford\n in the iDEA world.\n\n Not yet fully implemented though (still missing analytic derivatives and a\n proper line search algorithm).\n ', '__init__': <function DiagMinimizer.__init__>, 'total_energy': <function DiagMinimizer.total_energy>, 'h_step': <function DiagMinimizer.h_step>, '__dict__': <attribute '__dict__' of 'DiagMinimizer' objects>, '__weakref__': <attribute '__weakref__' of 'DiagMinimizer' objects>})¶

__init__(self, pm, total_energy=None)[source]¶

Initializes variables

Parameters:	pm: object input parameters total_energy: callable function f(pm, waves) that returns total energy nstates: int how many states to retain. currently, this must equal the number of occupied states (for unoccupied states need to re-diagonalize)

__module__ = 'iDEA.minimize'¶

__weakref__¶: list of weak references to the object (if defined)

h_step(self, H0, H1)[source]¶

Performs one minimisation step

Parameters:	H0: array_like input Hamiltonian to be mixed (banded form) H1: array_like output Hamiltonian to be mixed (banded form)
Returns:	H: array_like mixed hamiltonian (banded form)

total_energy(self, wfs)[source]¶

Compute total energy for given wave function

This method must be provided by the calling module and is initialized in the constructor.

iDEA.minimize.orthonormalize(v)[source]¶

Return orthonormalized set of vectors

Return orthonormal set of vectors that spans the same space as the input vectors.

Parameters:	v: array_like (n, m) array of m vectors in n-dimensional space

iDEA.mix module¶

Mixing schemes for self-consistent calculations

class iDEA.mix.PulayMixer(pm, order, preconditioner=None)[source]¶

Bases: object

Performs Pulay mixing

Can perform Pulay mixing with Kerker preconditioning as described on p.34 of [Kresse1996] Can also be combined with other preconditioners (see precondition.py).

__dict__ = mappingproxy({'__module__': 'iDEA.mix', '__doc__': 'Performs Pulay mixing\n\n Can perform Pulay mixing with Kerker preconditioning as described on p.34 of [Kresse1996]_\n Can also be combined with other preconditioners (see precondition.py).\n ', '__init__': <function PulayMixer.__init__>, 'update_arrays': <function PulayMixer.update_arrays>, 'compute_coefficients': <function PulayMixer.compute_coefficients>, 'mix': <function PulayMixer.mix>, 'precondition': <function PulayMixer.precondition>, '__dict__': <attribute '__dict__' of 'PulayMixer' objects>, '__weakref__': <attribute '__weakref__' of 'PulayMixer' objects>})¶

__init__(self, pm, order, preconditioner=None)[source]¶

Initializes variables

Parameters:	order: int order of Pulay mixing (how many densities to keep in memory) pm: object input parameters preconditioner: string May be None, ‘kerker’ or ‘rpa’

__module__ = 'iDEA.mix'¶

__weakref__¶: list of weak references to the object (if defined)

compute_coefficients(self, m, ncoef)[source]¶

Computes mixing coefficients

See [Kresse1996] equations (87) - (90)

\[\begin{split}A_{ij} = \langle R[\rho^j_{in}] | R[\rho_{in}^i \rangle \\ \bar{A}_{ij} = \langle \Delta R^j | \Delta R^i \rangle\end{split}\]

See [Kresse1996] equation (92)

Parameters:	m: int array index for non-delta quantities ncoef: int number of coefficients to compute

mix(self, den_in, den_out, eigv=None, eigf=None)[source]¶

Compute mix of densities

Computes new input density rho_in^{m+1}, where the index m corresponds to the index m used in [Kresse1996] on pp 33-34.

Parameters:	den_in: array_like input density den_out: array_like output density

precondition(self, f, eigv, eigf)[source]¶: Return preconditioned f

update_arrays(self, m, den_in, den_out)[source]¶

Updates densities and residuals

We need to store:

delta-quantities from i=1 up to m-1
den_in i=m-1, m
r i=m-1, m

In order to get Pulay started, we do one Kerker-only step (step 0).

Note: When self.step becomes larger than self.order, we overwrite data that is no longer needed.

Parameters:	m: int array index for non-delta quantities den_in: array_like input density den_out: array_like output density

iDEA.parameters module¶

iDEA.parameters.v_ext(x)[source]¶: Ground-state external potential

iDEA.parameters.v_pert(x)[source]¶: Perturbing potential (switched on at t=0)

iDEA.plot module¶

Plotting output quantities of iDEA

iDEA.plot.read_quantity(pm, name)[source]¶

Read a file from a pickle file in (/raw)

Parameters:	pm: parameters object parameters object name: string name of pickle file in (/raw) (e.g ‘gs_ext_den’) returns array_like data extracted from pickle file

iDEA.plot.to_anim(pm, names, data, td, dim, file_name=None, step=1)[source]¶

Outputs data to a .mp4 file in (/animations)

Parameters:

pm: parameters object: parameters object
names: list of string: names of the data to be saved (e.g ‘gs_ext_den’)
data: list of array_like: list of arrays to be plotted
td: bool: True: Time-dependent data, False: ground-state data
dim: int: number of dimentions of the data (0,1,2 or 3) (eg gs_ext_E would be 0, gs_ext_den would be 1)
file_name: string: name of output file (if None will be saved as default name e.g ‘gs_ext_den.dat’)
step: int: number of frames to skip when animating

iDEA.plot.to_data(pm, names, data, td, dim, file_name=None, timestep=0)[source]¶

Outputs data to a .dat file in (/data)

Parameters:

pm: parameters object: parameters object
names: list of string: names of the data to be saved (e.g ‘gs_ext_den’)
data: list of array_like: list of arrays to be plotted
td: bool: True: Time-dependent data, False: ground-state data
dim: int: number of dimentions of the data (0,1,2 or 3) (eg gs_ext_E would be 0, gs_ext_den would be 1)
file_name: string: name of output file (if None will be saved as default name e.g ‘gs_ext_den.dat’)
timestep: int: if td=True or data=3D specify the timestep to be saved

iDEA.plot.to_plot(pm, names, data, td, dim, file_name=None, timestep=0)[source]¶

Outputs data to a .pdf file in (/plots)

Parameters:

pm: parameters object: parameters object
names: list of string: names of the data to be saved (e.g ‘gs_ext_den’)
data: list of array_like: list of arrays to be plotted
td: bool: True: Time-dependent data, False: ground-state data
dim: int: number of dimentions of the data (0,1,2 or 3) (eg gs_ext_E would be 0, gs_ext_den would be 1)
file_name: string: name of output file (if None will be saved as default name e.g ‘gs_ext_den.pdf’)
timestep: int: if td=True or data=3D specify the timestep to be saved

iDEA.precondition module¶

Preconditioners for self-consistent calculations.

Can be used in conjunction with mixing schemes (see mix.py).

class iDEA.precondition.KerkerPreconditioner(pm)[source]¶

Bases: object

Performs Kerker preconditioning

Performs Kerker preconditioning, as described on p.34 of [Kresse1996]

__dict__ = mappingproxy({'__module__': 'iDEA.precondition', '__doc__': 'Performs Kerker preconditioning\n\n Performs Kerker preconditioning,\n as described on p.34 of [Kresse1996]_\n\n ', '__init__': <function KerkerPreconditioner.__init__>, 'precondition': <function KerkerPreconditioner.precondition>, '__dict__': <attribute '__dict__' of 'KerkerPreconditioner' objects>, '__weakref__': <attribute '__weakref__' of 'KerkerPreconditioner' objects>})¶

__init__(self, pm)[source]¶

Initializes variables

Parameters:	pm: object input parameters

__module__ = 'iDEA.precondition'¶

__weakref__¶: list of weak references to the object (if defined)

precondition(self, f, eigv, eigf)[source]¶: Return preconditioned f

class iDEA.precondition.RPAPreconditioner(pm)[source]¶

Bases: object

Performs preconditioning using RPA dielectric function

The static dielectric function as a function of x and x’ is computed in the Hartree approximation.

__dict__ = mappingproxy({'__module__': 'iDEA.precondition', '__doc__': "Performs preconditioning using RPA dielectric function\n\n The static dielectric function as a function of x and x'\n is computed in the Hartree approximation.\n\n ", '__init__': <function RPAPreconditioner.__init__>, 'precondition': <function RPAPreconditioner.precondition>, 'chi': <function RPAPreconditioner.chi>, '__dict__': <attribute '__dict__' of 'RPAPreconditioner' objects>, '__weakref__': <attribute '__weakref__' of 'RPAPreconditioner' objects>})¶

__init__(self, pm)[source]¶

Initializes variables

Parameters:	pm: object input parameters

__module__ = 'iDEA.precondition'¶

__weakref__¶: list of weak references to the object (if defined)

chi(self, eigv, eigf)[source]¶

Computes RPA polarizability

The static, non-local polarisability (aka density-potential response) in the Hartree approximation (often called RPA):

\[\chi^0(x,x') = \sum_j^{'} \sum_k^{''} \phi_j(x)\phi_k^*(x)\phi_j^*(x')\phi_k(x') \frac{2}{\varepsilon_j-\varepsilon_k}\]

where $\sum^{'}$ sums over occupied states and $\sum^{''}$ sums over empty states

Parameters:	eigv: array_like array of eigenvalues eigf: array_like array of eigenfunctions
Returns:	epsilon: array_like dielectric matrix in real space

precondition(self, r, eigv, eigf)[source]¶

Preconditioning using RPA dielectric matrix

\[\frac{\delta V(x)}{\delta \rho(x')} = DXC(x) \delta(x-x') + v(x-x')\]

Parameters:	r: array_like array of residuals to be preconditioned eigv: array_like array of eigenvalues eigf: array_like array of eigenfunctions

class iDEA.precondition.StubPreconditioner(pm)[source]¶

Bases: object

Performs no preconditioning

__dict__ = mappingproxy({'__module__': 'iDEA.precondition', '__doc__': 'Performs no preconditioning\n\n ', '__init__': <function StubPreconditioner.__init__>, 'precondition': <function StubPreconditioner.precondition>, '__dict__': <attribute '__dict__' of 'StubPreconditioner' objects>, '__weakref__': <attribute '__weakref__' of 'StubPreconditioner' objects>})¶

__init__(self, pm)[source]¶

Initializes variables

Parameters:	pm: object input parameters

__module__ = 'iDEA.precondition'¶

__weakref__¶: list of weak references to the object (if defined)

precondition(self, f, eigv, eigf)[source]¶: Return preconditioned f

iDEA.results module¶

Bundles and saves iDEA results

class iDEA.results.Results[source]¶

Bases: object

Container for results.

A convenient container for storing, reading and saving the results of a calculation.

Usage:

res = Results()
res.add(my_result, 'my_name')
res.my_name  # now contains my_result
res.save(pm)  # saves to disk + keeps track

res.add(my_result2, 'my_name2')
res.save(pm)  # saves only my_result2 to disk

__dict__ = mappingproxy({'__module__': 'iDEA.results', '__doc__': "Container for results.\n\n A convenient container for storing, reading and saving the results of a\n calculation.\n\n Usage::\n\n res = Results()\n res.add(my_result, 'my_name')\n res.my_name # now contains my_result\n res.save(pm) # saves to disk + keeps track\n\n res.add(my_result2, 'my_name2')\n res.save(pm) # saves only my_result2 to disk\n\n ", 'calc_dict': {'td': 'time-dependent', 'gs': 'ground state'}, 'method_dict': {'non': 'non-interacting', 'ext': 'exact', 'hf': 'Hartree-Fock', 'lda': 'LDA'}, 'quantity_dict': {'den': '$\\rho$', 'cur': '$j$', 'eigf': '$\\psi_j$', 'eigv': '$\\varepsilon_j$', 'vxt': '$V_{ext}$', 'vh': '$V_{H}$', 'vxc': '$V_{xc}$', 'vks': '$V_{KS}$', 'tden': '$\\rho$', 'W': '$W$', 'S': '$\\Sigma$', 'Sc': '$\\Sigma_{c}$', 'Sx': '$\\Sigma_{x}$', 'Sxc': '$\\Sigma_{xc}$'}, '__init__': <function Results.__init__>, '_not_saved': <property object>, 'label': <staticmethod object>, 'add': <function Results.add>, 'read': <staticmethod object>, 'add_pickled_data': <function Results.add_pickled_data>, 'save': <function Results.save>, 'save_hdf5': <function Results.save_hdf5>, '__dict__': <attribute '__dict__' of 'Results' objects>, '__weakref__': <attribute '__weakref__' of 'Results' objects>})¶

__init__(self)[source]¶: Initialize self. See help(type(self)) for accurate signature.

__module__ = 'iDEA.results'¶

__weakref__¶: list of weak references to the object (if defined)

_not_saved¶: Returns list of results not yet saved to disk.

add(self, results, name)[source]¶

Add results to the container.

Note: Existing results are overwritten.

add_pickled_data(self, name, pm, dir=None)[source]¶

Read results from pickle file and adds to results.

Parameters:	name : string name of results to be read (filepath = raw/name.db) pm : object iDEA.input.Input object dir : string directory where result is stored default: pm.output_dir + ‘/raw’

calc_dict = {'gs': 'ground state', 'td': 'time-dependent'}¶

static label(shortname)[source]¶

returns full label for shortname of result.

Expand shortname used for quantities saved by iDEA. E.g. ‘gs_non_den’ => ‘ground state $rho$ (non-interacting)’

method_dict = {'ext': 'exact', 'hf': 'Hartree-Fock', 'lda': 'LDA', 'non': 'non-interacting'}¶

quantity_dict = {'S': '$\\Sigma$', 'Sc': '$\\Sigma_{c}$', 'Sx': '$\\Sigma_{x}$', 'Sxc': '$\\Sigma_{xc}$', 'W': '$W$', 'cur': '$j$', 'den': '$\\rho$', 'eigf': '$\\psi_j$', 'eigv': '$\\varepsilon_j$', 'tden': '$\\rho$', 'vh': '$V_{H}$', 'vks': '$V_{KS}$', 'vxc': '$V_{xc}$', 'vxt': '$V_{ext}$'}¶

static read(name, pm, dir=None)[source]¶

Reads and returns results from pickle file

Parameters:	name : string name of result to be read (filepath = raw/name.db) pm : object iDEA.input.Input object dir : string directory where result is stored default: pm.output_dir + ‘/raw’ Returns data

save(self, pm, dir=None, list=None)[source]¶

Save results to disk.

Note: Saves only results that haven’t been saved before.

Parameters:	pm : object iDEA.input.Input object dir : string directory where to save results default: pm.output_dir + ‘/raw’ verbosity : string additional info will be printed for verbosity ‘high’ list : array_like if given, saves listed results if not set, saves results that haven’t been saved before

save_hdf5(self, pm, dir=None, list=None, f=None)[source]¶

Save results to HDF5 database.

This requires the h5py python package.

Parameters:	pm : object iDEA.input.Input object dir : string directory where to save results default: pm.output_dir + ‘/raw’ verbosity : string additional info will be printed for verbosity ‘high’ list : array_like if set, only the listed results will be saved f : HDF5 handle handle of HDF5 file (or group) to write to

iDEA.splash module¶

Prints iDEA logo as splash

iDEA.splash.draw(pm)[source]¶

iDEA.test_EXT1 module¶

Tests for 1-electron exact calculations in iDEA

class iDEA.test_EXT1.TestAtom(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests for an atomic-like potential

External potential is a softened atomic-like potential. Testing ground-state case. Testing 3-, 5- and 7-point stencil for the second-derivative.

__module__ = 'iDEA.test_EXT1'¶

setUp(self)[source]¶: Sets up atomic system

test_stencil_five(self)[source]¶: Test 5-point stencil

test_stencil_seven(self)[source]¶: Test 7-point stencil

test_stencil_three(self)[source]¶: Test 3-point stencil

class iDEA.test_EXT1.TestHarmonicOscillator(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests for the harmonic oscillator potential

External potential is the harmonic oscillator (this is the default in iDEA). Testing ground-state and time-dependence case.

__module__ = 'iDEA.test_EXT1'¶

setUp(self)[source]¶: Sets up harmonic oscillator system

test_system(self)[source]¶: Test ground-state and then real time propagation

iDEA.test_EXT2 module¶

Tests for 2-electron exact calculations in iDEA

class iDEA.test_EXT2.TestDoubleWell(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests for an asymmetric double-well potential

External potential is an asymmetric double-well potential (System 1 in Hodgson et al. 2016). Testing ground-state interacting case. Testing 3-, 5- and 7-point stencil for the second-derivative. Testing initial wavefunction by starting from the non-interacting orbitals.

__module__ = 'iDEA.test_EXT2'¶

setUp(self)[source]¶: Sets up double-well system

test_stencil_five(self)[source]¶: Test 5-point stencil

test_stencil_seven(self)[source]¶: Test 7-point stencil

test_stencil_three(self)[source]¶: Test 3-point stencil

class iDEA.test_EXT2.TestHarmonicOscillator(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests for the harmonic oscillator potential

External potential is the harmonic oscillator (this is the default in iDEA). Testing both ground-state non-interacting and ground-state interacting case. Testing time-dependent interacting case.

__module__ = 'iDEA.test_EXT2'¶

setUp(self)[source]¶: Sets up harmonic oscillator system

test_interacting_system_1(self)[source]¶: Test interacting system

test_non_interacting_system_1(self)[source]¶: Test non-interacting system

test_time_dependence(self)[source]¶: Test real time propagation

iDEA.test_EXT3 module¶

Tests for 3-electron exact calculations in iDEA

class iDEA.test_EXT3.TestDoubleWell(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests for an asymmetric double-well potential

External potential is an asymmetric double-well potential (System 1 in Hodgson et al. 2016). Testing ground-state interacting case. Testing 3-, 5- and 7-point stencil for the second-derivative. Testing initial wavefunction by starting from the non-interacting orbitals.

__module__ = 'iDEA.test_EXT3'¶

setUp(self)[source]¶: Sets up double-well system

test_stencil_five(self)[source]¶: Test 5-point stencil

test_stencil_seven(self)[source]¶: Test 7-point stencil

test_stencil_three(self)[source]¶: Test 3-point stencil

class iDEA.test_EXT3.TestHarmonicOscillator(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests for the harmonic oscillator potential

External potential is the harmonic oscillator (this is the default in iDEA). Testing both ground-state non-interacting and ground-state interacting case. Testing time-dependent interacting case.

__module__ = 'iDEA.test_EXT3'¶

setUp(self)[source]¶: Sets up harmonic oscillator system

test_interacting_system_1(self)[source]¶: Test interacting system

test_non_interacting_system_1(self)[source]¶: Test non-interacting system

test_time_dependence(self)[source]¶: Test real time propagation

iDEA.test_HF module¶

Tests for the local density approximation

class iDEA.test_HF.HFTestHarmonic(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests on the harmonic oscillator potential

__module__ = 'iDEA.test_HF'¶

setUp(self)[source]¶: Sets up harmonic oscillator system

test_1_electron(self)[source]¶: Ensures HF is exact for one electron

iDEA.test_HYB module¶

Tests for the local density approximation

class iDEA.test_HYB.HYBTestHarmonic(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests on the harmonic oscillator potential

__module__ = 'iDEA.test_HYB'¶

setUp(self)[source]¶: Sets up harmonic oscillator system

test_alpha(self)[source]¶: Ensures HYB is idential to HF if a=1.0

iDEA.test_LDA module¶

Tests for the local density approximation

class iDEA.test_LDA.LDATestHarmonic(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests on the harmonic oscillator potential

__module__ = 'iDEA.test_LDA'¶

setUp(self)[source]¶: Sets up harmonic oscillator system

test_banded_hamiltonian_1(self)[source]¶

Test construction of Hamiltonian

Hamiltonian is constructed in banded form for speed. This checks that the construction actually works.

test_kinetic_energy_1(self)[source]¶

Checks kinetic energy

Constructs Hamiltonian with KS-potential set to zero and computes expectation values.

test_total_energy_1(self)[source]¶

Compares total energy computed via two methods

One method uses the energy eigenvalues + correction terms. The other uses it only for the kinetic part, while the rest is computed usin energy functinals.

iDEA.test_NON module¶

Tests for non-interacting calculations in iDEA

class iDEA.test_NON.TestAtom(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests for an atomic-like potential

External potential is a softened atomic-like potential. Testing ground-state case. Testing 3-, 5- and 7-point stencil for the second-derivative.

__module__ = 'iDEA.test_NON'¶

setUp(self)[source]¶: Sets up atomic system

test_stencil_five(self)[source]¶: Test 5-point stencil

test_stencil_seven(self)[source]¶: Test 7-point stencil

test_stencil_three(self)[source]¶: Test 3-point stencil

class iDEA.test_NON.TestHarmonicOscillator(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests for the harmonic oscillator potential

External potential is the harmonic oscillator (this is the default in iDEA). Testing ground-state and time-dependence case.

__module__ = 'iDEA.test_NON'¶

setUp(self)[source]¶: Sets up harmonic oscillator system

test_system(self)[source]¶: Test ground-state and then real time propagation

iDEA.test_minimize module¶

Tests for direct minimizers

class iDEA.test_minimize.TestCG(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests for the conjugate gradient minimizer

__module__ = 'iDEA.test_minimize'¶

setUp(self)[source]¶: Sets up harmonic oscillator system

test_conjugate(self)[source]¶

Check that gradient is computed correctly

Use conjugate-gradient method on a quadratic function in n-dimensional space, where it is guaranteed to converge in n steps.

Note: This exact condition actually doesn’t apply here because we are orthonormalising the vectors (i.e. rotating) after each step. This renders this test a bit pointless… I am currently taking 2*n steps and am still far from machine precision.

Task: minimize the function

\[E(v) = \sum_{i=1}^2 v_i^* H v_i \triangle_{v_i} E = H v_i\]

for a fixed matrix H and a set of orthonormal vectors v_i.

Note: The dimension n of the vector space is the grid spacing times number of electrons.

test_orthonormalisation(self)[source]¶

Testing orthonormalisation of a set of vectors

minimize. This should do the same as Gram-Schmidt

test_steepest_dirs(self)[source]¶

Testing orthogonalisation in steepest descent

Just checking that the efficient numpy routines do the same as more straightforward loop-based techniques

class iDEA.test_minimize.TestCGLDA(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests for the conjugate gradient minimizer

On an actual LDA system

__module__ = 'iDEA.test_minimize'¶

setUp(self)[source]¶: Sets up harmonic oscillator system

test_energy_derivative(self)[source]¶: Compare analytical total energy derivative vs finite differences

\[\frac{dE}{d\lambda}(\lambda) \approx \frac{E(\lambda)+E(\lambda+\delta)}{\delta} \]

iDEA.test_mix module¶

Tests for mixing schemes

class iDEA.test_mix.TestKerker(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests for the Kerker preconditioner

__module__ = 'iDEA.test_mix'¶

setUp(self)[source]¶: Sets up harmonic oscillator system

test_screening_length_1(self)[source]¶

Testing screening length in Kerker

Check that for infinite screening length, simple mixing is recovered. [Kresse1996] p.34 …

class iDEA.test_mix.TestPulay(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests for the Pulay mixer

__module__ = 'iDEA.test_mix'¶

setUp(self)[source]¶: Sets up harmonic oscillator system

test_array_update_1(self)[source]¶

Testing internal variables of Pulay mixer

Just checking that the maths works as expected from [Kresse1996] p.34 …

class iDEA.test_mix.TestRPA(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests for the RPA preconditioner

__module__ = 'iDEA.test_mix'¶

setUp(self)[source]¶: Sets up harmonic oscillator system

test_chi_1(self)[source]¶

Testing potential-density response

Testing some basic symmetry properties of the potential-density response and the preconditioning matrices required for density/potential mixing.

iDEA.test_result module¶

Tests for the result class

class iDEA.test_result.resultsTest(methodName='runTest')[source]¶

Bases: unittest.case.TestCase

Tests results object

__module__ = 'iDEA.test_result'¶

setUp(self)[source]¶: Sets up harmonic oscillator system

test_save_1(self)[source]¶: Checks that saving works as expected

Module contents¶

interacting Dynamic Electrons Approach (iDEA)

The iDEA code allows to propagate the time-dependent Schroedinger equation for 2-3 electrons in one-dimensional real space. Compared to other models, such as the Anderson impurity model, this allows us to treat exchange and correlation throughout the system and provides additional flexibility in bridging the gap between model systems and ab initio descriptions.