{'fraction_non_alphanumeric': 0.035161808583887996, 'fraction_numerical': 0.004360798714711958, 'mean_word_length': 4.388572841949048, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the possibility to use non-orthogonal orbitals for Density Matrix Renormalization Group calculations in Quantum Chemistry

The generalization of Density Matrix Renormalization Group (DMRG) approach as implemented in quantum chemistry, to the case of non-orthogonal orbitals is carefully analyzed. This generalization is attractive from the physical point of view since it allows a better localization of the orbitals. The possible implementation difficulties and drawbacks are estimated. General formulae for hamiltonian matrix elements useful in DMRG calculations are given.


I. INTRODUCTION
The DMRG method, first introduced by White [1,2,3] in solid state physics, has been very successful in calculations of low dimensional quantum lattice systems, and has become recently an interesting tool in quantum chemistry [4,5,6,7,8,9,10]. The method is attractive because of the good scaling properties of the calculation time with the size of the system [6]. In the present work we first list the main steps of DMRG algorithm as implemented in quantum chemistry. All previous implementations were restricted to the orthonormalized molecular orbitals (MOs). While rendering all algorithms much simpler, the standard MOs may be not quite appropriate from the physical point of view, especially if we try to follow the original DMRG ideas that assume the separation of the system into well defined physical blocks. From this point of view the use of bond structures, e.g. valence bond orbitals [11] can be much more attractive as it allows for clear interpretation of resulting wave functions. Even if it is not clear that this can improve the results numerically, we consider below the theoretical aspects of the use of non-orthogonal orbitals as applied to the DMRG algorithms. The article is organized as follows: first we describe essential steps of the classical, orthogonal, DMRG algorithm. In section III we describe the generalization to the case of non-orthogonal orbitals. In appendices we give the full list of relevant formulae.

II. CLASSICAL, ORTHOGONAL, DMRG ALGORITHM
We do not describe here the details of the DMRG approach, it can be found in many reviews [3,4,9] . We also do not describe the details of implementation of DMRG in quantum chemistry, see [6,8]. Rather we concern only with practical, computational, aspect of the algorithm. The main points of DMRG procedure are: • The parametrization of the wave function In DMRG we divide the whole system into A and U subsystems (two subsets of the full set of MOs in quantum chemistry). The wave function is represented as where |A I > describes some state of subsystem A and similar for U . The tensor product A I ⊗ U J is automatically antysymmetrized due to the presence of the creation operators (see [8] for details). The conservation rules imply for nonzero C IJ where n α is fixed by the total number of electrons and total spin projection, the same being true for β. If there present additional spatial symmetries, the corresponding conservation rules must be also satisfied. One important drawback of the above expression is that the wave function given by Eq. (1) is not necessarily the eigenfunction of total spin operatorŜ 2 . We actually never considered the application of corresponding projection operator to get the spin-pure DMRG wavefunction, maybe it is not really as difficult as theŜ 2 is the combination of creation/annihilation operators whose matrix elements are in general available in DMRG (see below).

B. State pair types and matrix elements to be stored
To find the wave function given by Eq.(1) we must diagonalize Hamiltonian matrix in the space span by |A I ⊗ U J > state vectors. To build this matrix all we need are some matrix elements for pure A and U blocks that we call 'primitive' matrix elements. The set of these matrix elements depends on the relations between number of α and β electrons. In the following when we consider particular matrix element indices I and K refer to A subsystem while J and L -to U subsystem. Let's denote with the symbol a † a creation operator referring to block A and with u † a creation operator referring to block U. The electronic part of quantum chemistry Hamiltonian can be written aŝ whereĤ A ,Ĥ U contain only a, a † or u, u † operators and the interaction term is given bŷ where p, q indices refer to orbitals of A block, x, y to U . Then, given I and K, the Table I  With this definition, the 'primitive' matrix elements needed to build the full Hamiltonian matrix, depend on the 'case' given above. The complete list of these matrix elements is given in Appendix A, together with resulting matrix elements of the total Hamiltonian (see below). Here we just note that all ME's have the typical form of with one to three creation/annihilation operators, or the linear combination of such ME's.

C. Calculation of Hamiltonian matrix elements
After we have defined the stored 'primitive' matrix elements in the previous subsection, the matrix elements of Hamiltonian matrix look like where type and range of index i depend on the IK-type. The full set of matrix elements is given in Appendix A, where the bi-orthogonal orbitals are assumed, see next sections. For the standard orthogonal case the dual orbitals coincide with original orbitals so the tildes in the formulas of Appendix A can be simply omitted.

D. Adding orbitals
Adding of orbital(s) is the most important and basic step of the whole DMRG procedure. Given the set of A states (the same is valid for U , so we do not speak of it further), we add few orbitals to A and build the new bigger set of states including all possible occupations of added orbital pattern. Indices Ii, Jj, etc. refer to the new enlarged subsystems. The subsequent step is the recalculation of the 'primitive' matrix elements in this new set of states. The only practical difficulty here is that the number of cases to treat grows considerably. In our early implementations we recalculated the 'primitive' matrix elements after adding orbitals and stored them on disk. This led to significant reductions in performance. Therefore later, to save disk space and improve performance, we changed our algorithms as not to calculate and to store the extended set of 'primitive' matrix elements but rather to build the Hamiltonian matrix in the extended space explicitly using the non-extended 'primitive' matrix elements and some additional molecular integrals. Also, the highly optimized BLAS [12,13,14,15,16] routines are used as much as possible to make the code faster. Note that this part of code can be easily parallelized. The full list of relevant formulas is given in Appendix B. The tildes can be safely omitted again there.

E. Hamiltonian diagonalization in the orthogonal case
Once the Hamiltonian matrix elements in |A I ⊗U J > basis are available, the next task is to diagonalize this Hamiltonian matrix. We can not apply the explicit techniques like Jacobi diagonalization, because of quite big, tens or hundreds of thousands states, dimension of typical problem. Therefore we are forced to use the direct technique, which is based on the ability to calculate the so called σ-vector, i.e. action of Hamiltonian given vector x without explicit storage of Hamiltonian matrix H αβ . In the usual formulation with orthogonal orbitals, Davidson method is the most efficient. As well known, it consists of an iterative procedure that starts from some reasonable guess and then expands the set of trial vectors up to say 30-50 that provide good approximation to exact eigenvector. The procedure can be seen as a minimization of the Rayleigh quotient in which one applies the Newton method and replace D = diag H in the Hessian matrix.

F. New states construction
After diagonalizing Hamiltonian matrix in the extended set of vectors, the next and the last task for the whole DMRG procedure is the construction of reduced subset of states describing A and U in the extended orbital set. This is done by finding the best possible approximation to the calculated eigenstate of Hamiltonian matrix when the number of states retained is fixed to a number (M ) smaller than the total dimension spanned by |Ii ⊗ Jj >; this approximation is obtained by minimizing In practice it is done by performing the Singular Value Decomposition (SVD) of the matrix C IJ (see Eq.(2.1)) which gives the components of the wave function in the A I ⊗ U J basis, and then selecting the largest singular values and the corresponding states. After we know the expansion of selected M states in terms of |Ii > for A and |Jj > for U , we must recalculate and store the 'primitive' matrix elements for these new states. Again it is done 'on the fly', expanding the matrix elements to newly added orbitals and contracting at the same time with coefficients produced by SVD procedure. For reference purposes in Appendix C we list all expanded matrix elements. These formulas are closely related to that of Appendix B.

III. NON-ORTHOGONAL ORBITALS
In this section we consider how the algorithms of the previous section can be extended to the case of non-orthogonal molecular orbitals. We start by describing general techniques to work with non-orthogonal orbitals. In principle there exist a general Löwdin formula that expresses the matrix elements of one-and two-electron operators between two determinants build with non-orthogonal orbitals as the sum of appropriate matrix determinants. It is this necessity to evaluate many matrix determinants that renders this approach completely unusable in practice. The practical way to work with non-orthogonal orbitals consists in introducing the so called bi-orthogonal (dual) orbitals [17,18].

A. General remarks on the eigenvalue problem with non-orthogonal bases
First of all, we cannot use the variational approach for the simple reason that we are unable to evaluate directly the diagonal energy expression which we demonstrate in the following section [20]. In order to compute the variational expression (Eq.(2.9)) for the energy, we must be able to express the vector |Ψ > in the biorthogonal orbitals; if we were able to find this expansion, then the computation follows more or less the standard root with insignificant changes since the one-electron and two-electron integrals loose their transformation properties. However in practice to find and use such expansion is an almost impossible task. In fact when the orbitals of A and U do overlap, a given state of the block A would become a linear combination (of nearly full CI size) of different A and U states. Therefore the whole strategy of DMRG becomes unapplicable. Still we will show that some modern non-variational approach can be used with more or less success. Let us introduce some notations that will make clearer the following discussion. Given a set {φ i } of non orthogonal orbitals, the dual or biorthogonal basisφ i verifies the condition: Let S denote the overlap matrix < φ i |φ j > ; given two one-particle states |a >= i a i |φ i >, |b >= i b i |φ i >, we can define two different scalar products: the "physical" scalar product which will be denoted by Dirac brackets, and a simple "numerical" scalar product defined by: which will be denoted by round brackets. Furthermore for any state |a >= i a i |φ i > we define the dual state |ã >= i a i |φ i >; therefore we have the identity: The concept of dual state can be generalized: a dual many particle state is the state built with the same CI (or DMRG) expansion coefficients but with {φ i } orbitals replaced by dual ones [21].
A nice property of the biorthogonal formalism is that the "dual" annihilation operators defined bỹ obey the usual anticommutation relations [17]: The one-particle and two-particle integrals are computed with the φ i wave functions on the right and the dualφ i on the left. For instance: In non-variational approach instead of looking for the minimum of energy in the trial subspace we project the energy eigenvalue equation onto the 'dual' subspace. To find the current approximation to energy and vector given the trial vector set we solve the following non-symmetric eigenvalue problem where H ij =<x i |Ĥ|x j >= (x i |Ĥ|x j ); because of biorthogonality <x i |x j >= (x i |x j ) = δ ij . If expansion set becomes reasonably full for the description of target state, this non-variation approach must produce reasonably good result, which is also demonstrated by some experience with non-orthogonal Full CI calculations.
In principle, if we want to set up a variational calculations with non orthogonal orbitals, we need the knowledge of the transformation operatorT from normal to dual vectors : x =Tx Then the variational energy becomes An iterative minimization procedure for this energy expression can be easily written down; however, it is almost impossible to follow the variational approach. Even the simple square < x|x > of the physical norm of a many particle state |x > which appears in the denominator of Eq.(3.10), is a huge expansion of determinants, since the scalar product of two simple configurations gives rise to a determinant; neglecting the spin indexes, and denoting by |0 > the vacuum, we have: As we shall see in the next paragraph, the same kind of difficulties exist for DMRG calculations.
Once the non symmetric Hamiltonian matrix H ij is constructed, we are left with the problem of using Davidson method or a similar direct meyhod. In paragraph 2.5 we have seen that the Davidson method can be considered as a minimization procedure of the expectation value (2.9). This procedure looses any meaning for a non symmetric matrix, since the minimum of (2.9) does not coincides with the lowest eigenvalue. However the Davidson method can be seen from a different angle, and can be used for solving the eigenvalue problem of non symmetric matrices. Of course the orthonormalization can be performed only with the "numerical" scalar product (3.4).
As well known for non symmetric matrices, there are other good methods that can be tried. The simplest one is perhaps the two-sided or non-symmetric Lanczos method that treats on an equal footing the Krylov sequences {x, Hx, H 2 x...} and {x, H T x, (H T ) 2 x...} generated by the matrix H and its transposed H T . Let us denote by V and W the two supspaces generated by the two Krylov sequences. It is easy to find two biorthogonal bases {v j } and {w i } in V and W , in such a way that the matrix (w i |Hv j ) is tridiagonal [19].
More complicated methods, like the two-sided Jacobi-Davidson method, which makes use of biorthogonal basis like the two-sided Lanczos, can also be considered [19] . The use of orthogonal bases is generally more numerically stable, while two-sided methods (and biorthogonal basis) may provide faster convergence.

B. Variational calculations with non-orthogonal orbitals in DMRG?
Here we show that variational calculations are completely impractical for non-orthogonal orbitals. We consider the overlap matrix just. The Hamiltonian matrix would be even much more complicated. So the problem is to evaluate Actually the problem is not just to evaluate Eq(22) but do it efficiently. For that we must factorize this expression to something containing only pair quantities, as it is anyway impossible to store and use M 4 quantities, corresponding to Eq (22). This means that we would like to have something like It is rather easy to show that such factorization is impossible in general case. We show that by considering the most simple example when all states I, J, K, L are determinants, and the number of electrons in I and K states is the same (we consider say α electrons, β being similar). Note that while for orthogonal orbitals the block structure is automatic, i.e. number of electrons in I and K (and also in J and L) MUST be the same to get non-zero overlap, it is NOT the case for non-orthogonal orbitals. Then the overlap, according to Löwdin formula is given by the determinant where S-matrices are buildee by the corresponding orbital overlaps. Assuming the general case when matrix S IK is non-degenerate and using the Gauss formula we have which has a complicated mixed structure and cannot be represented in the factorized form Eq(23).

C. Non-variational calculations
As we have seen that the only way is to perform non-variational calculations, the rest is more or less straightforward. To evaluate the matrix elements of Hamiltonian we use the same expressions as for orthogonal orbitals with minor modifications due to loss of integral symmetry. The full list of matrix elements to build the Hamiltonian matrix if given in Appendices A and B. Having Hamiltonian matrix elements (or better corresponding σ-routines), we solve non-symmetric eigenvalue problem as described before. The only remaining point is the reduction of states. We cannot calculate the reduced density for the same reason for which we cannot calculate theoverlap matrices. Actually, the concept of reduced density matrix is not clear in the case of non-orthogonal A and U . The matrix is not the reduced density matrix of one block and does not verifies Tr[ρ] = 1 .
However we can follow a simpler approach: we can abandon the idea of approximating the wave function using the physical norm < x|x > and we can use the following norm induced by the simple expression (3.4): Of course the norm defined by Eq. (27) is a true norm, so it is positive, it is zero only when x is zero ecc. Then we can find the 'best' approximation using this norm, which as before leads to SVD decomposition of the wave function coefficient matrix C IJ . Again when the orbital overlap matrix resembles unit, this should give results close to those obtained using the true physical norm. However the minimization with the "numerical" norm can give good results in more general cases.
Let us consider in more detail the practical side of the calculations. Let us define OP (∼↔ no ∼) the operator that interchanges normal dual orbitals in all the one-and two-electron integrals. Clearly we have: [22] OP (∼↔ no ∼)[< IJ|Ĥ|KL >] =< KL|Ĥ|IJ > This means that all 'primitive' matrix elements not containing integrals are the same for IK and KI. Therefore in practice we can keep the same number of cases (7) and just calculate two copies (one normal and one dual) of those primitive MEs that contain integrals, i.e. H IK for case 1, and f x for cases 3 and 4. Then we can use Eq.(28) to calculate inverse cases. The attention we must pay is that when we have LJ instead of JL referring to primitive MEs, we must take the dual copies. Considering that in practice the most of storage is taken by f pq primitive matrix elements of the case 1, which remains unchanged, the bi-orthogonal case has essentially the same memory-and-disk requirements as the standard orthogonal procedure.

APPENDIX A: NON-EXTENDED DMRG MATRIX ELEMENTS
We start by explaining some notations we use in the below tables. The phases as as follows: the A is created first, and U next. Within the state the α-(spin-up) electrons are created first, and β-(spin-down) next. So the |IJ > state where I refers to A and J to U , can be written as As to letters, I and K will refer to A block and J and L -to U . Orbital labels p, q, r, s, . . . will refer to A orbitals while x, y, z, . . . -to U . The states are ordered as follows in the program: I < K if: • n β I < n β K • n β I = n β K and n α I < n α K • n β I = n β K and n α I = n α K and symm(I) < symm(K) Here we list the matrix elements needed to be kept and the corresponding contributions to the matrix elements of total Hamiltonian. As noted before we consider only the case I ≤ K. The opposite case can be obtained from the modified 'Hermitian' property (3.19). For the specific case 1, when J > L in the below expressions, we can use the property:  MEs to keep < IJ|Ĥ|KL > J ≤ L? In this Appendix we give the matrix elements of Hamiltonian in the extended basis, i.e. when we have added the orbitals to all four states I, J, K, L. The states over those added orbitals are simple determinants, which we denote as patterns i, j, k, l. The number of cases thus grows enormously if we count for all possible occupations combinations for states I, J, K, L and patterns i, j, k, l. The fact that states over i, j, k, l are complete, i.e. are simple determinants, allows us to reformulate the problem by shifting everything to A block, i.e. we call the 'combined' pattern i + j as new i and k + l as new k. The only thing that remains after we have diagonalized Hamiltonian matrix is to cast this combined representation back to original split one. This is just the phase and indexing problem that is done very efficiently. Thus we have now orbitals added only to A block, as i and k patterns and thus the number of cases becomes reasonable. The phase convention is as follows: The orbital labels a, b, c, . . . will refer to orbitals of patterns i and k. In below expressions the sums over those indices a, b, c, . . . is assumed when not given explicitly. Table III lists all the matrix elements. For the explanation of quantities f and g see Table II. The 'new' case IiKk means the relations between number of electrons in Ii and Kk blocks, in the same way as before. The overscore on the case index mean inverted case, i.e. Ii > Kk.

APPENDIX C: EXTENDED 'PRIMITIVE' MATRIX ELEMENTS
Here we give the working expressions for 'primitive' matrix elements in the extended basis Ii, Kk. The notations are the same as in Appendix B. As cited in the text, these matrix elements are not stored. Instead they are calculated 'on the fly' and then contracted with vectors coefficients produced by SVD procedure to give final matrix elements for the reduced M state space. The full list of nonzero matrix elements is given in the Table IV.