step3¶

libra_py.workflows.nbra.step3.apply_normalization(S, St)[source]¶

Transforms the input transition density matrix computed with potentially non-orthogonalized orbitals such that it would correspond to the properly orthonormalized ones

Parameters

S (CMATRIX(2N, 2N)) –
is a matrix of MO overlaps S_ij = <i|j>. It has a block structure as:

\[S = \begin{vmatrix} S_{aa} & S_{ab} \ S_{ba} & S_{bb} \end{vmatrix}\]

Here, S_xy are the overlaps of the MOs for spin channels x and y (alpha, beta) - only spatial components of the orbitals are taken into account here. Here, N - is the total number of orbitals (double occupancies)
St (CMATRIX(2N, 2N)) –
the transition density matrix St_ij = <i|d/dt|j>. It has a block structure as:

\[St = \begin{vmatrix} St_{aa} & St_{ab} \ St_{ba} & St_{bb} \end{vmatrix}\]

Here, St_xy are the transition density matrix for spin channels x and y (alpha, beta) - only spatial components of the orbitals are taken into account here. Here, N - is the total number of orbitals (double occupancies)

Returns

but the input matrix `St` is changed

Return type

None

libra_py.workflows.nbra.step3.apply_orthonormalization_general(S, St)[source]¶

Transforms the input transition density matrix computed with potentially non-orthogonalized orbitals such that it would correspond to the properly orthonormalized ones

Parameters

S (CMATRIX(N, N)) – is a matrix of MO overlaps S_ij = <i|j>
St (CMATRIX(N, N)) – the transition density matrix St_ij = <i|d/dt|j>

Returns

but the input matricies `S` and `St` are changed

Return type

None

libra_py.workflows.nbra.step3.apply_phase_correction(St)[source]¶

Performs the phase correction according to: Akimov, A. V. J. Phys. Chem. Lett, 2018, 9, 6096

Parameters

St (list of CMATRIX(N,N)) –

St_ij[n] = <i(n)|j(n+1)> transition density matrix for the timestep n, where N is the number of spin-orbitals in the active space. Spin-orbitals, not just orbitals! So it is composed as:

\[St = \begin{vmatrix} St_{aa} & St_{ab} \ St_{ba} & St_{bb} \end{vmatrix}\]

Returns

but changes the input St matrices

Return type

None

libra_py.workflows.nbra.step3.apply_phase_correction_general(St)[source]¶

Performs the phase correction according to: Akimov, A. V. J. Phys. Chem. Lett, 2018, 9, 6096

This function is for dat NOT in spin-block format

Parameters: St (list of CMATRIX(N,N)) – St_ij[n] = <i(n)|j(n+1)> transition density matrix for the timestep n, where N is the number of states in the active space. Spin-orbitals, not just orbitals! So it is composed as:
Returns: but changes the input St matrices
Return type: None

libra_py.workflows.nbra.step3.apply_state_reordering(St, E, params)[source]¶

Performs the state’s identity reordering in a given basis for all time steps. This is reflects in the corresponding changess of the TDM.

Parameters

St (list of CMATRIX(nstates, nstates)) – TDM for each timestep
E (list of CMATRIX(nstates, nstates)) – energies of all states at every step
params (dictionary) –
parameters controlling the reordering * params[“do_state_reordering”] ( int ): option to select the state reordering algorithm
Available options:
1: older version developed by Kosuke Sato, may not the working all the times

2: Munkres-Kuhn (Hungarian) method [default]
- params[“state_reordering_alpha”] ( double ): a parameter that controls how
  many states will be included in the reordering

Returns

but changes the input St object

Return type

None

libra_py.workflows.nbra.step3.apply_state_reordering_general(St, E, params)[source]¶

Performs the state’s identity reordering in a given basis for all time steps. This is reflects in the corresponding changess of the TDM.

This function is for dat NOT in spin-block format

Parameters

St (list of CMATRIX(nstates, nstates)) – TDM for each timestep
E (list of CMATRIX(nstates, nstates)) – energies of all states at every step
params (dictionary) –
parameters controlling the reordering * params[“do_state_reordering”] ( int ): option to select the state reordering algorithm
Available options:
1: older version developed by Kosuke Sato, may not the working all the times

2: Munkres-Kuhn (Hungarian) method [default]
- params[“state_reordering_alpha”] ( double ): a parameter that controls how
  many states will be included in the reordering

Returns

but changes the input St object

Return type

None

libra_py.workflows.nbra.step3.build_SD_basis(data_dim, cbm_alpha_index, alpha_include, cbm_beta_index, beta_include, excitation_type)[source]¶

Builds a Slater Determinant basis based on the indexing notation scheme used in Libra

Parameters

data_dim (int) – how many rows or columns in the vibronic Hamiltonian matrix. This will be an even number becuase the number of alpha orbtials should equal the number of beta orbitals
cbm (alpha/beta)_index (int) – index of VBM (or HOMO) in the matrix of the vibronic Hamiltonian (row or column index). Note, this index is from 1
(alpha/beta)_include (int) – how many orbitals to include from the cbm_(alpha/beta)_index
excitation_type (int) – 0: Make SDs with beta electrons excited 1: Make SDs with alpha electrons excited 2: Make two sets of SDs, one for beta and one for alpha electrons excited

libra_py.workflows.nbra.step3.compute_Hvib(basis, St_ks, E_ks, dE, dt)[source]¶

Compute the vibronic Hamiltonian matrix

Parameters

basis (list of lists of integers) –
defines the basis of Slater Determinants, such that: basis[iSD][iks] is the indicator of the spin-orbital occupied by the electron iks in the Slater Determinant iSD

Example

The following example defines a ground state SD (the lowest KS of the active space) and two single excitations, which are different from each other by two spin flips of the electrons The convention is to start indexing from 1 (corresponds to index 0 in the KS matrices) Positive - for alpha electrons, negative - for beta electrons Need to be consistent: [ -1, 2 ] and [ 2, -1 ] are treated differently, this is needed for spin-adaptation

>> basis = [ [ 1,-1 ], [ 1,-2 ], [ 2,-1 ] ]

The next example is for a system of 4 electrons and hole excitations >> basis = [ [ 1,-1, 2, -2 ], [ 3, -1, 2, -2 ], [ 1, -3, 2, -2 ] ]
St_ks (CMATRIX(2*norbs, 2*norbs)) – transition density matrix in the KS spin-orbitals basis, where norb - the number of double-occupied orbitals.
E_ks (CMATRIX(2*norbs, 2*norbs)) – the orbital energies in the KS spin-orbitals basis, where norb - the number of double-occupied orbitals.
dE (list of doubles) –
define corrections of the SD state energies in comparison to the energy give by the sum energies of the occupied spin-orbitals. The convention is: dE[iSD] is the correction to energy of the SD with index iSD. This is a constant correction - same for all energies in the set [units: Ha]

Example

For instance, for the SD examples above, the corrections could be: >> dE = [0.0, 0.01, 0.05]
dt (double) – the timestep for MD integrations [units: a.u.]

Returns

CMATRIX(nstates, nstates) The Vibronic Hamiltonian

libra_py.workflows.nbra.step3.compute_ci_energies_midpoint(ci_energies, num_excited_states, istep, fstep)[source]¶

This function compute the excitation energies energies at the midpoint from a list of excitation energies at each step. At each step, there are many electronic states. This function takes a list as an input, and is meant to be used in the NBRA workflow calculatiosn where lists may be more convenient than matricies.

This funciton is made to be used within the NBRA Libra workflow, where things such as ci_energies have been extracted from TD-DFT calculations. As of 11/30/2020, compatable ES programs include CP2K, DFTB+ and Gaussian.

Energies are assumed to be energies from TDDFT calculatons. This function gives zero as the ground state total energy

Parameters

ci_energies (list of lists) – energies of the MB states
num_excited_states (int) – number of excited states
istep (int) – step at which to start counting
fstep (int) – stap at which to stop counting

Returns

energies in Ha. Ground state energy is set to zero

Return type

ci_midpoint_energies (list of CMATRIX)

libra_py.workflows.nbra.step3.do_phase_corr(cum_phase1, St, cum_phase2, phase_i)[source]¶

This function changes the St matrix according to the previous cumulative phases and the current phase correction as:

St = <bra|ket>

St -> St = F_n * St * (F_{n+1})^+ = F_n * St * (F_{n})^+ * (f_{n+1})^+

Parameters

cum_phase1 (CMATRIX(nstates, 1)) – cumulative phase corrections up to step n (F_n) for bra-vectors
St (CMATRIX(nstates, nstates)) – input/output TDM to be processed: could be alpha-alpha, beta-beta, alpha-beta, or beta-alpha sub-blocks
cum_phase2 (CMATRIX(nstates, 1)) – cumulative phase corrections up to step n (F_n) for ket-vectors
phase_i (CMATRIX(nstates, 1)) – the current step phase corrections (f_{n+1}) for a given pair of vectors

Returns

but changes the input matrix St

Return type

None

libra_py.workflows.nbra.step3.get_Lowdin(S)[source]¶

Find the S_i_half for the S matrix - alpha and beta components

Parameters

S (CMATRIX(2N, 2N)) –

is a matrix of MO overlaps. It has a block structure as:

\[S = \begin{vmatrix} S_{aa} & S_{ab} \ S_{ba} & S_{bb} \end{vmatrix}\]

Here, S_xy are the overlaps of the MOs for spin channels x and y (alpha, beta) - only spatial components of the orbitals are taken into account here. Here, N - is the total number of orbitals (double occupancies)

Returns

(S_aa_i_half, S_bb_i_half), where:

S_aa_i_half ( CMATRIX(N,N) ): S_aa^{-1/2} - inverse square root matrix for the alpha-alpha block

S_bb_i_half ( CMATRIX(N,N) ): S_bb^{-1/2} - inverse square root matrix for the beta-beta block

Return type

tuple

libra_py.workflows.nbra.step3.get_Lowdin_general(S)[source]¶

Find the S_i_half for the S matrix :param S: is a matrix of MO overlaps. :type S: CMATRIX(N, N)

Returns

S_i_half, where:

S_i_half ( CMATRIX(N,N) ): S^{-1/2} - inverse square root matrix

Return type

tuple

libra_py.workflows.nbra.step3.get_step2_data(_params)[source]¶

A light function to obtain the step2 data

Parameters

params (dictionary) –

Control paramerter of this type of simulation. Can include the follwing keys:

params[“basis”] ( string ): describes if one is using either the spin-diabatic (non spin-orbit coupling)
or is using the spin-adiabatic (spin orbit-coupling)

libra_py.workflows.nbra.step3.make_T_matricies(ci_coefficients, ci_basis_states, spin_components, sd_states_unique_sorted, nstates, istep, fstep, outdir, verbose=1)[source]¶

This function makes the “T”ransformation matricies that convert between the SD basis to the CI-like (or many-body (MB)) basis.

This funciton is made to be used within the NBRA Libra workflow, where things such as ci_coefficients, ci_basis_states, spin_components, and sd_states_unique_sorted have been extracted from TD-DFT calculations. As of 11/30/2020, compatable ES programs include CP2K, DFTB+ and Gaussian.

Parameters

ci_coefficients (list of lists of lists) – coefficients for the many-body states for each step
ci_basis_states (list of lists) – All SD basis states that comprise the many-body excitations for each step
spin_components (list of lists) – the spin components of the excitation (alpha or beta excitaiton?) for all states and all steps
sd_basis_states_unique (list) – 1 of each of the SP transitions (and its spin) that made up the considered CI states
nstates (int) – number of excited MB states
istep (int) – step at which to start counting
fstep (int) – stap at which to stop counting
outdir (string) – output directory for the T matricies
verbose (int) – want to see some messages?

Returns

CMATRIX at each timestep where the rows are SDs and the cols are MB states. The columns contain the coefficients of the MB expansion for each MB state

Return type

SD2CI (list of CMATRIX)

libra_py.workflows.nbra.step3.make_cost_mat(orb_mat_inp, en_mat_inp, alpha)[source]¶

Makes the cost matrix from a given TDM and information on states’ energies

Parameters

orb_mat_inp (CMATRIX(nstates,nstates) or MATRIX(nstates,nstate)) – the transition density matrix TDM in a given basis. Here, `nstates` - the number of states (e.g. the number of doubly-occupied orbitals )
en_mat_inp (MATRIX(nstates, nstates)) – Matrix of energies in a given basis [units: a.u.]
alpha (float) – Parameter controlling the range of the orbitals that can participate in the reordering. Setting is to 0 makes all orbitals be considered for reordering Setting it to a large number makes the effective number of orbitals participating in the reordering smaller - this can be used to turn off the reordering. [units: a.u.^-1]

Returns

the matrix of the cost values for different pairs of states

Return type

MATRIX(nstates, nstates)

libra_py.workflows.nbra.step3.output_sorted_Hvibs(Hvib, orbital_index_energy_pairs)[source]¶

This function outputs the vibronic Hamiltonians in the SD basis according to their sorted order

Parameters

Hvib (list of lists of CMATRIX objects) – vibronic Hamiltonian in the Slater determinant basis
orbital_index_energy_pairs (list of lists of lists of lists) – orbital index and energy pair lists for each step and nuclear trajectory Ex. orbital_index_energy_pairs[i][j][k][0] = kth slater determinant index at the jth step on the ith nuclear trajectory Ex. orbital_index_energy_pairs[i][j][k][1] = kth slater determinant energy at the jth step on the ith nuclear trajectory

libra_py.workflows.nbra.step3.print_SD_basis(SD_basis)[source]¶

Just a light function to print the SD basis

Parameters

SD_basis (list of lists of ints) –

Slater determinant basis in terms of Kohn-Sham orbital indicies

Possible ground state configurations

Ex. 1. SD_basis[0] = [ 5, -15 ] Ex. 2. SD_basis[0] = [ 2, 3, 4, 5, -12, -13, -14, -15 ]

Possible ground state configurations

Ex. 1. SD_basis[N > 0] = [ 6, -15 ] Ex. 2. SD_basis[N > 0] = [ 2, 3, 4, 6, -12, -13, -14, -15 ]

Where we have 10 alpha and beta spin orbitals, and the cbm index for the alpha spin-channel is 5 and the cbm index for the beta spin-channel is 15.

libra_py.workflows.nbra.step3.pyxaid2libra(Hvib_pyxaid, params)[source]¶

libra_py.workflows.nbra.step3.run(S_dia_ks, St_dia_ks, E_dia_ks, params)[source]¶

The procedure to converts the results of QE calculations (KS orbital energies and time-overlaps = transition density matrices in the KS basis) to the generic Hvib matrices, which (optionally) account for:

enforces orthogonalization of the input KS states
state reordering
phase corrections
multi-electron wavefunction (Slater determinants) and spin-adaptation

Parameters

S_dia_ks (list of lists of CMATRIX objects) – overlaps of the KS orbitals along trajectories for each data set. Such that S_dia_ks[idata][istep].get(i,j) is <i(istep)|j(istep)> for the trajectory (=dataset) `idata`.
St_dia_ks (list of lists of CMATRIX objects) – time-overlaps (=transition density matrices) in the basis of the KS orbitals along trajectories for each data set. Such that St_dia_ks[idata][istep].get(i,j) is <i(istep)|j(istep+1)> for the trajectory (=dataset) `idata`.
E_dia_ks (list of lists of CMATRIX objects) – energies the KS orbitals at the mid-points along trajectories for each data set. Such that E_dia_ks[idata][istep].get(i,i) is 0.5*(E_i(istep) + E_i(istep+1)) for the trajectory (=dataset) `idata`
params (dictionary) –
Control paramerter of this type of simulation. Can include the follwing keys:
- params[“SD_basis”] ( list of lists of ints ): define the Slater Determinants basis
  The convention is: params[“SD_basis”][iSD][iks] is the indicator of the spin-orbital occupied by the electron iks in the Slater Determinant iSD [required!]
  
  Example:
  
  The following example defines a ground state SD (the lowest KS of the active space) and two single excitations, which are different from each other by two spin flips of the electrons The convention is to start indexing from 1 (corresponds to index 0 in the KS matrices) Positive - for alpha electrons, negative - for beta electrons Need to be consistent: [ -1, 2 ] and [ 2, -1 ] are treated differently, this is needed for spin-adaptation
  
  >> params[“SD_basis”] = [ [ 1,-1 ], [ 1,-2 ], [ 2,-1 ] ]
  
  The next example is for a system of 4 electrons and hole excitations >> params[“SD_basis”] = [ [ 1,-1, 2, -2 ], [ 3, -1, 2, -2 ], [ 1, -3, 2, -2 ] ]
- params[“SD_energy_corr”] ( list of doubles ): define corrections of the SD state energies in comparison to
  the energy give by the sum energies of the occupied spin-orbitals. The convention is: params[“SD_energy_corr”][iSD] is the correction to energy of the SD with index iSD. This is a constant correction - same for all energies in the set [units: Ha] [required!]
  
  Example:
  For instance, for the SD examples above, the corrections could be: >> params[“SD”] = [0.0, 0.01, 0.05]
- params[“CI_basis”] ( list of lists of complex number ): configuration interaction coefficients
  that define a superpositions to SDs that are considered the states of interest, e.g. spin-adapted configurations The convention is: params[“CI_basis”][iCI][iSD] is a coefficient of `iSD`-th SD in the expansion of the CI with index `iCI`. These coefficients don’t have to account for the overal CI’s normalization - the normalization will be done on the go. [required!]
  
  Example:
  
  For the SD example above we can construct the following combinations: >> params[“CI_basis”] = [ [1.0, 0.0, 0.0 ],
  
  [0.0, 1.0,-1.0 ], [0.0, 1.0, 1.0 ]
  
  ]
- params[“output_set_paths”] ( list of strings ): the directory pathes where the resulting files
  are to be written (if so!). If you don’t plan on writing the files, just provide a list of empty strings or whatever else - they will not be used in that case. The number of the strings should be equal to the number of the input data sets, e.g. to len(St_dia_ks) [required!]
- params[“dt”] ( double ): nuclear dynamics integration time step [units: a.u. of time, default: 41.0]
- params[“do_orthogonalization”] ( int ): the option to do Lowdin orthogonalization of the orbitals - using
  the “raw” overlaps (at the same time). This option is needed because the wavefunction output by QE are not exactly orthonormal (because of the use of pseudopotentials). So before we use them (implicitly) in the rest of the calculations here, we may need to account for this non-ideality effect. Options:
  - 0: don’t do the orthogonalization - this is the same as in Pyxaid [default]
  - 1: do the orthogonalization
- params[“do_state_reordering”] ( int ): the option to control the state reordering
  Options:
  - 0: no state reordering - same as in Pyxaid
  - 1: older method (is not robust, may or may not work)
  - 2: Hungarian algorithm [default]
- params[“state_reordering_alpha”] ( double ): the parameter that controls the width of
  the energy interval within wich the state reordering is in effect. Zero value means all available orbitals, larger positive value decreases the width of the window. This parameter is not in effect unless the Hungarian algorithm is selected [default: 0.0]
- params[“do_phase_correction”] ( int ): option to do the phase correction
  - 0 - don’t do
  - 1 - do it [default]
- params[“do_output”] ( int ): whether to print out the Hvib matrices ( = 1) to the files or not ( = 0).
- params[“Hvib_re_prefix”] ( string ): common prefix of the output files with real part of the vibronic
  Hamiltonian at all times [default: “Hvib_”]
- params[“Hvib_re_suffix”] ( string ): common suffix of the output files with real part of the vibronic
  Hamiltonian at all times [default: “_re”]
- params[“Hvib_im_prefix”] ( string ): common prefix of the output files with imaginary part of the vibronic
  Hamiltonian at all times [default: “Hvib_”]
- params[“Hvib_im_suffix”] ( string ): common suffix of the output files with imaginary part of the vibronic
  Hamiltonian at all times [default: “_im”]

Returns

Hvib, such that:: Hvib[idata][istep] is a CMATRIX(N,N) containing the vibronic Hamiltonian for the trajectory (dataset) `idata` and for the timestep `istep`. Here, N is the number of states included in the active space.

Return type

list of lists of CMATRIX(N,N)

libra_py.workflows.nbra.step3.sac_matrices(coeff, basis, S_ks)[source]¶

This function makes the Phi-to-Chi (P2C) transformation matrix. Normalization factros for the Chi states are computed based on the overlaps of Phi states. < Chi_i | Chi_j > = 1 = N_i * N_j * < Phi_i - Phi_i’ | Phi_j - Phi_j’ > = 1

coeff [List of lists] - P2C as initialized by the user basis [Phi basis] - as initialized by the user S_ks [CMATRIX] - Time overlap matrix of elementary KS orbtials, from step2

libra_py.workflows.nbra.step3.scale_H_vib(hvib, en_gap, dNAC, sc_nac_method)[source]¶

This function scales the energies and NACs in the vibrionic Hamiltonian in the Chi basis.

Parameters

hvib (list of CMATRIX objects) – CMATRIXlist of vibronic hamiltonians in the Chi basis
en_gap (float) – The desired energy gap (E_1 - E_0), for the Chi basis
dNAC (list of lists of (list, float)) –
The scaling terms by which specific nacs will be scaled datatype = list of lists of (list, float)

[ [ [i,j], val ], … ]

n and n+1 are the col (and thereby row) indicies of the nacs to be scaled by the value val
sc_nac_method (int) –
The method used to scale NACs in the Chi basis, chosen by the user. If sc_nac_method = 1, then the NACs are scaled by the ivnerse of the magnitude of the change in energy, according to Lin et al.

Reference: Lin, Y. & Akimov, A. V. J. Phys. Chem. A (2016)

libra_py.workflows.nbra.step3.sort_SD_energies(Hvib)[source]¶

This function goes into the Hvib (SD basis) files and sorts the energies For each Hvib file, we are going to obtain a list of lists of orbital index and energy pair These orbital index and energy pairs for each step will be sorted based on energies

Parameters: Hvib (list of lists of CMATRIX objects) – vibronic Hamiltonian in the Slater determinant basis

libra_py.workflows.nbra.step3.sort_unique_SD_basis(E_ks, sd_states_unique, sd_states_reindexed, istep, fstep, sorting_type)[source]¶

This function computes the energies of the SP transitions (according to the sum of 1 electron terms) - no J or K It then may sort the order of the sd_states either based on their energy at each timestep

Parameters

E_ks (list of CMATRIX) –
KS orbital energies at each timestep. Spin block style Ex) [ alp*alp alp*bet ]

[ bet*alp bet*bet ]
sd_states_unique (list of lists) – all SP transitions and which spin it was Ex) [ [ [‘28 29’], [‘alp’] ]. [ [‘28 30’], [‘alp’] ] ]
sd_states_reindexed (list of lists) – sd_states_unique but in internal Libra notation Ex) [ [1,-1,3,-2], [3,-1,2,-2] ]
sorting_type ((string)) – “energy” - sort by energy “identity” - sort by identity
istep (int) – step from which to start counting
fstep (int) – step at which to stop counting

Returns

SD energies at each timestep sd_states_unique_sorted (list of lists): All SP transitions and which spin it is, but now sorted either by identity (no sorting) or energy sd_states_reindexed_sorted (list of lists): The sd_states_unique_sorted, but in Libra’s notation reindex_nsteps (list of lists): The energy ordering of the SD for each step in terms of the index of the SD from the initial step

Return type

E_sd (list of CMATRIX)