pytdscf package

Subpackages

Submodules

pytdscf.ase_handler module

Multiprocessing mesh electronic structure calculation caller using ASE library

class pytdscf.ase_handler.DVR_Mesh(dvr_prims: list[DVRPrimitivesMixin], atoms: Atoms, disp_vec: ndarray, unit: str = 'angstrom')[source]

Bases: object

DVR grid coordinate

Parameters:

  • dvr_prims (List[DVRPrimitivesMixin]) – DVR primitive list

  • atoms (List[List[str, Tuple[float,float,float]]] or ase.atoms.Atoms) – reference coordinate. Format is the same as PySCF ones

  • disp_vec (np.ndarray) – displacement vectors (row vectors) in angstrom

  • unit (Optional[bool]) – Input reference coordinate unit. Defaults to ‘angstrom’

displace: dict
done_jobs: deque
execute_multiproc(calculator: Calculator, max_workers: int | None = None, timeout: float = 60.0, jobname: str | None = None, reset_calc: bool = False, judge_func: Callable | None = None)[source]

Execute electronic structure calculation by multiprocessing

Parameters:

  • calculator (Calculator) – calculator for each geomtry

  • max_workers (Optional[int]) – maximum workers in multi-processing. Defaults to None. If None, use cpu_count - 1.

  • timeout (float) – Timeout calculation in second. Defaults to 60.0

  • jobname (Optional[str]) – jobname

  • reset_calc (Optional[bool]) – set new calculator in any case. Defaults to False.

  • judge_func (Optional[Callable[[Any],bool]]) – judge function whether re-calculation is needed. Defaults to None.

geometry: dict
grid_id: dict
jobname: str
remain_jobs: deque
reset_calc: bool
save_geoms(jobname: str, nMR: int | None = None, overwrite: bool | None = None) → dict[str, dict[str, int]][source]

Generate cartesian coordinate geometry for each grid mesh.

Parameters:

  • nMR (Optional[int]) – Tne number of mode representation. limits n dimensional mesh. Defaults to None, thus, ngrid**ndof coords will be generated.

  • overwrite (Optional[bool]) – overwrite detabase

Returns:

DVR Mesh coordinates. E.g. [(0,1)][(2,3)] gives 2nd, 3rd grids of 0-mode, 1-mode coordinate.

Return type:

Dict[str, Dict[str, int]]

thrown_jobs: deque

pytdscf.dvr_operator_cls module

Discrete Variable Representation Operator Class

class pytdscf.dvr_operator_cls.PotentialFunction(DOFs: list[int], polynomial_coef: dict[tuple[int, ...], float], cut_off: float = -1.0, div_factorial: bool = True, dipole: bool = False, efield: tuple[float, float, float] = (0.01, 0.01, 0.01))[source]

Bases: object

Analytical Model Potential from polynomial-based PES

Parameters:

  • DOFs (List[int]) – Degree of freedoms used for potential function. 1-index.

  • polynomial_coef (Dict[Tuple[int], float]) – Polynomial PES coefficient.

  • cut_off (float) – potential coefficient cut-off. defaults to no-cutoff.

  • div_factorial (bool) – Whether or not divided by factorial term. Defaults to True.

  • dipole (Optional[bool]) – Use dipole moment surface (3D vector).

  • efiled (Optional[tuple[float, float, float]]) – Electronic field. (only dipole)

class pytdscf.dvr_operator_cls.TensorOperator(*, shape: tuple[int, ...] | None = None, tensor: ndarray | None = None, only_diag: bool = False, legs: tuple[int, ...] | None = None, name: str | None = None, mpo: list[ndarray] | None = None)[source]

Bases: object

Tensor Operator class

shape

Tensor shape

Type:

Tuple[int]

tensor_orig

Original tensor

Type:

np.ndarray

tensor_decomposed

Decomposed tensor

Type:

np.ndarray

tensor_full

if only_diag==False, tensor_orig is zero-filled.

Type:

np.ndarray

only_diag

Only diagonal terms are non-zero. Defaults to False.

Type:

bool

legs

Tensor legs

Type:

Tuple[int]

name

Tensor Operator name. Defaults to None.

Type:

str

Tensor Example:

\[\begin{split}\begin{array}{c} j_0^\prime & & j_1^\prime & & & & j_3^\prime \\ | & & | & & & & | \\ W(0) & - & W(1) & - & - & - & W(3) \\ | & & | & & & & | \\ j_0 & & j_1 & & & & j_3 \\ \end{array}\\\end{split}\]
\[j_0 = 0, 1, .. n-1, j_1 = 0, 1, .. m-1, j_3 = 0, 1, .. l-1\]

tensor = np.array([W0, W1, W3]).

When only_diag=False, above example is shape = (n, n, m, m, l, l),

legs=(0, 0, 1, 1, 3, 3), (even index is bra, odd index is ket).

Otherwise, only_diag=True, shape=(n,m,l), legs=(0, 1, 3).

W_1[beta_01,j_1^prime,j_1,beta_13] = \(W{\small \beta_{0,1}}_{j_1}^{j_1^\prime}{\small \beta_{1,3}}\)

add_dammy_legs(add_legs: tuple[int, ...], add_shape: tuple[int, ...]) → TensorOperator[source]

Dammy legs addition

Parameters:

  • add_legs (Tuple[int]) – additional tensor legs. e.g. add \(j_1 ` is ``add_legs=(1,)`\).

  • add_shape (Tuple[int]) – additional tensor shape. e.g. add \(j_1 = 0,1,\ldots,3\) is add_shape=(4,).

Returns:

Filled Tensor

Return type:

TensorOperator

e.g. from legs = (0,1) to legs=(0,1,2), where dammy legs is 2.

bond_dimension: list[int]
decompose(bond_dimension: list[int] | int | None = None, decompose_type: str = 'SVD', rate: float | None = None, square_sum: bool = True, overwrite: bool = False) → list[ndarray][source]

MPO Decomposition

Parameters:

  • bond_dimension (List[int] or int) – MPO bond dimension

  • decompose_type (str) – Tensor Train Decomposition Type 'QRD' or 'SVD'. Defaults to 'SVD'.

  • rate (float) – SVD decomposition contribution rate truncation. Defaults to None.

  • square_sum (bool) – Whether rate is based on the sum of squares of the eigenvalues, or simply sum of the eigenvalues. Square sum can calculate easily by trace of matrix, but need high rate percentage. We recommend rate=0.99999999999, if square_sum=True.

Returns:

MPO

Return type:

List[np.ndarray]

property dtype: dtype
estimate_error(tensor: ndarray | None = None)[source]

Estimate Frobenius norm between self and given-tensor or renormalized tensor

Parameters:

tensor (np.ndarray) – given tensor (option)

Returns:

Error %

Return type:

float

full_bond_dimension: list[int]
get_bond_dimension(bond_dimension: list[int] | int | None = None)[source]

get MPO bond-dimension

Parameters:

bond_dimension (List[int] | int) – MPO bond dimension

get_frobeinus_norm() → float[source]

Get Frobenius norm of tensor operator

get_num_of_elements() → tuple[int, int | None][source]

Get number of decomposed tensor operator element

Returns:

Number of tensor elements (before, after) decomposition

Return type:

Tuple[int, int]

get_tensor_full() → ndarray[source]

fill off-diagonal term zero if only_diag=True

legs: tuple[int, ...]
classmethod load(name) → TensorOperator[source]

Load Tensor Operator Object from binary file in ./tensor_operators directory

name: str
only_diag: bool
restore_from_decoposed() → ndarray[source]

Restore tensor_orig from MPO (tensor_decomposed). This routine is almost the same as mps_cls.buid_CICoef

save(name: str | None = None) → None[source]

Save Tensor Operator Object as binary file in ./tensor_operators directory

Parameters:

name (str) – file name not including extension

shape: tuple[int, ...]
tensor_decomposed: list[ndarray]
tensor_full: ndarray
tensor_orig: ndarray
pytdscf.dvr_operator_cls.construct_fulldimensional(dvr_prims: list[DVRPrimitivesMixin], func: Callable | None = None, db: str | None = None, ref_ene: float = 0.0, dipole: bool = False, efield: tuple[float, float, float] = (0.01, 0.01, 0.01)) → dict[tuple[int, ...], TensorOperator][source]

Construct full-dimensional Operator from DVR grid-based PES

Parameters:

  • dvr_prims (List[DVRPrimitivesMixin]) – DVR functions. Sorted in Database order.

  • func (Optional,Callable) – Potential energy function

  • db (Optional,str) – Electronic structure calculation database path, such as ‘foo/hoge.db’

  • ref_ene (Optional, float) – Opt coord energy in a.u., if you need subtract from function or database energy. Defaults to 0.0

  • dipole (Optional, bool) – Get dipole moment from database. Defaults to False

  • efield (Optional, List[float, float, float]) – Electronic field for inner products with dipole moment. Defaults to [1.0,1.0,1.0]

Returns:

full-dimensional tensor operator

Return type:

Dict[Tuple[int], TensorOperator]

pytdscf.dvr_operator_cls.construct_kinetic_operator(dvr_prims: list[DVRPrimitivesMixin], coefs: list[float] | None = None, forms: str = 'mpo') → dict[tuple[tuple[int, int], ...], TensorOperator][source]

Note

The off-diagonal terms of the kinetic operator are not zero in DVR basis.

Kinetic energy operator (KEO) in linear coordinate is represented by following matrix product operator (MPO):

\[\begin{split}\begin{pmatrix} -\frac{\hat{P}_1^2}{2} & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ -\frac{\hat{P}_2^2}{2} & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ -\frac{\hat{P}_3^2}{2} & 1 \end{pmatrix} \begin{pmatrix} 1 \\ -\frac{\hat{P}_4^2}{2} \end{pmatrix}\end{split}\]
Parameters:

  • dvr_prims (List[DVRPrimitivesMixin]) – DVR functions

  • coefs (List[float]) – Coefficients of kinetic operator. Defaults to [1.0, 1.0, ...], i.e., \(\sum_i -\frac{1}{2}\frac{d^2}{dQ_i^2}\) is given. For example, when one uses dimensionless coordinate and rotational coordinate, the kinetic operator is given by \(\hat{T} = -\frac{\omega}{2} \frac{d^2}{dx^2} - \frac{1}{2I} \frac{d^2}{d\theta^2}\),then coefs should be given by \([\omega, 1/I]\).

  • forms (str) – Either ‘sop’ or ‘mpo’. Defaults to ‘mpo’.

pytdscf.dvr_operator_cls.construct_nMR_recursive(dvr_prims: list[DVRPrimitivesMixin], nMR: int = 3, ndof: int | None = None, func: dict[tuple[int, ...], Callable] | None = None, db: str | None = None, df: DataFrame | None = None, active_dofs: list[int] | None = None, site_order: dict[int, int] | None = None, zero_indices: list[int] | None = None, return_tensor: bool = False, include_const_in_mpo: bool = False, ref_ene: float | None = None, dipole: bool = False, efield: tuple[float, float, float] = (1.0, 1.0, 1.0), rate: float = 1.0, k: int = 200, nsweep: int = 4) → list[ndarray] | dict[tuple[int, ...], TensorOperator][source]

Construct n-Mode Representation Operator

n-Mode Representation (nMR) are used for reducing grid points for ab initio calculation. To avoid duplicate addition of the same coordinates, use principle of inclusion and exclusion. (It is easy to understand by Venn-diagram)

Parameters:

  • dvr_prims (List[DVRPrimitivesMixin]) – DVR functions

  • nMR (Optional, int) – Mode Representation Number. Defaults to 3.

  • ndof (Optional, int) – number of mode including in database.

  • func (Optional, Dict[Tuple[int], Callable]) – E.g. Potential energy function

  • db (Optional, str) – Electronic structure calculation database path, such as ‘foo/hoge.db’

  • df (Optional, pl.DataFrame) – Electronic structure calculation dataframe.

  • active_dofs (Optional, List[int]) – Active modes Defaults to ALL.

  • site_order (Optional, Dict[int, int]) – MPO site order. Defaults to DOF index order.

  • zero_indices (Optional, List[int]) – nMR criteria coordinate grid.

  • return_tensor (bool) – Return before decomposition. Defaults to False.

  • include_const_in_mpo (bool) – Include scalar constant value (such as refernce energy) in MPO. Defaults to False.

  • ref_ene (float) – Opt coord energy in a.u., if you need subtract from function or database energy. Defaults to 0.0

  • dipole (bool) – Get dipole moment from database. Defaults to False

  • efield (List[float]) –

    Electronic field for inner products with dipole moment.Defaults to [1.0e-02 , 1.0e-02, 1.0e-02]

  • rate (float) – SVD contribution rate. Defaults to 1.0

  • k (int) – SVD bond dimension. Defaults to 200

Returns:

nMR MPO

Return type:

List[np.ndarray]

Example:

\[\begin{split}V &:= \sum_{\boldsymbol{x}\in \mathrm{coord.}} V(x_\alpha,x_\beta) \\ &= V_0 + V_{\mathrm{1MR}} + V_{\mathrm{2MR}}\\\end{split}\]

where,

\[\begin{split}V_0 &= V(0,0)\\ V_{\mathrm{1MR}} &= \sum_{x_\alpha} \left(V(x_\alpha,0) - V_0\right) + \sum_{x_\beta} \left(V(0, x_\beta) - V_0\right)\\ V_{\mathrm{2MR}} &= \sum_{x_\alpha}\sum_{x_\beta} \ \left(V(x_\alpha,x_\beta) - V(x_\alpha,0)- V(0, x_\beta) - V_0\right)\end{split}\]
pytdscf.dvr_operator_cls.database_to_dataframe(db: str, reference_id: int | None = None, reference_grids: tuple[int, ...] | None = None) → DataFrame[source]

ASE Database to Polars DataFrame

DataFrame takes more memory than Database, but it is faster than Database.

Parameters:

  • db (str) – ASE Database path such as “foo/hoge.db”. Extension must be “.db”

  • reference_id (int) – Database id of reference geometry

  • reference_grids (Optional, List[int]) – Grid indices of reference geometry

Returns:

Polars DataFrame

Return type:

pl.DataFrame

  • Dipole is saved in Debye.

  • Energy is saved in Hartree.

pytdscf.dvr_operator_cls.deepcopy(item)[source]

copy.deepcopy() is too lazy

pytdscf.dvr_operator_cls.tensor_dict_to_mpo(tensor_dict: dict[tuple[int, ...], ndarray], rate: float = 1.0, nsweep: int = 4) → list[ndarray][source]

Convert tensor dictionary to MPO

Parameters:

  • tensor_dict (Dict[Tuple[int, ...], np.ndarray]) – Tensor dictionary

  • rate (float) – SVD contribution rate. Defaults to 1.0. If rate < 1.0, MPO bond dimension is reduced. Typically, rate = 0.999999999 is enough.

  • nsweep (int) – Number of sweep in SVD compression. Defaults to 4

Returns:

MPO

Return type:

List[np.ndarray]

Notes

scalar term tensor_dict[()] is not included in MPO.

pytdscf.dvr_operator_cls.tensor_dict_to_tensor_op(tensor_dict: dict[tuple[int, ...], ndarray]) → dict[tuple[int, ...], TensorOperator][source]

Convert tensor dictionary to TensorOperator dictionary

Parameters:

tensor_dict (Dict[Tuple[int, ...], np.ndarray]) – Tensor dictionary

Returns:

TensorOperator dictionary

Return type:

Dict[Tuple[int, …], TensorOperator]

pytdscf.hamiltonian_cls module

The operator modules consists Hamiltonian.

class pytdscf.hamiltonian_cls.HamiltonianMixin(name: str, nstate: int, ndof: int, matJ: list[list[complex | float]] | None = None)[source]

Bases: object

Hamiltonian abstract class.

name

The name of operator.

Type:

str

onesite_name

The name of onesite term. Defaults tp 'onesite'.

Type:

str

nstete

The number of electronic states.

Type:

int

ndof

The degree of freedoms.

Type:

int

coupleJ

The couplings of electronic states. coupleJ[i][j] denotes the couplings between i-states (bra) and j-states (ket). This contains scalar operator.

Type:

List[List[complex]]

class pytdscf.hamiltonian_cls.PolynomialHamiltonian(ndof: int, nstate: int = 1, name: str = 'hamiltonian', matJ=None)[source]

Bases: HamiltonianMixin

PolynomialHamiltonian package, contains some model Hamiltonian. Sometimes, this class also manage observable operators, such as dipole operator

coupleJ

The couplings of electronic states. coupleJ[i][j] denotes the couplings between i-states (bra) and j-states (ket). This contains scalar operator.

Type:

List[List[complex]]

onesite

The onesite operators. onesite[i][j][k] denotes k-th onesite operator (TermOneSiteForm) between i-states (bra) and j-states(ket).

Type:

List[List[List[TermProductForm or TermOneSiteForm]]]

general

The multisite operators. The definition is almost the same as onesite.

Type:

List[List[List[TermProductForm or TermOneSiteForm]]]

is_hermitian()[source]
set_ConIns_potential(basinfo)[source]

Intra-state terms i><i

set_HO_potential(basinfo, *, enable_onesite=True) → None[source]

Setting Harmonic Oscillator polynomial-based potential

set_HO_potential_ham1(basinfo)[source]
set_LVC(bas_info, first_order_coupling: dict[tuple[int, int], dict[int, float]])[source]

Setting Linear Vibronic Coupling (LVC) Model

Parameters:

  • bas_info (BasisInfo) – Basis information

  • first_order_coupling (Dict[Tuple[int, int], Dict[int, float]]) – First order coupling if {(0,1): {2: 0.1}} is given, then the coupling between 0-state and 1-state is 0.1*Q_2

set_henon_heiles(omega: float, lam: float, f: int, omega_unit: str = 'cm-1', lam_unit: str = 'a.u.') → list[list[TermProductForm]][source]

Setting Henon-Heiles potential

In dimensionless form, the Hamiltonian is given by

\[\hat{H} = \frac{\omega}{2}\sum_{i=1}^{f} \left( - \frac{\partial^2}{\partial q_i^2} + q_i^2 \right) \ + \lambda \left( \sum_{i=1}^{f-1} q_i^2 q_{i+1} - \frac{1}{3} q_{i+1}^3 \right)\]

But PyTDSCF adopts mass-weighted coordinate, thus the Hamiltonian is given by

\[\hat{H} = \frac{1}{2}\sum_{i=1}^{f} \left( - \frac{\partial^2}{\partial Q_i^2} + \omega^2 Q_i^2 \right) \ + \lambda \omega^{\frac{3}{2}} \left( \sum_{i=1}^{f-1} Q_i^2 Q_{i+1} - \frac{1}{3} Q_{i+1}^3 \right)\]
Parameters:

  • omega (float) – Harmonic frequency in a.u.

  • lam (float) – Coupling constant in a.u.

  • f (int) – Number of degrees of freedom

Returns:

List of Hamiltonian terms

Return type:

List[List[TermProductForm]]

set_henon_heiles_2D_4th(lam: float = 0.2) → list[list[TermProductForm]][source]
\[H(x,y) = -\frac{1}{2}\left(\frac{d^2}{dx^2} + \frac{d^2}{dy^2}\right) \ + \frac{1}{2}(x^2 + y^2) + \lambda \left(xy^2 - \frac{1}{3}x^3\right)\ + \lambda^2 \left(\frac{1}{16}(x^4 + y^4) + \frac{1}{8}x^2y^2\right)\]
class pytdscf.hamiltonian_cls.TensorHamiltonian(ndof: int, potential: list[list[dict[tuple[int | tuple[int, int], ...], TensorOperator]]], name: str = 'hamiltonian', kinetic: list[list[dict[tuple[tuple[int, int], ...], TensorOperator] | None]] | None = None, decompose_type: Literal['QRD', 'SVD'] = 'QRD', rate: float | None = None, bond_dimension: list[int] | int | None = None, backend: Literal['jax', 'numpy'] = 'jax')[source]

Bases: HamiltonianMixin

Hamiltonian in tensor formulation.

mpo

MPO between i-state and j-state

Type:

List[List[MatrixProductOperators]]

distribute_mpo_cores()[source]
mpo: list[list[MatrixProductOperators | None]]
class pytdscf.hamiltonian_cls.TermOneSiteForm(coef: float | complex, op_dof: int, op_key: str)[source]

Bases: object

Operator with legs in only one degree of freedom

coef

The coefficient of the product operator, e.g. 1.0.

Type:

float or complex

op_dofs

The MPS lattice (0-)index of operator, e.g. 0.

Type:

int

op_keys

The type of operators for each legs, e.g. 'q', 'd^2'.

Type:

str

class pytdscf.hamiltonian_cls.TermProductForm(coef: float | complex, op_dofs: list[int], op_keys: list[str])[source]

Bases: object

Hold a product form of operators, such as \(1.0 q_1 q_2^2\).

coef

The coefficient of the product operator, e.g. 1.0.

Type:

float or complex

op_dofs

The MPS lattice (0-)index of operator, e.g. [0, 1].

Type:

List[int]

op_keys

The type of operators for each legs, e.g. ['q', 'q^2'].

Type:

List[str]

mode_ops

The pair of op_dofs and op_keys.

Type:

Dict[Tuple[int], str]

blockop_key_sites

Operator legs divided by central site. E.g. 3 DOFs, 'q_0q^2_1' gives {left:['ovlp','q_0', 'q_0q^2_1'], centr:['q_0', 'q^2_1', 'ovlp'], right:['q^2_1', 'ovlp', 'ovlp']}.

Type:

Dict[str, List[str]]

static convert_key(op_dofs: list[int], op_keys: list[str]) → str[source]

convert product operator to string key

Parameters:

  • op_dofs (List[int]) – The MPS lattice (0-)index of operator, e.g. [0, 1].

  • op_keys (List[str]) – The type of operators for each legs, e.g. ['q', 'q^2'].

Returns:

string key, e.g. q_0 q^2_1.

Return type:

str

Examples

>>> TermProductForm.convert_key([1, 2], ['q', 'q^2'])
    'q_1 q^2_2'
is_op_ovlp(block_lcr_list: list[str], psite: int) → bool[source]

check MPS lattice block is 'ovlp' operator or not.

Parameters:

  • block_lcr_list (List[str]) – 'left', or 'centr', or 'right'

  • psite (int) – (0-) index of center site. LCR of C.

Returns:

Return True if all operator types with legs in the given block are 'ovlp' when regard centered site as the psite-th site.

Return type:

bool

set_blockop_key(ndof: int, *, print_out: bool = False)[source]

For the complementary operator, the product operator is set to an attribute blockop_key_sites that is classified according to the center site.

Parameters:

  • ndof (int) – The length of MPS lattice, that is degree of freedoms.

  • print_out (bool) – Print blockop_key_sites or not. Defaults to False.

term_key() → str[source]

call convert_key attributes in this class.

Returns:

string key, e.g. q_1 q^2_2.

Return type:

str

pytdscf.hamiltonian_cls.read_potential_nMR(potential_emu: dict[tuple[int, ...], float | complex], *, active_modes: list[int] | None = None, name: str = 'hamiltonian', cut_off: float | None = None, dipole_emu: dict[tuple[int, ...], tuple[float, float, float]] | None = None, print_out: bool = False, active_momentum=None, div_factorial: bool = True, efield: tuple[float, float, float] = (1.0, 1.0, 1.0)) → PolynomialHamiltonian[source]

Construct polynomial potential, such as n-Mode Representation form. multi-state nMR is not yet implemented.

Parameters:

  • potential_emu (Dict[Tuple[int], float or complex]) –

    The coeffcients of operators. The index is start from 1. Factorial factor is not included in the coefficients. E.g. V(q_1) = 0.123 * 1/(2!) * q_1^2 gives {(1,1):0.123}.

    units are all a.u.:

    k_orig{ij} in [Hartree/emu/bohr**2]

    k_orig{ijk} in [Hartree/emu**i(3/2)/BOHR**3]

    k_orig{ijkl} in [Hartree/emu**(2)/BOHR**4]

  • active_modes (List[int]) – The acctive degree of freedoms. Note that, The order of the MPS lattice follows this list. The index is start from 1.

  • name (str) – The name of operator.

  • cut_off (float) – The threshold of truncatig terms by coeffcients after factorial division.

  • dipole_emu (Dict[Tuple[int], float or complex]) –

    Dipole moment operator to generate a vibrational excited states. Defaults to None. The definition is almost the same as potential_emu, but this value is 3d vector (list).

    mu{ij} in [Debye/emu/bohr**2]

    mu{ijk} in [Debye/emu**(3/2)/BOHR**3]

    mu{ijkl} in [Debye/emu**(2)/BOHR**4]

  • print_out (bool) – Defaults to False

  • div_factorial (bool) – Whether or not divided by factorial term. Defaults to True.

  • active_momentum (List[Tuple[int, float]]) – [(DOF, coef)]. Defaults to None When active_momentum are set, unselected kinetic energy operators are excluded. Default option include all kinetic energy operators.

  • efield (List[float]) – The electric field in [Hartree Bohr^-1 e^-1]. Defaults to [1.0, 1.0, 1.0]

Returns:

operator

Return type:

PolynomialHamiltonian

pytdscf.hamiltonian_cls.truncate_terms(terms: list[TermProductForm], cut_off: float = 0.0) → list[TermProductForm][source]

Truncate operators by a certain cut_off.

Parameters:

  • terms (List[TermProductForm]) – Terms list of sum of products form before truncation.

  • cut_off (float) – The threshold of truncation. unit is a.u.

Returns:

Terms list of sum of products form after truncation.

Return type:

terms (List[TermProductForm])

pytdscf.model_cls module

The model class for the Time-Dependent Schrodinger Equation (TDSE) calculation.

class pytdscf.model_cls.BasInfo(prim_info, spf_info=None, ndof_per_sites=None)[source]

Bases: object

The Wave function basis information class

Parameters:

  • prim_info (List[List[PrimBas_HO or DVRPrimitivesMixin]]) – prim_info[istate][idof] gives PrimBas_HO.

  • spf_info (List[List[int]]) – spf_info[istate][idof] gives the number of SPF.

  • ndof_per_sites (bool) – Defaults to None. N.Y.I.

prim_info

prim_info[istate][idof] gives PrimBas_HO.

Type:

List[List[PrimBas_HO or DVRPrimitivesMixin]]

spf_info

spf_info[istate][idof] gives the number of SPF.

Type:

List[List[int]]

nstate

The number of electronic states.

Type:

int

ndof

The degree of freedoms.

Type:

int

nspf

The number of SPF.

Type:

int

nspf_list

The number of SPFs.

Type:

List[int]

nspf_rangelist

The indices of SPFs.

Type:

List[List[int]]

get_ndof() → int[source]

Get ndof attributes

Returns:

Degree of freedom

Return type:

int

get_ndof_per_sites() → list[int][source]

Get ndof_per_sites attributes

get_ngrid(istate: int, idof: int) → int[source]
get_nprim(istate: int, idof: int) → int[source]

Get number of primitive basis

Parameters:

  • istate (int) – index of electronic states

  • idof (int) – index of degree of freedom

Returns:

Number of primitive basis

Return type:

int

get_nspf(istate: int, idof: int) → int[source]

Get number of SPF

Parameters:

  • istate (int) – index of electronic states

  • idof (int) – index of degree of freedom

Returns:

Number of SPF

Return type:

int

get_nspf_list(istate: int) → list[int][source]

Get number of SPFs list nspf_list attributes

Parameters:

istate (int) – index of electronic states

Returns:

Number of SPFs. e.g. [2, 2]

Return type:

list(int)

get_nspf_rangelist(istate: int) → list[list[int]][source]

Get number of SPFs list nspf_rangelist attributes

Parameters:

istate (int) – index of electronic states

Returns:

the indices of SPFs. e.g. [[0,1,2],[0,1,2]]

Return type:

List[List[int]]

get_nstate() → int[source]

Get nstate attributes

Returns:

Number of electronic states

Return type:

int

get_primbas(istate: int, idof: int) → PrimBas_HO[source]

Get prim_info[istate][idof] attributes :param istate: index of electronic states :type istate: int :param idof: index of degree of freedom :type idof: int

Returns:

The primitive basis of istate, idof.

Return type:

PrimBas_HO

class pytdscf.model_cls.Model(basinfo: BasInfo, operators: dict[str, HamiltonianMixin], *, build_td_hamiltonian: PolynomialHamiltonian | None = None, space: Literal['hilbert', 'liouville'] = 'hilbert')[source]

Bases: object

The wavefunction and operator information class

Parameters:

  • basinfo (BasInfo) – The wavefunction basis information

  • operators (Dict) – The operators. operators[name] gives a operator named name, such as hamiltonian.

  • build_td_hamiltonian (PolynomialHamiltonian) – Time dependent PolynomialHamiltonian. Defaults to None.

init_weight_VIBSTATE

Initial weight of VIBSTATE. List length is [nstate, ndof].

Type:

List[List[float]]

init_weight_GRID

Initial weight of GRID. List length is [nstate, ndof].

Type:

List[List[float]]

init_weight_ESTATE

Initial weight of ESTATE. List length is nstate.

Type:

List[float]

ints_prim_file

The file name of the primitive integrals

Type:

str

m_aux_max

The maximum number of auxiliary basis.

Type:

int

basinfo

The wavefunction basis information

Type:

BasInfo

hamiltonian

PolynomialHamiltonian

Type:

PolynomialHamiltonian

observables

Observable operators, such as PolynomialHamiltonian, dipole moment, occupation number etc.

Type:

Dict

build_td_hamiltonian

Time-dependent PolynomialHamiltonian.

Type:

PolynomialHamiltonian

get_ndof() → int[source]
Returns:

Degree of freedoms

Return type:

int

get_ndof_per_sites()[source]

N.Y.I ?

get_nprim(istate: int, idof: int) → int[source]
Parameters:

  • istate (int) – index of electronic states

  • idof (int) – index of degree of freedom

Returns:

The number of primitive basis in istate, idof.

Return type:

int

get_nspf(istate: int, idof: int) → int[source]
Parameters:

  • istate (int) – index of electronic states

  • idof (int) – index of degree of freedom

Returns:

The number of SPF in istate, idof.

Return type:

int

get_nspf_list(istate: int) → list[int][source]
Parameters:

istate (int) – index of electronic states

Returns:

Number of SPFs. e.g. [2, 2]

Return type:

List(int)

get_nspf_rangelist(istate: int) → list[list[int]][source]
Parameters:

istate (int) – index of electronic states

Returns:

the indices of SPFs. e.g. [[0,1,2],[0,1,2]]

Return type:

List[List[int]]

get_nstate() → int[source]
Returns:

Number of electronic states

Return type:

int

get_primbas(istate: int, idof: int) → PrimBas_HO[source]
Parameters:

  • istate (int) – index of electronic states

  • idof (int) – index of degree of freedom

Returns:

The primitive basis in istate, idof.

Return type:

primints_cls.PrimBas_HO

init_HartreeProduct: list[list[list[float]]] | None = None
init_weight_ESTATE: list[float] | None = None
init_weight_VIBSTATE: list[list[float]] | None = None
ints_prim_file: str | None = None
m_aux_max: int | None = None

pytdscf.properties module

Property handling module

class pytdscf.properties.Properties(wf: WFunc, model: Model, time=0.0, t2_trick=True, wf_init=None, reduced_density=None)[source]

Bases: object

Structure class of calculated property

wf

wave function

Type:

WFunc

model

model

Type:

Model

time

time in atomic units

Type:

float

t2_trick

whether to use so-called T/2 trick

Type:

bool

autocorr

auto-correlation function

Type:

complex

energy

energy

Type:

float

pops

populations of each state

Type:

List[float]

expectations

observables of given operator

Type:

List[float]

export_properties(*, autocorr_per_step=1, populations_per_step=1, observables_per_step=1, bonddim_per_step=1)[source]
get_properties(*, autocorr=True, energy=True, norm=True, populations=True, observables=True, bonddim=True, autocorr_per_step=1, energy_per_step=1, norm_per_step=1, populations_per_step=1, observables_per_step=1, bonddim_per_step=1)[source]
update(delta_t)[source]

update time

pytdscf.simulator_cls module

The main simulator module of PyTDSCF

This module consists of Simulator class.

class pytdscf.simulator_cls.Simulator(jobname: str, model: Model, ci_type: Literal['mps', 'mctdh', 'ci'] = 'mps', backend: Literal['jax', 'numpy'] = 'jax', proj_gs: bool = False, t2_trick: bool = True, verbose: int = 2)[source]

Bases: object

The simulator of the PyTDSCF

set parameter of the restart, propagate, operate dipole, save_file etc …

Parameters:

  • jobname (str) – the jobname

  • model (model_cls.Model) – run parameter (basis, hamiltonian, observable, bond dimension, initial_weight etc.)

  • t2_trick (bool, optional) – Use so-called t/2 trick in auto-correlation. Note that it requires initial state to be real. Defaults to True.

  • ci_type (str, optional)

  • backend (str, optional) – JAX or Numpy. Defaults to 'jax'. When polynomial operator, FBR basis and small bond-dimension is used, 'Numpy' is recommended.

  • proj_gs (bool, optional) – Initial state is projected from the ground state. Defaults to False. If proj_gs=True, one must be set attribute model.primbas_gs: List[Primbas_HO].

backup_interval: int
get_initial_wavefunction(ints_prim: PrimInts) → WFunc[source]
get_primitive_integrals() → PrimInts[source]
maxstep: int
operate(maxstep: int = 10, restart: bool = False, savefile_ext: str = '_operate', loadfile_ext: str = '_gs', verbose: int = 2) → tuple[float, WFunc][source]

Apply operator such as dipole operator to the wavefunction

Parameters:

  • maxstep (int, optional) – Maximum number of iteration. Defaults to 10.

  • restart (bool, optional) – Restart from the previous wavefunction. Defaults to False.

  • backend (str, optional) – JAX or Numpy. Defaults to 'jax'.

  • savefile_ext (str, optional) – Extension of the save file. Defaults to '_operate'.

  • loadfile_ext (str, optional) – Extension of the load file. Defaults to '_gs'. When restart=False, loadfile_ext is ignored.

  • verbose (int, optional) – Verbosity level. Defaults to 2.

Returns:

norm of O|Ψ> and Wavefunction after operation.

Return type:

Tuple[float, WFunc]

propagate(stepsize: float = 0.1, maxstep: int = 5000, restart: bool = False, savefile_ext: str = '', loadfile_ext: str = '_operate', backup_interval: int = 1000, autocorr: bool = True, energy: bool = True, norm: bool = True, populations: bool = True, observables: bool = False, reduced_density: tuple[list[tuple[int, ...]], int] | None = None, Δt: float | None = None, thresh_sil: float = 1e-09, autocorr_per_step: int = 1, observables_per_step: int = 1, energy_per_step: int = 1, norm_per_step: int = 1, populations_per_step: int = 1, parallel_split_indices: list[tuple[int, int]] | None = None, adaptive: bool = False, adaptive_Dmax: int = 20, adaptive_dD: int = 5, adaptive_p_proj: float = 0.0001, adaptive_p_svd: float = 1e-07, integrator: Literal['lanczos', 'arnoldi'] = 'lanczos', step_size_is_fs: bool = True) → tuple[float, WFunc][source]

Propagation

Parameters:

  • stepsize (float, optional) – Step size in fs. Defaults to 0.1.

  • maxstep (int, optional) – Maximum number of steps. Defaults to 5000., i.e. 500 fs.

  • restart (bool, optional) – Restart from the previous wavefunction. Defaults to False.

  • savefile_ext (str, optional) – Extension of the save file. Defaults to ''.

  • loadfile_ext (str, optional) – Extension of the load file. Defaults to '_operate'. When restart=False, loadfile_ext is ignored.

  • backup_interval (int, optional) – Number of steps at which, the wavefunction is saved. Defaults to 1000.

  • autocorr (bool, optional) – Calculate autocorrelation function. Defaults to True.

  • energy (bool, optional) – Calculate energy. Defaults to True.

  • norm (bool, optional) – Calculate norm. Defaults to True.

  • populations (bool, optional) – Calculate populations. Defaults to True.

  • observables (bool, optional) – Calculate observables. Defaults to False.

  • reduced_density (Dict[Tuple[int, ...], int], optional) – Calculate reduced density of the given modes. For example, ([(0, 1),], 10) means calculate the diagonal elements of reduced density of the \(\rho_{01}=\mathrm{Tr}_{p\notin \{0,1\}}\left|\Psi^{(\alpha)}\rangle\langle\Psi^(\alpha)\right|\) in per 10 steps. Note that it requires enough disk space. Defaults to None. It is better if the target modes are close to rightmost, i.e., 0. (Because this program calculate property in the most right-canonized form of MPS.) If you want coherence, i.e., off-diagonal elements of density matrix, you need to set like ([(0, 0), ], 10).

  • Δt (float, optional) – Same as stepsize

  • thresh_sil (float) – Convergence threshold of short iterative Lanczos. Defaults to 1.e-09.

  • step_size_is_fs (bool, optional) – If True, stepsize is in fs. Defaults to True.

Returns:

Energy during propagation (it conserves) and Wavefunction after propagation.

Return type:

Tuple[float, WFunc]

relax(stepsize: float = 0.1, maxstep: int = 20, improved: bool = True, restart: bool = False, savefile_ext: str = '_gs', loadfile_ext: str = '', backup_interval: int = 10, norm: bool = True, populations: bool = True, observables: bool = False, integrator: Literal['lanczos', 'arnoldi'] = 'lanczos') → tuple[float, WFunc][source]

Relaxation

Parameters:

  • stepsize (float, optional) – Step size in fs. Defaults to 0.1. This is used only when imaginary time propagation is used.

  • maxstep (int, optional) – Maximum number of steps. Defaults to 20.

  • improved (bool, optional) – Use improved relaxation. Defaults to True.

  • restart (bool, optional) – Restart from the previous wavefunction. Defaults to False.

  • savefile_ext (str, optional) – Extension of the save file. Defaults to '_gs'.

  • loadfile_ext (str, optional) – Extension of the load file. Defaults to ''. When restart=False, loadfile_ext is ignored.

  • backup_interval (int, optional) – Number of steps at which, the wavefunction is saved. Defaults to 10.

  • norm (bool, optional) – Calculate norm. Defaults to True.

  • populations (bool, optional) – Calculate populations. Defaults to True.

  • observables (bool, optional) – Calculate observables. Defaults to False.

Returns:

Energy after relaxation in Eh, and Wavefunction after relaxation.

Return type:

Tuple[float, WFunc]

save_wavefunction(wf: WFunc, log: bool = False)[source]
stepsize: float

pytdscf.spectra module

Plot spectra from FFT of auto-correlation dat file

pytdscf.spectra.export_spectrum(wave_number, intensity, filename: str = 'spectrum.dat')[source]

Export spectrum to dat file

Parameters:

  • wave_number (np.ndarray) – wave number in cm-1

  • intensity (np.ndarray) – intensity

  • filename (str, optional) – path to dat file. Defaults to “spectrum.dat”.

pytdscf.spectra.ifft_autocorr(time_fs: ndarray, autocorr: ndarray, E_shift: float = 0.0, window: str = 'cos2', power: bool = False) → tuple[ndarray, ndarray][source]

Inverse FFT of auto-correlation data

Parameters:

  • time_fs (np.ndarray) – time in fs

  • autocorr (np.ndarray) – auto-correlation data

  • E_shift (float, optional) – energy shift in eV. we often use ZPE. Defaults to 0.0.

  • window (str, optional) – window function. Defaults to “cos2”.

  • power (bool, optional) – if True, intensity is power spectrum. Defaults to False, thus IR absorption spectrum in arb. unit. and energy is shifted.

Returns:

(wave number in cm-1, intensity)

Return type:

Tuple[np.ndarray, np.ndarray]

pytdscf.spectra.load_autocorr(dat_file: str) → tuple[ndarray, ndarray][source]

Load auto-correlation data from dat file

Parameters:

dat_file (str) – path to dat file

Returns:

(time in fs, auto-correlation data)

Return type:

Tuple[np.ndarray, np.ndarray]

pytdscf.spectra.plot_autocorr(time_fs: ndarray, autocorr: ndarray, gui: bool = True)[source]

Plot auto-correlation data

Parameters:

  • time_fs (np.ndarray) – time in fs

  • autocorr (np.ndarray) – auto-correlation data

  • gui (bool, optional) – if True, plt.show() use GUI. Defaults to True.

pytdscf.spectra.plot_spectrum(wave_number: ndarray, intensity: ndarray, lower_bound: float = 0.0, upper_bound: float = 4000.0, filename: str = 'spectrum.pdf', export: bool = True, show_in_eV: bool = False, gui: bool = True, normalize: bool = True)[source]

Plot spectrum

Parameters:

  • wave_number (np.ndarray) – wave number in cm-1

  • intensity (np.ndarray) – intensity

  • lower_bound (float, optional) – lower bound of wave number in cm-1. Defaults to 0.0.

  • upper_bound (float, optional) – upper bound of wave number in cm-1. Defaults to 4000.0.

  • filename (str) – path to figure. Defaults to “spectrum.pdf”.

  • export (bool) – export spectrum to dat file. Defaults to True.

  • show_in_eV (bool) – show in eV. Defaults to False.

  • gui (bool) – if True, plt.show() use GUI. Defaults to True.

  • normalize (bool) – normalize intensity in range of between lower and upper bound. Defaults to True.

pytdscf.units module

Unit Convert Class

xxx_in_yyy = z

means

1 [xxx] = z [yyy]

Example

>>> from pytdscf.units import au_in_angstrom
>>> print(au_in_angstrom)
0.529177210903

Units is based on SciPy constants (CODATA 2018)

pytdscf.units.au_in_cm1

1 a.u. in cm^-1

Type:

float

pytdscf.units.au_in_fs

1 a.u. in fs

Type:

float

pytdscf.units.au_in_eV

1 a.u. in eV

Type:

float

pytdscf.units.au_in_dalton

1 a.u. in Dalton

Type:

float

pytdscf.units.au_in_angstrom

1 a.u. in Angstrom

Type:

float

pytdscf.units.au_in_debye

1 a.u. in Debye

Type:

float

pytdscf.wavefunction module

Wave function handling module

class pytdscf.wavefunction.WFunc(ci_coef: MPSCoef, spf_coef: SPFCoef, ints_prim: PrimInts)[source]

Bases: object

Wave Function (WF) class

This class holds WF and some operation to WF.

ci_coef

Decomposed A-vector object. (In MCTDH level, this object class is ci_cls.CICoef, not mps_cls.MPSCoef.)

Type:

mps_cls.MPSCoef

spf_coef

primitive coefficient of SPFs.

Type:

spf_cls.SPFCoef

ints_prim

the integral matrix element of primitives. <χ_i|V|χ_j> is time-independent.

Type:

primints_cls.PrimInts

ints_spf

the integral matrix element of SPFs. <φ_i|V|φ_j> = Σ_k Σ_l c^*_k c_l <χ_k|V|χ_l> is time-dependent.

Type:

spf_cls.SPFInts

apply_dipole(matOp) → float[source]

Get excited states by operating dipole operator to WF.

Parameters:

matOp (hamiltonian_cls.HamiltonianMixin) – Dipole operator

autocorr()[source]

get auto-correlation <Psi(t/2)*|Psi(t/2)> of WF

Returns:

auto-correlation value

Return type:

complex

bonddim() → list[int] | None[source]

get bond dimension of WF

Returns:

bond dimension of WF

Return type:

List[int]

expectation(matOp: HamiltonianMixin)[source]

get expectation value of operator

Parameters:

matOp (hamiltonian_cls.HamiltonianMixin) – operator

Returns:

expectation value (real number)

Return type:

float

get_reduced_densities(remain_nleg: tuple[int, ...]) → list[ndarray][source]

Calculate reduced density matrix of given degree of freedom pair. If (0,1) is given, calculate reduced density matrix of 0th and 1st degree of freedom.

Parameters:

remain_nleg (Tuple[int, ...]) – degree of freedom pair

Returns:

reduced density matrix of given degree of freedom pair for each state.

Return type:

List[np.ndarray]

ints_spf: SPFInts | None
norm()[source]

get WF norm

Returns:

norm of WF

Return type:

float

pop_states()[source]

get population of electronic states

Returns:

[norm of each electronic states]

Return type:

List[float]

propagate(matH: HamiltonianMixin, stepsize: float)[source]

propagate WF with VMF

Parameters:
Returns:

(g value, spf_occ)

where spf_occ = [mean field operator(mfop) overlap of each states]

Return type:

Tuple[float, List]

propagate_CMF(matH: HamiltonianMixin, stepsize_guess: float)[source]

propagate WF with CMF (Constant Mean Field)

Parameters:
Returns:

(g value, spf_occ, actual stepsize, next stepsize guess),

where spf_occ = [mean field operator(mfop) overlap of each states]

Return type:

tuple[float, list, float, float]

propagate_SM(matH: HamiltonianMixin, stepsize: float, istep: int)[source]

propagate Standard Method Wave function (SPF == PRIM)

Parameters:

  • matH (hamiltonian_cls.HamiltonianMixin) – Polynomial SOP Hamiltonian or grid MPO Hamiltonian

  • stepsize (float) – the width of time step [a.u.]

  • calc_spf_occ (Optional[bool]) – Calculate spf_occ

Returns:

spf_occ,

where spf_occ = [mean field operator(mfop) overlap of each states]

Return type:

List[float] or None

Module contents

class pytdscf.BasInfo(prim_info, spf_info=None, ndof_per_sites=None)[source]

Bases: object

The Wave function basis information class

Parameters:

  • prim_info (List[List[PrimBas_HO or DVRPrimitivesMixin]]) – prim_info[istate][idof] gives PrimBas_HO.

  • spf_info (List[List[int]]) – spf_info[istate][idof] gives the number of SPF.

  • ndof_per_sites (bool) – Defaults to None. N.Y.I.

prim_info

prim_info[istate][idof] gives PrimBas_HO.

Type:

List[List[PrimBas_HO or DVRPrimitivesMixin]]

spf_info

spf_info[istate][idof] gives the number of SPF.

Type:

List[List[int]]

nstate

The number of electronic states.

Type:

int

ndof

The degree of freedoms.

Type:

int

nspf

The number of SPF.

Type:

int

nspf_list

The number of SPFs.

Type:

List[int]

nspf_rangelist

The indices of SPFs.

Type:

List[List[int]]

get_ndof() → int[source]

Get ndof attributes

Returns:

Degree of freedom

Return type:

int

get_ndof_per_sites() → list[int][source]

Get ndof_per_sites attributes

get_ngrid(istate: int, idof: int) → int[source]
get_nprim(istate: int, idof: int) → int[source]

Get number of primitive basis

Parameters:

  • istate (int) – index of electronic states

  • idof (int) – index of degree of freedom

Returns:

Number of primitive basis

Return type:

int

get_nspf(istate: int, idof: int) → int[source]

Get number of SPF

Parameters:

  • istate (int) – index of electronic states

  • idof (int) – index of degree of freedom

Returns:

Number of SPF

Return type:

int

get_nspf_list(istate: int) → list[int][source]

Get number of SPFs list nspf_list attributes

Parameters:

istate (int) – index of electronic states

Returns:

Number of SPFs. e.g. [2, 2]

Return type:

list(int)

get_nspf_rangelist(istate: int) → list[list[int]][source]

Get number of SPFs list nspf_rangelist attributes

Parameters:

istate (int) – index of electronic states

Returns:

the indices of SPFs. e.g. [[0,1,2],[0,1,2]]

Return type:

List[List[int]]

get_nstate() → int[source]

Get nstate attributes

Returns:

Number of electronic states

Return type:

int

get_primbas(istate: int, idof: int) → PrimBas_HO[source]

Get prim_info[istate][idof] attributes :param istate: index of electronic states :type istate: int :param idof: index of degree of freedom :type idof: int

Returns:

The primitive basis of istate, idof.

Return type:

PrimBas_HO

class pytdscf.Boson(nstate: int)[source]

Bases: object

Boson basis class

Parameters:

nstate (int) – number of states

get_annihilation_matrix(margin: int = 0) → ndarray[tuple[int, ...], dtype[float64]][source]

Returns the annihilation matrix for the boson basis

\[b |n> = \sqrt{n} |n-1>\]

Example

>>> exciton = Boson(4)
>>> exciton.get_annihilation_matrix()
array([[0., 1., 0., 0.],
       [0., 0., 1.41421356, 0.],
       [0., 0., 0., 1.73205081],
       [0., 0., 0., 0.]])
>>> state = np.array([[1.0],
                      [1.0],
                      [1.0],
                      [0.0]])
>>> np.allclose(exciton.get_annihilation_matrix() @ state, np.array([[0.0], [1.0], [math.sqrt(2)], [math.sqrt(3)]]))
True
get_creation_matrix(margin: int = 0) → ndarray[tuple[int, ...], dtype[float64]][source]

Returns the creation matrix for the boson basis

\[b^\dagger |n> = \sqrt{n+1} |n+1>\]

Example

>>> boson = Boson(3)
>>> boson.get_creation_matrix()
array([[0., 0., 0.],
       [1., 0., 0.],
       [1.41421356, 0., 0.]])
>>> state = np.array([[1.0],
                      [1.0],
                      [0.0]])
>>> np.allclose(boson.get_creation_matrix() @ state, np.array([[0.0], [1.0], [math.sqrt(2)]]))
True
get_number_matrix() → ndarray[tuple[int, ...], dtype[float64]][source]

Returns the number matrix for the boson basis

\[b^\dagger b |n> = n |n>\]

Example

>>> boson = Boson(3)
>>> boson.get_number_matrix()
array([[0., 0., 0.],
       [0., 1., 0.],
       [0., 0., 2.]])
>>> state = np.array([[1.0],
                      [1.0],
                      [0.0]])
>>> np.allclose(boson.get_number_matrix() @ state, np.array([[0.0], [1.0], [0.0]]))
True
get_p2_matrix() → ndarray[tuple[int, ...], dtype[float64]][source]

Returns the p^2 matrix for the boson basis

\[p^2 = -1/2 (b^\dagger b + b b^\dagger - b^\dagger ^2 - b^2)\]
get_p_matrix(margin: int = 0) → ndarray[tuple[int, ...], dtype[complex128]][source]

Returns the p matrix for the boson basis

\[p = 1/\sqrt{2} i (b - b^\dagger) = i/\sqrt{2} (b^\dagger - b)\]
get_q2_matrix() → ndarray[tuple[int, ...], dtype[float64]][source]

Returns the p^2 matrix for the boson basis

\[p^2 = 1/2 (b^\dagger b + b b^\dagger + b^\dagger ^2 + b^2)\]
get_q_matrix(margin: int = 0) → ndarray[tuple[int, ...], dtype[float64]][source]

Returns the q matrix for the boson basis

q = 1/sqrt(2) (b + b^dagger)

names: list[str]
property nprim: int
nstate: int
class pytdscf.Exciton(nstate: int, names: list[str] | None = None)[source]

Bases: object

Exciton basis class

Parameters:

  • nstate (int) – number of states

  • names (list[str], optional) – names of the states such as [“down”, “up”]. Defaults to None, in which case the names are set to [“S0”, “S1”, …].

get_annihilation_matrix() → ndarray[tuple[int, ...], dtype[float64]][source]

Returns the annihilation matrix for the exciton basis

Example

>>> exciton = Exciton(2)
>>> exciton.get_annihilation_matrix()
array([[0., 1.],
       [0., 0.]])
>>> state = np.array([[1.0],
                      [0.0]])
>>> np.allclose(exciton.get_annihilation_matrix() @ state, np.array([[0.0], [0.0]]))
True
>>> state = np.array([[0.0],
                      [1.0]])
>>> np.allclose(exciton.get_annihilation_matrix() @ state, np.array([[1.0], [0.0]]))
get_creation_matrix() → ndarray[tuple[int, ...], dtype[float64]][source]

Returns the creation matrix for the exciton basis

Example

>>> exciton = Exciton(2)
>>> exciton.get_creation_matrix()
array([[0., 0.],
       [1., 0.]])
>>> state = np.array([[1.0],
                      [0.0]])
>>> np.allclose(exciton.get_creation_matrix() @ state, np.array([[0.0], [1.0]]))
True
>>> state = np.array([[0.0],
                      [1.0]])
>>> np.allclose(exciton.get_creation_matrix() @ state, np.array([[0.0], [0.0]]))
True
names: list[str]
property nprim: int
nstate: int
class pytdscf.Exponential(ngrid: int, length: float, x0: float = 0.0, doAnalytical: bool = True)[source]

Bases: DVRPrimitivesMixin

Exponential DVR is equivalent to FFT. It is suitable for periodic system.

See also

Primitive Function :

\[\varphi_{j}(x)= \frac{1}{\sqrt{L}} \exp \left(\mathrm{i} \frac{2 \pi j (x-x_0)}{L}\right) \quad \left(-\frac{N-1}{2} \leq j \leq \frac{N-1}{2}\right)\]

They are periodic with period \(L=x_{N}-x_{0}\):

\[\varphi_{j}(x+L)=\varphi_{j}(x)\]

Note that \(N\), i.e., ngrid must be odd number.

Note

Naturally, the DVR function \(\chi_\alpha(x)\) is given by the multiplication of delta function at equidistant grid \(\in \{x_0, x_0+\frac{L}{N}, x_0+\frac{2L}{N}, \cdots, x_0+\frac{(N-1)L}{N}\}\) point and primitive function \(\varphi_{j}(x)\) not by the Unitary transformation of the primitive function \(\varphi_{j}(x)\). (The operator \(\hat{z}\) introduced in MCTDH review Phys.Rep. 324, 1 (2000) appendix B.4.5 might be wrong.)

Parameters:

  • ngrid (int) – Number of grid. Must be odd number. (The center grid is give by index of n=ngrid//2)

  • length (float) – Length (unit depends on your potential. e.g. radian, bohr, etc.)

  • x0 (float) – left terminal point. Basis functions are defined in \([x_0, x_0+L)\).

  • doAnalytical (bool) – If True, use analytical formula. At the moment, numerical integration is not implemented.

diagnalize_pos_rep_matrix() → None[source]

Not diagonalizing the position representation matrix but manually setting the grid points and transformation matrix.

\[\chi_\alpha(x) = \sum_{j=1}^N \varphi_j(x) U_{j\alpha} = \sum_{j=1}^N \langle\varphi_j(x_\alpha)|\varphi_j(x)\rangle \varphi_j(x)\]
fbr_func(n: int, x: float) → complex[source]

Primitive functions \(\varphi_n(x)\)

Parameters:

  • n (int) – index (0 <= n < ngrid)

  • x (float) – position in a.u.

Returns:

\(\frac{1}{\sqrt{L}} \exp \left(\mathrm{i} \frac{2 \pi n (x-x_0)}{L}\right)\)

Return type:

complex

j = 0, ±1, ±2, …, ±(N-1)/2

get_1st_derivative_matrix_dvr() → ndarray[source]

Exponential 1st derivative matrix given by analytical formula

\[\begin{split}D_{\alpha \beta}^{(1), \mathrm{DVR}} = \begin{cases} 0 & \text { if } \alpha=\beta \\ \frac{\pi}{L} \frac{(-1)^{\alpha-\beta}}{\sin \left(\frac{\pi(\alpha-\beta)}{N}\right)} & \text { if } \alpha \neq \beta \end{cases} \quad (\alpha, \beta=0, \cdots, N-1)\end{split}\]
get_1st_derivative_matrix_fbr() → ndarray[source]

Numerical integral of \(\langle\varphi_j|\frac{d}{dq}|\varphi_k\rangle\)

If analytical integral is known, implemented in inheritance class

get_2nd_derivative_matrix_dvr() → ndarray[source]

Exponential 2nd derivative matrix given by analytical formula

\[\begin{split}D_{\alpha \beta}^{(2), \mathrm{DVR}} = \begin{cases} -\frac{\pi^{2}}{3 L^{2}}\left(N^{2}-1\right) & \text { if } \alpha=\beta \\ -\frac{2\pi^{2}}{L^{2}} (-1)^{\alpha-\beta} \frac{\cos\left(\frac{\pi(\alpha-\beta)}{N}\right)} {\sin ^{2}\left(\frac{\pi(\alpha-\beta)}{N}\right)} & \text { if } \alpha \neq \beta \end{cases} \quad (\alpha, \beta=0, \cdots, N-1)\end{split}\]
get_2nd_derivative_matrix_fbr() → ndarray[source]

Numerical integral of \(\langle\varphi_j|\frac{d^2}{dq^2}|\varphi_k\rangle\)

If analytical integral is known, implemented in inheritance class

get_pos_rep_matrix() → ndarray[source]

Exponential position matrix

I think this method is not necessary.

class pytdscf.HarmonicOscillator(ngrid: int, omega: float, q_eq: float | None = 0.0, units: str | None = 'cm-1', dimnsionless: bool | None = False, doAnalytical: bool | None = True)[source]

Bases: DVRPrimitivesMixin

Harmonic Oscillator DVR functions

See also MCTDH review Phys.Rep. 324, 1 (2000) appendix B https://doi.org/10.1016/S0370-1573(99)00047-2

Normalization factor :

\[A_n = \frac{1}{\sqrt{n! 2^n}} \left(\frac{m\omega}{\pi\hbar}\right)^{\frac{1}{4}}\ \xrightarrow[\rm impl]{} \frac{1}{\sqrt{n! 2^n}} \left(\frac{\omega_i}{\pi}\right)^{\frac{1}{4}}\]

Dimensionless coordinate :

\[\zeta = \sqrt{\frac{m\omega} {\hbar}}(x-a) \xrightarrow[\rm impl]{} \sqrt{\omega_i}q_i\]

Primitive Function :

\[\varphi_n = A_n H_n(\zeta) \exp\left(- \frac{\zeta^2}{2}\right)\quad (n=0,2,\ldots,N-1)\]
ngrid

# of grid

Type:

int

nprim

# of primitive function. same as ngrid

Type:

int

omega

frequency in a.u.

Type:

float

q_eq

eq_point in mass-weighted coordinate

Type:

float

dimensionless

Input is dimensionless coordinate or not.

Type:

bool

doAnalytical

Use analytical integral or diagonalization.

Type:

bool

diagnalize_pos_rep_matrix()[source]

Analytical formulation has not yet derived.

fbr_func(n: int, q)[source]

Primitive functions :math:\varphi_n(q)

Parameters:

n (int) – index (0 <= n < ngrid)

Returns:

\(\varphi_n(q)\)

Return type:

float

get_1st_derivative_matrix_dvr()[source]

\(\langle\chi_\alpha|\frac{d}{dq}|\chi_\beta\rangle\)

get_1st_derivative_matrix_fbr()[source]

Analytical Formulation

\[D_{j k}^{(1)}=\ \left\langle\varphi_{j} | \frac{\mathrm{d}}{\mathrm{d} x} | \varphi_{k}\right\rangle=\ -\sqrt{\frac{m \omega}{2}}\left(\sqrt{j+1} \delta_{j, k-1}-\sqrt{j} \delta_{j, k+1}\right)\]
get_2nd_derivative_matrix_dvr()[source]

\(\langle\chi_\alpha|\frac{d^2}{dq^2}|\chi_\beta\rangle\)

get_2nd_derivative_matrix_fbr()[source]

Analytical Formulation

\[D_{j k}^{(2)}=\ \frac{m\omega}{2}\left(\sqrt{(j-1)j}\delta_{j,k+2}-(2j+1)\delta_{j,k} + \sqrt{(j+2)(j+1)}\delta_{j,k-2}\right)\]
get_ovi_CS_HO(p: float, q: float, type: str = 'DVR') → ndarray[source]

Get overlap integral 1D-array between coherent state <p,q,ω,| and |HO(ω)>

Parameters:

  • p (float) – momentum of coherent state in mass-weighted a.u.

  • q (float) – position of coherent state in mass-weighted a.u.

  • type (str) – Whether “DVR” or “FBR”. Default is “DVR”.

Returns:

overlap integral 1D-array between coherent state <p,q,ω,| and |HO(ω)>

Return type:

np.ndarray

get_pos_rep_matrix()[source]

HO has analytical formulation

\[\left\langle\varphi_{j}|\hat{x}| \varphi_{k}\right\rangle=\ \sqrt{\frac{j+1}{2 m \omega}} \delta_{j, k-1}+x_{\mathrm{eq}}\ \delta_{j k}+\sqrt{\frac{j}{2 m \omega}} \delta_{j, k+1}\]
class pytdscf.Model(basinfo: BasInfo, operators: dict[str, HamiltonianMixin], *, build_td_hamiltonian: PolynomialHamiltonian | None = None, space: Literal['hilbert', 'liouville'] = 'hilbert')[source]

Bases: object

The wavefunction and operator information class

Parameters:

  • basinfo (BasInfo) – The wavefunction basis information

  • operators (Dict) – The operators. operators[name] gives a operator named name, such as hamiltonian.

  • build_td_hamiltonian (PolynomialHamiltonian) – Time dependent PolynomialHamiltonian. Defaults to None.

init_weight_VIBSTATE

Initial weight of VIBSTATE. List length is [nstate, ndof].

Type:

List[List[float]]

init_weight_GRID

Initial weight of GRID. List length is [nstate, ndof].

Type:

List[List[float]]

init_weight_ESTATE

Initial weight of ESTATE. List length is nstate.

Type:

List[float]

ints_prim_file

The file name of the primitive integrals

Type:

str

m_aux_max

The maximum number of auxiliary basis.

Type:

int

basinfo

The wavefunction basis information

Type:

BasInfo

hamiltonian

PolynomialHamiltonian

Type:

PolynomialHamiltonian

observables

Observable operators, such as PolynomialHamiltonian, dipole moment, occupation number etc.

Type:

Dict

build_td_hamiltonian

Time-dependent PolynomialHamiltonian.

Type:

PolynomialHamiltonian

get_ndof() → int[source]
Returns:

Degree of freedoms

Return type:

int

get_ndof_per_sites()[source]

N.Y.I ?

get_nprim(istate: int, idof: int) → int[source]
Parameters:

  • istate (int) – index of electronic states

  • idof (int) – index of degree of freedom

Returns:

The number of primitive basis in istate, idof.

Return type:

int

get_nspf(istate: int, idof: int) → int[source]
Parameters:

  • istate (int) – index of electronic states

  • idof (int) – index of degree of freedom

Returns:

The number of SPF in istate, idof.

Return type:

int

get_nspf_list(istate: int) → list[int][source]
Parameters:

istate (int) – index of electronic states

Returns:

Number of SPFs. e.g. [2, 2]

Return type:

List(int)

get_nspf_rangelist(istate: int) → list[list[int]][source]
Parameters:

istate (int) – index of electronic states

Returns:

the indices of SPFs. e.g. [[0,1,2],[0,1,2]]

Return type:

List[List[int]]

get_nstate() → int[source]
Returns:

Number of electronic states

Return type:

int

get_primbas(istate: int, idof: int) → PrimBas_HO[source]
Parameters:

  • istate (int) – index of electronic states

  • idof (int) – index of degree of freedom

Returns:

The primitive basis in istate, idof.

Return type:

primints_cls.PrimBas_HO

init_HartreeProduct: list[list[list[float]]] | None = None
init_weight_ESTATE: list[float] | None = None
init_weight_VIBSTATE: list[list[float]] | None = None
ints_prim_file: str | None = None
m_aux_max: int | None = None
space: Literal['hilbert', 'liouville']
class pytdscf.PrimBas_HO(origin: float, freq_cm1: float, nprim: int, origin_is_dimless: bool = True)[source]

Bases: object

The Harmonic Oscillator eigenfunction primitive basis.

This class holds information on Gauss Hermite type SPF basis functions. The n-th primitive represents

Normalization factor :

\[A_n = \frac{1}{\sqrt{n! 2^n}} \left(\frac{m\omega}{\pi\hbar}\right)^{\frac{1}{4}}\]

Dimensionless coordinate :

\[\zeta = \sqrt{\frac{m\omega} {\hbar}}(x-a)\]

Primitive Function :

\[\chi_n = A_n H_n(\zeta) \exp\left(- \frac{\zeta^2}{2}\right)\]
Parameters:

  • origin (float) – The center (equilibrium) dimensionless coordinate \(\sqrt{\frac{m\omega}{\hbar}} a\) of Hermite polynomial.

  • freq_cm1 (float) – The frequency \(\omega\) of Hermite polynomial. The unit is cm-1.

  • nprim (int) – The number (max order) of primitive basis on a certain SPF.

  • origin_is_dimless (bool, optional) – If True, given self.origin is dimensionless coordinate. Defaults to True.

origin

The center (equilibrium) dimensionless coordinate \(\sqrt{\frac{m\omega}{\hbar}} a\) of Hermite polynomial.

Type:

float

freq_cm1

The frequency \(\omega\) of Hermite polynomial. The unit is cm-1.

Type:

float

nprim

The number (max order) of primitive basis on a certain SPF.

Type:

int

freq_au

self.freq_cm1 in a.u.

Type:

float

origin_mwc

Mass weighted coordinate self.origin in a.u.

Type:

float

origin

origin is mass-weghted

todvr()[source]
class pytdscf.Simulator(jobname: str, model: Model, ci_type: Literal['mps', 'mctdh', 'ci'] = 'mps', backend: Literal['jax', 'numpy'] = 'jax', proj_gs: bool = False, t2_trick: bool = True, verbose: int = 2)[source]

Bases: object

The simulator of the PyTDSCF

set parameter of the restart, propagate, operate dipole, save_file etc …

Parameters:

  • jobname (str) – the jobname

  • model (model_cls.Model) – run parameter (basis, hamiltonian, observable, bond dimension, initial_weight etc.)

  • t2_trick (bool, optional) – Use so-called t/2 trick in auto-correlation. Note that it requires initial state to be real. Defaults to True.

  • ci_type (str, optional)

  • backend (str, optional) – JAX or Numpy. Defaults to 'jax'. When polynomial operator, FBR basis and small bond-dimension is used, 'Numpy' is recommended.

  • proj_gs (bool, optional) – Initial state is projected from the ground state. Defaults to False. If proj_gs=True, one must be set attribute model.primbas_gs: List[Primbas_HO].

backup_interval: int
get_initial_wavefunction(ints_prim: PrimInts) → WFunc[source]
get_primitive_integrals() → PrimInts[source]
maxstep: int
operate(maxstep: int = 10, restart: bool = False, savefile_ext: str = '_operate', loadfile_ext: str = '_gs', verbose: int = 2) → tuple[float, WFunc][source]

Apply operator such as dipole operator to the wavefunction

Parameters:

  • maxstep (int, optional) – Maximum number of iteration. Defaults to 10.

  • restart (bool, optional) – Restart from the previous wavefunction. Defaults to False.

  • backend (str, optional) – JAX or Numpy. Defaults to 'jax'.

  • savefile_ext (str, optional) – Extension of the save file. Defaults to '_operate'.

  • loadfile_ext (str, optional) – Extension of the load file. Defaults to '_gs'. When restart=False, loadfile_ext is ignored.

  • verbose (int, optional) – Verbosity level. Defaults to 2.

Returns:

norm of O|Ψ> and Wavefunction after operation.

Return type:

Tuple[float, WFunc]

propagate(stepsize: float = 0.1, maxstep: int = 5000, restart: bool = False, savefile_ext: str = '', loadfile_ext: str = '_operate', backup_interval: int = 1000, autocorr: bool = True, energy: bool = True, norm: bool = True, populations: bool = True, observables: bool = False, reduced_density: tuple[list[tuple[int, ...]], int] | None = None, Δt: float | None = None, thresh_sil: float = 1e-09, autocorr_per_step: int = 1, observables_per_step: int = 1, energy_per_step: int = 1, norm_per_step: int = 1, populations_per_step: int = 1, parallel_split_indices: list[tuple[int, int]] | None = None, adaptive: bool = False, adaptive_Dmax: int = 20, adaptive_dD: int = 5, adaptive_p_proj: float = 0.0001, adaptive_p_svd: float = 1e-07, integrator: Literal['lanczos', 'arnoldi'] = 'lanczos', step_size_is_fs: bool = True) → tuple[float, WFunc][source]

Propagation

Parameters:

  • stepsize (float, optional) – Step size in fs. Defaults to 0.1.

  • maxstep (int, optional) – Maximum number of steps. Defaults to 5000., i.e. 500 fs.

  • restart (bool, optional) – Restart from the previous wavefunction. Defaults to False.

  • savefile_ext (str, optional) – Extension of the save file. Defaults to ''.

  • loadfile_ext (str, optional) – Extension of the load file. Defaults to '_operate'. When restart=False, loadfile_ext is ignored.

  • backup_interval (int, optional) – Number of steps at which, the wavefunction is saved. Defaults to 1000.

  • autocorr (bool, optional) – Calculate autocorrelation function. Defaults to True.

  • energy (bool, optional) – Calculate energy. Defaults to True.

  • norm (bool, optional) – Calculate norm. Defaults to True.

  • populations (bool, optional) – Calculate populations. Defaults to True.

  • observables (bool, optional) – Calculate observables. Defaults to False.

  • reduced_density (Dict[Tuple[int, ...], int], optional) – Calculate reduced density of the given modes. For example, ([(0, 1),], 10) means calculate the diagonal elements of reduced density of the \(\rho_{01}=\mathrm{Tr}_{p\notin \{0,1\}}\left|\Psi^{(\alpha)}\rangle\langle\Psi^(\alpha)\right|\) in per 10 steps. Note that it requires enough disk space. Defaults to None. It is better if the target modes are close to rightmost, i.e., 0. (Because this program calculate property in the most right-canonized form of MPS.) If you want coherence, i.e., off-diagonal elements of density matrix, you need to set like ([(0, 0), ], 10).

  • Δt (float, optional) – Same as stepsize

  • thresh_sil (float) – Convergence threshold of short iterative Lanczos. Defaults to 1.e-09.

  • step_size_is_fs (bool, optional) – If True, stepsize is in fs. Defaults to True.

Returns:

Energy during propagation (it conserves) and Wavefunction after propagation.

Return type:

Tuple[float, WFunc]

relax(stepsize: float = 0.1, maxstep: int = 20, improved: bool = True, restart: bool = False, savefile_ext: str = '_gs', loadfile_ext: str = '', backup_interval: int = 10, norm: bool = True, populations: bool = True, observables: bool = False, integrator: Literal['lanczos', 'arnoldi'] = 'lanczos') → tuple[float, WFunc][source]

Relaxation

Parameters:

  • stepsize (float, optional) – Step size in fs. Defaults to 0.1. This is used only when imaginary time propagation is used.

  • maxstep (int, optional) – Maximum number of steps. Defaults to 20.

  • improved (bool, optional) – Use improved relaxation. Defaults to True.

  • restart (bool, optional) – Restart from the previous wavefunction. Defaults to False.

  • savefile_ext (str, optional) – Extension of the save file. Defaults to '_gs'.

  • loadfile_ext (str, optional) – Extension of the load file. Defaults to ''. When restart=False, loadfile_ext is ignored.

  • backup_interval (int, optional) – Number of steps at which, the wavefunction is saved. Defaults to 10.

  • norm (bool, optional) – Calculate norm. Defaults to True.

  • populations (bool, optional) – Calculate populations. Defaults to True.

  • observables (bool, optional) – Calculate observables. Defaults to False.

Returns:

Energy after relaxation in Eh, and Wavefunction after relaxation.

Return type:

Tuple[float, WFunc]

save_wavefunction(wf: WFunc, log: bool = False)[source]
stepsize: float
class pytdscf.Sine(ngrid: int, length: float, x0: float = 0.0, units: str = 'angstrom', doAnalytical: bool = True, include_terminal: bool = True)[source]

Bases: DVRPrimitivesMixin

Sine DVR functions Note that Sine DVR position matrix is not tridiagonal! Index starts from j=1, alpha=1

Terminal points (x_{0} and x_{N+1}) do not belog to the grid.

See also MCTDH review Phys.Rep. 324, 1 (2000) appendix B https://doi.org/10.1016/S0370-1573(99)00047-2

Primitive Function :

\[\begin{split}\varphi_{j}(x)= \begin{cases}\sqrt{2 / L} \sin \left(j \pi\left(x-x_{0}\right) / L\right) \ & \text { for } x_{0} \leq x \leq x_{N+1} \\ 0 & \text { else }\end{cases} \quad (j=1,2,\ldots,N)\end{split}\]
ngrid

Number of grid. Note that which does not contain terminal point.

Type:

int

nprim

Number of primitive function. sama as ngrid.

Type:

int

length

Length in a.u.

Type:

float

x0

start point in a.u.

Type:

float

doAnalytical

Use analytical integral or diagonalization.

Type:

bool

include_terminal

Whether include terminal grid.

Type:

bool

diagnalize_pos_rep_matrix()[source]

Analytical diagonalization has been derived

\[U_{j \alpha}=\sqrt{\frac{2}{N+1}} \sin \left(\frac{j \alpha \pi}{N+1}\right)\]
\[z_{\alpha}=\cos \left(\frac{\alpha \pi}{N+1}\right)\]

This leads to the DVR grid points

\[x_{\alpha}=x_{0}+\frac{L}{\pi} \arccos \left(z_{\alpha}\right)\ =x_{0}+\alpha \frac{L}{N+1}=x_{0}+\alpha \Delta x\]
fbr_func(n: int, x: float)[source]

Primitive functions :math:\varphi_n(x)

Parameters:

n (int) – index (0 <= n < ngrid)

Returns:

\(\varphi_n(x)\)

Return type:

float

get_1st_derivative_matrix_dvr()[source]

\(\langle\chi_\alpha|\frac{d}{dq}|\chi_\beta\rangle\)

get_1st_derivative_matrix_fbr()[source]

Analytical Formulation

\[D_{j k}^{(1)}=\ \left\langle\varphi_{j} |\frac{\mathrm{d}}{\mathrm{d} x} | \varphi_{k}\right\rangle={\rm mod} (j-k, 2) \frac{4}{L} \frac{j k}{j^{2}-k^{2}},\ \quad j \neq k\]
get_2nd_derivative_matrix_dvr()[source]

Analytical forumulation exists.

\[\begin{split}D_{\alpha \beta}^{(2), \mathrm{DVR}}=-\left(\frac{\pi}{\Delta x}\right)^{2}\left\{\begin{array}{l} -\frac{1}{3}+\frac{1}{6(N+1)^{2}}-\frac{1}{2(N+1)^{2} \sin ^{2}\left(\frac{\alpha \pi}{N+1}\right)}, \quad \alpha=\beta \\ \frac{2(-1)^{\alpha-\beta}}{(N+1)^{2}} \frac{\sin \left(\frac{\alpha \pi}{N+1}\right) \sin \left(\frac{\beta \pi}{N+1}\right)}{\left(\cos \left(\frac{\alpha \pi}{N+1}\right)-\cos \left(\frac{\beta \pi}{N+1}\right)\right)^{2}},\quad \alpha \neq \beta \end{array} \right.\end{split}\]
get_2nd_derivative_matrix_fbr()[source]

Analytical Formulation (start from j=1 !)

\[D_{j k}^{(1)}=\ \left\langle\varphi_{j} |\frac{\mathrm{d}}{\mathrm{d} x} | \varphi_{k}\right\rangle={\rm mod} (j-k, 2) \frac{4}{L} \frac{j k}{j^{2}-k^{2}},\ \quad j \neq k\]
get_pos_rep_matrix()[source]

Sine position matrix

\[Q_{j k}=\left\langle\varphi_{j}|\hat{z}| \ \varphi_{k}\right\rangle=\frac{1}{2}\left(\delta_{j, k+1}+\delta_{j, k-1}\right)\]

where transformed variable

\[z=\cos \left(\pi\left(x-x_{0}\right) / L\right)\]

is introduced.

class pytdscf.TensorHamiltonian(ndof: int, potential: list[list[dict[tuple[int | tuple[int, int], ...], TensorOperator]]], name: str = 'hamiltonian', kinetic: list[list[dict[tuple[tuple[int, int], ...], TensorOperator] | None]] | None = None, decompose_type: Literal['QRD', 'SVD'] = 'QRD', rate: float | None = None, bond_dimension: list[int] | int | None = None, backend: Literal['jax', 'numpy'] = 'jax')[source]

Bases: HamiltonianMixin

Hamiltonian in tensor formulation.

mpo

MPO between i-state and j-state

Type:

List[List[MatrixProductOperators]]

distribute_mpo_cores()[source]
mpo: list[list[MatrixProductOperators | None]]
class pytdscf.TensorOperator(*, shape: tuple[int, ...] | None = None, tensor: ndarray | None = None, only_diag: bool = False, legs: tuple[int, ...] | None = None, name: str | None = None, mpo: list[ndarray] | None = None)[source]

Bases: object

Tensor Operator class

shape

Tensor shape

Type:

Tuple[int]

tensor_orig

Original tensor

Type:

np.ndarray

tensor_decomposed

Decomposed tensor

Type:

np.ndarray

tensor_full

if only_diag==False, tensor_orig is zero-filled.

Type:

np.ndarray

only_diag

Only diagonal terms are non-zero. Defaults to False.

Type:

bool

legs

Tensor legs

Type:

Tuple[int]

name

Tensor Operator name. Defaults to None.

Type:

str

Tensor Example:

\[\begin{split}\begin{array}{c} j_0^\prime & & j_1^\prime & & & & j_3^\prime \\ | & & | & & & & | \\ W(0) & - & W(1) & - & - & - & W(3) \\ | & & | & & & & | \\ j_0 & & j_1 & & & & j_3 \\ \end{array}\\\end{split}\]
\[j_0 = 0, 1, .. n-1, j_1 = 0, 1, .. m-1, j_3 = 0, 1, .. l-1\]

tensor = np.array([W0, W1, W3]).

When only_diag=False, above example is shape = (n, n, m, m, l, l),

legs=(0, 0, 1, 1, 3, 3), (even index is bra, odd index is ket).

Otherwise, only_diag=True, shape=(n,m,l), legs=(0, 1, 3).

W_1[beta_01,j_1^prime,j_1,beta_13] = \(W{\small \beta_{0,1}}_{j_1}^{j_1^\prime}{\small \beta_{1,3}}\)

add_dammy_legs(add_legs: tuple[int, ...], add_shape: tuple[int, ...]) → TensorOperator[source]

Dammy legs addition

Parameters:

  • add_legs (Tuple[int]) – additional tensor legs. e.g. add \(j_1 ` is ``add_legs=(1,)`\).

  • add_shape (Tuple[int]) – additional tensor shape. e.g. add \(j_1 = 0,1,\ldots,3\) is add_shape=(4,).

Returns:

Filled Tensor

Return type:

TensorOperator

e.g. from legs = (0,1) to legs=(0,1,2), where dammy legs is 2.

bond_dimension: list[int]
decompose(bond_dimension: list[int] | int | None = None, decompose_type: str = 'SVD', rate: float | None = None, square_sum: bool = True, overwrite: bool = False) → list[ndarray][source]

MPO Decomposition

Parameters:

  • bond_dimension (List[int] or int) – MPO bond dimension

  • decompose_type (str) – Tensor Train Decomposition Type 'QRD' or 'SVD'. Defaults to 'SVD'.

  • rate (float) – SVD decomposition contribution rate truncation. Defaults to None.

  • square_sum (bool) – Whether rate is based on the sum of squares of the eigenvalues, or simply sum of the eigenvalues. Square sum can calculate easily by trace of matrix, but need high rate percentage. We recommend rate=0.99999999999, if square_sum=True.

Returns:

MPO

Return type:

List[np.ndarray]

property dtype: dtype
estimate_error(tensor: ndarray | None = None)[source]

Estimate Frobenius norm between self and given-tensor or renormalized tensor

Parameters:

tensor (np.ndarray) – given tensor (option)

Returns:

Error %

Return type:

float

full_bond_dimension: list[int]
get_bond_dimension(bond_dimension: list[int] | int | None = None)[source]

get MPO bond-dimension

Parameters:

bond_dimension (List[int] | int) – MPO bond dimension

get_frobeinus_norm() → float[source]

Get Frobenius norm of tensor operator

get_num_of_elements() → tuple[int, int | None][source]

Get number of decomposed tensor operator element

Returns:

Number of tensor elements (before, after) decomposition

Return type:

Tuple[int, int]

get_tensor_full() → ndarray[source]

fill off-diagonal term zero if only_diag=True

legs: tuple[int, ...]
classmethod load(name) → TensorOperator[source]

Load Tensor Operator Object from binary file in ./tensor_operators directory

name: str
only_diag: bool
restore_from_decoposed() → ndarray[source]

Restore tensor_orig from MPO (tensor_decomposed). This routine is almost the same as mps_cls.buid_CICoef

save(name: str | None = None) → None[source]

Save Tensor Operator Object as binary file in ./tensor_operators directory

Parameters:

name (str) – file name not including extension

shape: tuple[int, ...]
tensor_decomposed: list[ndarray]
tensor_full: ndarray
tensor_orig: ndarray
pytdscf.construct_fulldimensional(dvr_prims: list[DVRPrimitivesMixin], func: Callable | None = None, db: str | None = None, ref_ene: float = 0.0, dipole: bool = False, efield: tuple[float, float, float] = (0.01, 0.01, 0.01)) → dict[tuple[int, ...], TensorOperator][source]

Construct full-dimensional Operator from DVR grid-based PES

Parameters:

  • dvr_prims (List[DVRPrimitivesMixin]) – DVR functions. Sorted in Database order.

  • func (Optional,Callable) – Potential energy function

  • db (Optional,str) – Electronic structure calculation database path, such as ‘foo/hoge.db’

  • ref_ene (Optional, float) – Opt coord energy in a.u., if you need subtract from function or database energy. Defaults to 0.0

  • dipole (Optional, bool) – Get dipole moment from database. Defaults to False

  • efield (Optional, List[float, float, float]) – Electronic field for inner products with dipole moment. Defaults to [1.0,1.0,1.0]

Returns:

full-dimensional tensor operator

Return type:

Dict[Tuple[int], TensorOperator]

pytdscf.construct_kinetic_operator(dvr_prims: list[DVRPrimitivesMixin], coefs: list[float] | None = None, forms: str = 'mpo') → dict[tuple[tuple[int, int], ...], TensorOperator][source]

Note

The off-diagonal terms of the kinetic operator are not zero in DVR basis.

Kinetic energy operator (KEO) in linear coordinate is represented by following matrix product operator (MPO):

\[\begin{split}\begin{pmatrix} -\frac{\hat{P}_1^2}{2} & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ -\frac{\hat{P}_2^2}{2} & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ -\frac{\hat{P}_3^2}{2} & 1 \end{pmatrix} \begin{pmatrix} 1 \\ -\frac{\hat{P}_4^2}{2} \end{pmatrix}\end{split}\]
Parameters:

  • dvr_prims (List[DVRPrimitivesMixin]) – DVR functions

  • coefs (List[float]) – Coefficients of kinetic operator. Defaults to [1.0, 1.0, ...], i.e., \(\sum_i -\frac{1}{2}\frac{d^2}{dQ_i^2}\) is given. For example, when one uses dimensionless coordinate and rotational coordinate, the kinetic operator is given by \(\hat{T} = -\frac{\omega}{2} \frac{d^2}{dx^2} - \frac{1}{2I} \frac{d^2}{d\theta^2}\),then coefs should be given by \([\omega, 1/I]\).

  • forms (str) – Either ‘sop’ or ‘mpo’. Defaults to ‘mpo’.

pytdscf.construct_nMR_recursive(dvr_prims: list[DVRPrimitivesMixin], nMR: int = 3, ndof: int | None = None, func: dict[tuple[int, ...], Callable] | None = None, db: str | None = None, df: DataFrame | None = None, active_dofs: list[int] | None = None, site_order: dict[int, int] | None = None, zero_indices: list[int] | None = None, return_tensor: bool = False, include_const_in_mpo: bool = False, ref_ene: float | None = None, dipole: bool = False, efield: tuple[float, float, float] = (1.0, 1.0, 1.0), rate: float = 1.0, k: int = 200, nsweep: int = 4) → list[ndarray] | dict[tuple[int, ...], TensorOperator][source]

Construct n-Mode Representation Operator

n-Mode Representation (nMR) are used for reducing grid points for ab initio calculation. To avoid duplicate addition of the same coordinates, use principle of inclusion and exclusion. (It is easy to understand by Venn-diagram)

Parameters:

  • dvr_prims (List[DVRPrimitivesMixin]) – DVR functions

  • nMR (Optional, int) – Mode Representation Number. Defaults to 3.

  • ndof (Optional, int) – number of mode including in database.

  • func (Optional, Dict[Tuple[int], Callable]) – E.g. Potential energy function

  • db (Optional, str) – Electronic structure calculation database path, such as ‘foo/hoge.db’

  • df (Optional, pl.DataFrame) – Electronic structure calculation dataframe.

  • active_dofs (Optional, List[int]) – Active modes Defaults to ALL.

  • site_order (Optional, Dict[int, int]) – MPO site order. Defaults to DOF index order.

  • zero_indices (Optional, List[int]) – nMR criteria coordinate grid.

  • return_tensor (bool) – Return before decomposition. Defaults to False.

  • include_const_in_mpo (bool) – Include scalar constant value (such as refernce energy) in MPO. Defaults to False.

  • ref_ene (float) – Opt coord energy in a.u., if you need subtract from function or database energy. Defaults to 0.0

  • dipole (bool) – Get dipole moment from database. Defaults to False

  • efield (List[float]) –

    Electronic field for inner products with dipole moment.Defaults to [1.0e-02 , 1.0e-02, 1.0e-02]

  • rate (float) – SVD contribution rate. Defaults to 1.0

  • k (int) – SVD bond dimension. Defaults to 200

Returns:

nMR MPO

Return type:

List[np.ndarray]

Example:

\[\begin{split}V &:= \sum_{\boldsymbol{x}\in \mathrm{coord.}} V(x_\alpha,x_\beta) \\ &= V_0 + V_{\mathrm{1MR}} + V_{\mathrm{2MR}}\\\end{split}\]

where,

\[\begin{split}V_0 &= V(0,0)\\ V_{\mathrm{1MR}} &= \sum_{x_\alpha} \left(V(x_\alpha,0) - V_0\right) + \sum_{x_\beta} \left(V(0, x_\beta) - V_0\right)\\ V_{\mathrm{2MR}} &= \sum_{x_\alpha}\sum_{x_\beta} \ \left(V(x_\alpha,x_\beta) - V(x_\alpha,0)- V(0, x_\beta) - V_0\right)\end{split}\]
pytdscf.read_potential_nMR(potential_emu: dict[tuple[int, ...], float | complex], *, active_modes: list[int] | None = None, name: str = 'hamiltonian', cut_off: float | None = None, dipole_emu: dict[tuple[int, ...], tuple[float, float, float]] | None = None, print_out: bool = False, active_momentum=None, div_factorial: bool = True, efield: tuple[float, float, float] = (1.0, 1.0, 1.0)) → PolynomialHamiltonian[source]

Construct polynomial potential, such as n-Mode Representation form. multi-state nMR is not yet implemented.

Parameters:

  • potential_emu (Dict[Tuple[int], float or complex]) –

    The coeffcients of operators. The index is start from 1. Factorial factor is not included in the coefficients. E.g. V(q_1) = 0.123 * 1/(2!) * q_1^2 gives {(1,1):0.123}.

    units are all a.u.:

    k_orig{ij} in [Hartree/emu/bohr**2]

    k_orig{ijk} in [Hartree/emu**i(3/2)/BOHR**3]

    k_orig{ijkl} in [Hartree/emu**(2)/BOHR**4]

  • active_modes (List[int]) – The acctive degree of freedoms. Note that, The order of the MPS lattice follows this list. The index is start from 1.

  • name (str) – The name of operator.

  • cut_off (float) – The threshold of truncatig terms by coeffcients after factorial division.

  • dipole_emu (Dict[Tuple[int], float or complex]) –

    Dipole moment operator to generate a vibrational excited states. Defaults to None. The definition is almost the same as potential_emu, but this value is 3d vector (list).

    mu{ij} in [Debye/emu/bohr**2]

    mu{ijk} in [Debye/emu**(3/2)/BOHR**3]

    mu{ijkl} in [Debye/emu**(2)/BOHR**4]

  • print_out (bool) – Defaults to False

  • div_factorial (bool) – Whether or not divided by factorial term. Defaults to True.

  • active_momentum (List[Tuple[int, float]]) – [(DOF, coef)]. Defaults to None When active_momentum are set, unselected kinetic energy operators are excluded. Default option include all kinetic energy operators.

  • efield (List[float]) – The electric field in [Hartree Bohr^-1 e^-1]. Defaults to [1.0, 1.0, 1.0]

Returns:

operator

Return type:

PolynomialHamiltonian