4.2. Matrix classes¶
Methods implemented in PyBEST operate on and return custom multidimensional array classes that store the actual arrays as their private attributes (_array). The linalg module contains two types of matrix representations: dense matrices and Cholesky matrices. Algebraic operations are implemented as methods.
4.2.1. DenseNIndex (N = One, Two, Three, Four) classes¶
In the dense representation, all the elements are stored in memory. The dense matrix can be created using DenseLinalgFactory or directly by using the constructor
my_matrix = DenseThreeIndex(n1, n2, n3)
where arguments n1, n2, and n3 are three integers defining the dimensionality of the three-index object. DenseNIndex objects are usually returned as products of mathematical operations.
4.2.2. CholeskyFourIndex class¶
The Cholesky Decomposition representation operates on the four-index object A that is internally stored in the form of two three-index object which fulfills A = LL*. The interface is analogous to the DenseFourIndex class interface, so CholeskyFourIndex instances are treated as four-index arrays.
4.3. LinalgFactory classes¶
The DenseNIndex and CholeskyFourIndex matrices can be created using linear algebra factory instances - DenseLinalgFactory and CholeskyLinagsFactory. See examples to compare different methods of creating new N-index multidimensional arrays.
4.4. Using PyBEST classes to perform algebraic operations¶
4.4.1. Creation of dense multidimensional arrays¶
from pybest.linalg.dense import DenseLinalgFactory, DenseTwoIndex
# Create 9x9 and 9x9x9 dense matrices using factory instance.
lf = DenseLinalgFactory(9)
matrix2d_0 = lf.create_two_index()
matrix3d_0 = lf.create_three_index()
# Create 3x9 dense matrix using direct initialization.
matrix2d_1 = DenseTwoIndex(3, 9)
# Create an empty matrix with the same shape as matrix2d_0.
matrix2d_2 = matrix2d_0.new()
# Create a copy of the matrix.
matrix2d_3 = matrix2d_0.copy()
4.4.2. Filling matrix with numbers¶
# Fill matrix with random numbers.
my_matrix.randomize()
# Fill diagonal elements of the subblock [3:6, 3:6] (defined by begin and
# end keywords) of the matrix.
my_matrix.assign_diagonal(2, begin=3, end=6).
# Set all elements of a matrix to 0
my_matrix.assign(0)
# Set all elements to 3 in a subblock [2:4,1:3] of my_matrix.
my_matrix.assign(3, begin0=2, end0=4, begin1=1, end1=3)
print(my_matrix.get_element[2,2]) # Prints 3
print(my_matrix.get_element[0,0]) # Prints 0
# Fill matrix with the content of another matrix.
my_matrix.assign(other_matrix)
# Fill a subblock (defined using begin0, end0, begin1, end1 keywords) of
# my_matrix with content of a subblock of another matrix (defined by
# begin2, end2, begin3, end3 keywords).
my_matrix.assign(other_matrix, begin0=2, begin2=2)
4.4.3. General NIndexObject multiplication - contraction¶
Evaluates the Einstein summation convention on the operands that are CholeskyFourIndex or DenseNIndex (where N = One, Two, Three, Four) instances.
4.4.3.1. Arguments¶
- subscript
(string) specifies the subscripts for summation as comma separated list of subscript labels, for example: ‘abcd,efcd->abfe’.
- operands
(DenseNIndex) these are the other arrays for the operation. If out keyword is not specified and the result of operation is not a scalar value, the last operand is treated as out.
4.4.3.2. Keyword arguments¶
- out
(DenseNIndex) The product of operation is added to out.
- factor
(float) The product of operation is multiplied by a factor before it is added to out.
- clear
(boolean) The out is cleared before the product of operation is added to out if clear is set to True.
- select
(string) Specifies the contraction algorithm: One of:
‘opt_einsum’ - opt_einsum.contract,
‘einsum’ - numpy.einsum function,
‘td’ - operations with utilization of numpy.tensordot function,
‘blas’ - operations written in c++,
None - first tries ‘td’ or ‘blas’ routine, then “einsum”
‘td’ is usually much faster but it requires more memory and it is available only for specific cases.
- optimize
(False, True, ‘optimal’) Specifies contraction path for operation if select is set to ‘einsum’. For details, see numpy.einsum_path.
- begin0, end0, begin1, …
(int) The lower and upper bounds for summation.
# Matrix-vector multiplication.
lf = DenseLinalgFactory(4)
my_vector = lf.create_one_index()
my_matrix = lf.create_two_index()
my_result = my_matrix.contract("ab,b", my_vector)
# Four-index object contracted with the subblock of the two-index object.
# output_2 is a transposed output_1.
output_1 = array_4d.contract("abcd,cd->ab", array_2d)
output_2 = array_4d.contract("abcd,cd->ba", array_2d)
# Multiply a result of matrix-matrix multiplication by factor -1 and add it
# to another matrix.
matrix_0.contract("ab,ba", matrix_1, out=matrix_2, factor=-1)
# There is no difference in syntax if you contract DenseFourIndex or
# CholeskyFourIndex object.
output_1 = dense4ind_object.contract("abcd,cd->ab", array_2d)
output_2 = cholesky4ind_object.contract("abcd,cd->ab", array_2d)
Keyword arguments begin0, end0, begin1, … can be used to perform contraction on a subblock of an N-dimentional array. In the following example, we perform the operation on the subblock of the electron repulsion integrals that are represented as CholeskyFourIndex object. Please note, that there is no need to use all of the begin0, end0, begin1, … keyword arguments since their default values allows to include all indices.
# nocc = number of occupied orbitals
nocc = 5
# subblock_ooov = defines boundaries of the subblock of the array
subblock_ooov = {"end0": nocc, "end1": nocc, "end2": nocc, "begin3": nocc}
# eri = electron repulsion integrals (CholeskyFourIndex instance)
exchange_ooov = eri.contract("abcd->abcd", **subblock_ooov)
exchange_ooov.iadd_transpose((0, 2, 1, 3), factor=-2)
4.4.4. General NIndexObject slicing¶
The syntax for the slicing operation is analogous as for the contract method.
4.4.4.1. Arguments¶
- subscripts
(string) specifies the subscripts for summation as comma-separated list of subscript labels, for example: ‘abcd,efcd->abfe’.
4.4.4.2. Keyword arguments¶
- out
(DenseNIndex) The product of operation is added to out.
- factor
(float) The product of operation is multiplied by a factor before it is added to out.
- clear
(boolean) The out is cleared before the product of operation is added to out if clear is set to True.
- begin0, end0, begin1, …
(int) The lower and upper bounds for summation.
# Get a slice
array_3d = array_4d.contract("abcb->abc")
# Get a slice from a given subblock
subblock = {"begin0": 1, "end0": 9, "begin3": 1, "end3": 9}
array_2d = array_4d.contract("abba->ab", **subblock)
# Get a slice and add it to array_out after setting its values to zero
array_4d.contract("abcb->abc", out=array_out, clear=True)
4.4.5. Adding arrays¶
The kinetic energy operator and nuclear attraction energy operator are represented by DenseTwoIndex objects. A nonrelativistic one-body energy operator can be obtained by summing these two operators.
kinetic_energy_operator = compute_kinetic(obasis)
nuclear_attraction_energy_operator = compute_nuclear(obasis)
one_body_energy_operator = kinetic_energy_operator.copy()
one_body_energy_operator.iadd(nuclear_attraction_energy_operator)
4.4.6. Scaling by a scalar factor¶
any_nindex_object.iscale(0.5)