Iterative linearized density matrix

The chapter is organized as follows The quantum-classical Liouville dynamics scheme is first outlined and a rigorous surface hopping trajectory algorithm for its implementation is presented. The iterative linearized density matrix propagation approach is then described and an approach for its implementation is presented. In the Model Simulations section the comparable performance of the two methods is documented for the generalized spin-boson model and numerical convergence issues are mentioned. In the Conclusions we review the perspectives of this study. 

In this section we present results using the two approaches described in the previous sections the Trotter factorized QCL (TQCL), and iterative linearized density matrix (ILDM) propagation schemes, to study the spin-boson model consisting of a two level system that is bi-linearly coupled to a bath with Mh harmonic modes. This popular model of a quantum system embedded in an environment is described by the following general hamiltonian ... 

Both the Fock matrix—through the density matrix—and the orbitals depend on the molecular orbital expansion coefficients. Thus, Equation 31 is not linear and must be solved iteratively. The procedure which does so is called the Self-Consistent Field... 

This iterative procedure depends linearly on the number of fragments and on the size of the target macromolecule M, as long as the parent molecules Mk are confined to some limited size. The storage of the information on the macromolecular basis set has relatively small computer memory requirements. The computation of the macromolecular electron density from this basis set information and the final macromolecular density matrix P(K) obtained from the finite iterative process (56) can rely on relation (32). As a consequence of the sparsity macromolecular density matrix P(AT), the computational task has linear computer time requirement with respect to the number of fragments, hence, with respect to the size of the target macromolecule M. 

The concept of purification is well known in the linear-scaling literature for one-particle theories like Hartree-Fock and density functional theory, where it denotes the iterative process by which an arbitrary one-particle density matrix is projected onto an idempotent 1-RDM [2,59-61]. An RDM is said to be pure A-representable if it arises from the integration of an Al-particle density matrix T T, where T (the preimage) is an Al-particle wavefiinction [3-5]. Any idempotent 1-RDM is N-representable with a unique Slater-determinant preimage. Within the linear-scaling literature the 1-RDM may be directly computed with unconstrained optimization, where iterative purification imposes the A-representabUity conditions [59-61]. Recently, we have shown that these methods for computing the 1 -RDM directly... 

Equation 5.82, a slight modification of Eq. 5.78, is the key equation in calculating the ab initio Fock matrix (you need memorize this equation only to the extent that the Fock matrix element consists of //corc, P, and the two-electron integrals). Each density matrix element Ptu represents the coefficients c for a particular pair of basis functions (f>, and (f> , summed over all the occupied MO s > /, (i 1,2,., n). We use the density matrix here just as a convenient way to express the Fock matrix elements, and to formulate the calculation of properties arising from electron distribution (Section 5.5.4), although there is far more to the density matrix concept [27]. Equation 5.82 enables the MO wavefunctions ij/ (which are linear combinations of the c s and s) and their energy levels e to be calculated by iterative diagonalization of the Fock matrix. 

Usually, the parent molecules M are confined to some limited size that allows rapid determination of the parent molecule density matrices within a conventional ab initio Hartree-Fock-Roothaan-Hall scheme, followed by the determination of the fragment density matrices and the assembly of the macro-molecular density matrix using the method described above. The entire iterative procedure depends linearly on the number of fragments, that is, on the size of the target macromolecule M. When compared to the conventional ab initio type methods of computer time requirements growing with the third or fourth power of the number of electrons, the linear scaling property of the ADMA method is advantageous. 

As described in Sec. 3.1, each Hartree-Fock iteration involves the construction of the Fock matrix for a given density matrix, followed by the diagonalization of the Fock matrix to generate a set of improved spin orbitals and thus an improved density matrix. Formally, the construction of the Fock matrix requires a number of operations proportional to K4, where K is the number of atoms (because the number of two-electron integrals scales as Al4). For large systems, however, this quartic scaling with K (i.e., with system size) can be reduced to linear by special techniques, as will now be discussed. 

As you see, self-consistency is achieved after about six iterations based on the criterion that the coefficients in the linear combination of equation 6.2, the starting LCAO-MO and the output coefficients, the LCAO-MO after the calculation cycle should be the same. You should devise a criterion based on the convergence of the density matrix elements and apply this to the calculation using the scfl worksheet. 

In iteration n the A matrix has dimension (n -1- 1) x (n -t 1), where n usually is less than 20. The coefficients c can be obtained by directly inverting the A matrix and multiplying it onto the b vector, i.e. in the subspace of the iterations the linear equations are solved by direct inversion , thus the name DIIS. Having obtained the coefficients that minimize the error function at iteration n, the same set of coefficients is used for generating an extrapolated Fock matrix (F ) at iteration n, which is used in place of F for generating the new density matrix. 

Using the variation condition one can derive a linear equation for Xq which, in turn, can be used to evaluate the first-order correction to the density matrix. This leads to an iterative scheme involving matrix manipulations in the local space. It is important to note that Kirtman s treatment describes not only inductive but delocalization effects, as well. A practical limitation is that one should work with orthogonal AO basis sets therefore the calculations have been done with semiempirical ZDO model Hamiltonians [91], or with explicitly Lowdin-orthogonalized basis sets in ab initio calculations. 

The system of Hartree-Fock equations (4.21) is nonhnear. To solve it, an iterative method is usually used. In the course of the pth iteration the electron-density matrix p(p)( r,r ) and hence the operators J and K are considered to be fixed. The system (4.21) then transforms into one linear equation with a fixed self-adjoint operator... 

The convergence properties of the density matrix-based equations, i.e., the number of iterations to converge P , are similar to the ones encountered for a solution in the MO space, so that the advantage of using sparse multiplications within the density-based approach allows us to reduce the scaling property of the computational effort in an efficient manner. In this way, NMR chemical shift calculations with linear-scaling effort become possible and systems with 1000 and more atoms can be treated at the HF or DFT level on today s computers.Extensions to other molecular properties can be formulated in a similar fashion. 

All matrix elements in the Newton-Raphson methods may be constructed from the one- and two-particle density matrices and transition density matrices. The linear equation solutions may be found using either direct methods or iterative methods. For large CSF expansions, such micro-iterative procedures may be used to advantage. If a micro-iterative procedure is chosen that requires only matrix-vector products to be formed, expansion-vector-dependent effective Hamiltonian operators and transition density matrices may be constructed for the efficient computation of these products. Sufficient information is included in the Newton-Raphson optimization procedures, through the gradient and Hessian elements, to ensure second-order convergence in some neighborhood of the final solution. 

---