An efficient algorithm for solving nonlinear equations with a minimal number of trial vectors: applications to atomic-orbital based coupled-cluster theory

The conjugate residual with optimal trial vectors (CROP) algorithm is developed. In this algorithm, the optimal trial vectors of the iterations are used as basis vectors in the iterative subspace. For linear equations and nonlinear equations with a small-to-medium nonlinearity, the iterative subspace may be truncated to a three-dimensional subspace with no or little loss of convergence rate, and the norm of the residual decreases in each iteration. The efficiency of the algorithm is demonstrated by solving the equations of coupled-cluster theory with single and double excitations in the atomic orbital basis. By performing calculations on H(2)O with various bond lengths, the algorithm is tested for varying degrees of nonlinearity. In general, the CROP algorithm with a three-dimensional subspace exhibits fast and stable convergence and outperforms the standard direct inversion in iterative subspace method.     

Development and application of the analytical energy gradient for the normalized elimination of the small component method

The analytical energy gradient of the normalized elimination of the small component (NESC) method is derived for the first time and implemented for the routine calculation of NESC geometries and other first order molecular properties. Essential for the derivation is the correct calculation of the transformation matrix U relating the small component to the pseudolarge component of the wavefunction. The exact form of ?U/?λ is derived and its contribution to the analytical energy gradient is investigated. The influence of a finite nucleus model and that of the picture change is determined. Different ways of speeding up the calculation of the NESC gradient are tested. It is shown that first order properties can routinely be calculated in combination with Hartree-Fock, density functional theory (DFT), coupled cluster theory, or any electron correlation corrected quantum chemical method, provided the NESC Hamiltonian is determined in an efficient, but nevertheless accurate way. The general applicability of the analytical NESC gradient is demonstrated by benchmark calculations for NESC/CCSD (coupled cluster with all single and double excitation) and NESC/DFT involving up to 800 basis functions.     

New algorithms for the computation of column dynamics of multicomponent liquid phase adsorption

New and efficient numerical algorithms were developed for simulating column dynamics of multicomponent liquid phase adsorption. Simple and realistic models are used for the simulation. Langmuir form of isotherm and linear driving force rate expressions are employed in the model equations. Algorithms were formulated for three different rate control mechanisms, namely, film diffusion control, particle diffusion control and combined film and particle diffusion control. The algorithms derived are explicit with the exception of the requirement of solving a nonlinear equation in one single variable which is the concentration of a reference species. Thus the tedious iterative calculation procedure for solving simultaneous nonlinear equations in a multicomponent fixed bed system is avoided. Example calculations indicated very good numerical accuracy as verified from an independent check by means of an overall mass balance.     

Full S matrix calculation via a single real-symmetric Lanczos recursion: the Lanczos artificial boundary inhomogeneity method

We present an efficient and robust method for the calculation of all S matrix elements (elastic, inelastic, and reactive) over an arbitrary energy range from a single real-symmetric Lanczos recursion. Our new method transforms the fundamental equations associated with Light's artificial boundary inhomogeneity approach from the primary representation (original grid or basis representation of the Hamiltonian or its function) into a single tridiagonal Lanczos representation, thereby affording an iterative version of the original algorithm with greatly superior scaling properties. The method has important advantages over existing iterative quantum dynamical scattering methods: (a) the numerically intensive matrix propagation proceeds with real symmetric algebra, which is inherently more stable than its complex symmetric counterpart; (b) no complex absorbing potential or real damping operator is required, saving much of the exterior grid space which is commonly needed to support these operators and also removing the associated parameter dependence. Test calculations are presented for the collinear H+H(2) reaction, revealing excellent performance characteristics.     

Manager-worker-based model for the parallelization of quantum Monte Carlo on heterogeneous and homogeneous networks

A manager-worker-based parallelization algorithm for Quantum Monte Carlo (QMC-MW) is presented and compared with the pure iterative parallelization algorithm, which is in common use. The new manager-worker algorithm performs automatic load balancing, allowing it to perform near the theoretical maximal speed even on heterogeneous parallel computers. Furthermore, the new algorithm performs as well as the pure iterative algorithm on homogeneous parallel computers. When combined with the dynamic distributable decorrelation algorithm (DDDA) [Feldmann et al., J Comput Chem 28, 2309 (2007)], the new manager-worker algorithm allows QMC calculations to be terminated at a prespecified level of convergence rather than upon a prespecified number of steps (the common practice). This allows a guaranteed level of precision at the least cost. Additionally, we show (by both analytic derivation and experimental verification) that standard QMC implementations are not "perfectly parallel" as is often claimed.     

Cholesky decomposition of the two-electron integral matrix in electronic structure calculations

A standard Cholesky decomposition of the two-electron integral matrix leads to integral tables which have a huge number of very small elements. By neglecting these small elements, it is demonstrated that the recursive part of the Cholesky algorithm is no longer a bottleneck in the procedure. It is shown that a very efficient algorithm can be constructed when family type basis sets are adopted. For subsequent calculations, it is argued that two-electron integrals represented by Cholesky integral tables have the same potential for simplifications as density fitting. Compared to density fitting, a Cholesky decomposition of the two-electron matrix is not subjected to the problem of defining an auxiliary basis for obtaining a fixed accuracy in a calculation since the accuracy simply derives from the choice of a threshold for the decomposition procedure. A particularly robust algorithm for solving the restricted Hartree-Fock (RHF) equations can be speeded up if one has access to an ordered set of integral tables. In a test calculation on a linear chain of beryllium atoms, the advocated RHF algorithm nicely converged, but where the standard direct inversion in iterative space method converged very slowly to an excited state.     

Linear scaling local correlation approach for solving the coupled cluster equations of large systems

A linear scaling local correlation approach is proposed for approximately solving the coupled cluster doubles (CCD) equations of large systems in a basis of orthogonal localized molecular orbitals (LMOs). By restricting double excitations from spatially close occupied LMOs into their associated virtual LMOs, the number of significant excitation amplitudes scales only linearly with molecular size in large molecules. Significant amplitudes are obtained to a very good approximation by solving the CCD equations of various subsystems, each of which is made up of a cluster associated with the orbital indices of a subset of significant amplitudes and the local environmental domain of the cluster. The combined effect of these two approximations leads to a linear scaling algorithm for large systems. By using typical thresholds, which are designed to target an energy accuracy, our numerical calculations for a wide range of molecules using the 6-31G or 6-31G* basis set demonstrate that the present local correlation approach recovers more than 98.5% of the conventional CCD correlation energy.     

Relativistic calculation of nuclear magnetic shielding using normalized elimination of the small component

The normalized elimination of the small component (NESC) theory, recently proposed by Filatov and Cremer, is extended to include magnetic interactions and applied to the calculation of the nuclear magnetic shielding in HX (X=F, Cl, Br, I) systems. The NESC calculations are performed at the levels of the zeroth-order regular approximation (ZORA) and the second-order regular approximation (SORA). The calculations show that the NESC-ZORA results are very close to the NESC-SORA results, except for the shielding of the I nucleus. Both the NESC-ZORA and NESC-SORA calculations yield very similar results to the previously reported values obtained using the relativistic infinite-order two-component coupled Hartree-Fock method. The difference between NESC-ZORA and NESC-SORA results is significant for the shieldings of iodine.     

An efficient implementation of the "cluster-in-molecule" approach for local electron correlation calculations

An efficient implementation of the "cluster-in-molecule" (CIM) approach is presented for performing local electron correlation calculations in a basis of orthogonal occupied and virtual localized molecular orbitals (LMOs). The main idea of this approach is that significant excitation amplitudes can be approximately obtained by solving the coupled cluster (or Moller-Plesset perturbation theory) equations of a series of "clusters," each of which contains a subset of occupied and virtual LMOs. In the present implementation, we have proposed a simple approach for constructing virtual LMOs of clusters, and new ways of constructing clusters and extracting the correlation contributions from calculations on clusters, which are more efficient than those suggested in the original work. More importantly, linear scaling of computational time of the CIM approach is achieved by evaluating the transformed two-electron integrals over LMOs using simple truncation techniques in limited operations (independent of the molecular size). With typical thresholds, for a variety of molecules our test calculations demonstrate that more than 99% of the conventional MP2 or coupled cluster with doubles correlation energies can be recovered in the present CIM approach.     

Wavelength selection for multivariate calibration using a genetic algorithm: a novel initialization strategy

Genetic algorithms and other procedures mimicking natural processes are being increasingly used for variable selection, to improve the predictive ability of partial least-squares multivariate calibration. Two issues are critical for the success of genetic algorithms: initialization (setting the first candidates for solving the problem at hand) and overfitting (the tendency to produce excellent results when training, but poor predictions toward fresh samples). A new procedure is presented for sensor selection problems, involving iterative reinitialization based on a statistical analysis of the included sensors. It is shown to give excellent results without the requirement of preparing independent test data sets. Monte Carlo simulations using a theoretical three-component example illustrate how partial least-squares regression greatly benefits from variable selection when the analyte of interest is diluted, and how the new initialization method compares with other strategies. The new genetic algorithm was applied to five experimental data sets. The target parameters were the concentrations of diluted analytes in four pharmaceutical mixtures studied by UV-visible spectrophotometry and the octane number in gasolines analyzed by near-infrared spectroscopy.     

On a partial differential equation method for determining the free energies and coexisting phase compositions of ternary mixtures from light scattering data

In this paper we present a method for determining the free energies of ternary mixtures from light scattering data. We use an approximation that is appropriate for liquid mixtures, which we formulate as a second-order nonlinear partial differential equation. This partial differential equation (PDE) relates the Hessian of the intensive free energy to the efficiency of light scattering in the forward direction. This basic equation applies in regions of the phase diagram in which the mixtures are thermodynamically stable. In regions in which the mixtures are unstable or metastable, the appropriate PDE is the nonlinear equation for the convex hull. We formulate this equation along with continuity conditions for the transition between the two equations at cloud point loci. We show how to discretize this problem to obtain a finite-difference approximation to it, and we present an iterative method for solving the discretized problem. We present the results of calculations that were done with a computer program that implements our method. These calculations show that our method is capable of reconstructing test free energy functions from simulated light scattering data. If the cloud point loci are known, the method also finds the tie lines and tie triangles that describe thermodynamic equilibrium between two or among three liquid phases. A robust method for solving this PDE problem, such as the one presented here, can be a basis for optical, noninvasive means of characterizing the thermodynamics of multicomponent mixtures.     

KSSOLV-GPU: an Efficient GPU-Enabled MATLAB Toolbox for Solving the Kohn-Sham Equations within Density Functional Theory in Plane-Wave Basis Set?

KSSOLV (Kohn-Sham Solver) is a MATLAB (Matrix Laboratory) toolbox for solving the Kohn-Sham density functional theory (KS-DFT) with the plane-wave basis set. In the KS-DFT calculations, the most expensive part is commonly the diagonalization of Kohn-Sham Hamiltonian in the self-consistent field (SCF) scheme. To enable a personal computer to perform medium-sized KS-DFT calculations that contain hundreds of atoms, we present a hybrid CPU-GPU implementation to accelerate the iterative diagonalization algorithms implemented in KSSOLV by using the MATLAB built-in Parallel Computing Toolbox. We compare the performance of KSSOLV-GPU on three types of GPU, including RTX3090, V100, and A100, with conventional CPU implementation of KSSOLV respectively and numerical results demonstrate that hybrid CPU-GPU implementation can achieve a speedup of about 10 times compared with sequential CPU calculations for bulk silicon systems containing up to 128 atoms.     

Globally convergent trust-region methods for self-consistent field electronic structure calculations

As far as more complex systems are being accessible for quantum chemical calculations, the reliability of the algorithms used becomes increasingly important. Trust-region strategies comprise a large family of optimization algorithms that incorporates both robustness and applicability for a great variety of problems. The objective of this work is to provide a basic algorithm and an adequate theoretical framework for the application of globally convergent trust-region methods to electronic structure calculations. Closed shell restricted Hartree-Fock calculations are addressed as finite-dimensional nonlinear programming problems with weighted orthogonality constraints. A Levenberg-Marquardt-like modification of a trust-region algorithm for constrained optimization is developed for solving this problem. It is proved that this algorithm is globally convergent. The subproblems that ensure global convergence are easy-to-compute projections and are dependent only on the structure of the constraints, thus being extendable to other problems. Numerical experiments are presented, which confirm the theoretical predictions. The structure of the algorithm is such that accelerations can be easily associated without affecting the convergence properties.     

Quantum and electromagnetic propagation with the conjugate symmetric Lanczos method

The conjugate symmetric Lanczos (CSL) method is introduced for the solution of the time-dependent Schrodinger equation. This remarkably simple and efficient time-domain algorithm is a low-order polynomial expansion of the quantum propagator for time-independent Hamiltonians and derives from the time-reversal symmetry of the Schrodinger equation. The CSL algorithm gives forward solutions by simply complex conjugating backward polynomial expansion coefficients. Interestingly, the expansion coefficients are the same for each uniform time step, a fact that is only spoiled by basis incompleteness and finite precision. This is true for the Krylov basis and, with further investigation, is also found to be true for the Lanczos basis, important for efficient orthogonal projection-based algorithms. The CSL method errors roughly track those of the short iterative Lanczos method while requiring fewer matrix-vector products than the Chebyshev method. With the CSL method, only a few vectors need to be stored at a time, there is no need to estimate the Hamiltonian spectral range, and only matrix-vector and vector-vector products are required. Applications using localized wavelet bases are made to harmonic oscillator and anharmonic Morse oscillator systems as well as electrodynamic pulse propagation using the Hamiltonian form of Maxwell's equations. For gold with a Drude dielectric function, the latter is non-Hermitian, requiring consideration of corrections to the CSL algorithm.     

New algorithms for iterative matrix‐free eigensolvers in quantum chemistry

New algorithms for iterative diagonalization procedures that solve for a small set of eigen‐states of a large matrix are described. The performance of the algorithms is illustrated by calculations of low and high‐lying ionized and electronically excited states using equation‐of‐motion coupled‐cluster methods with single and double substitutions (EOM‐IP‐CCSD and EOM‐EE‐CCSD). We present two algorithms suitable for calculating excited states that are close to a specified energy shift (interior eigenvalues). One solver is based on the Davidson algorithm, a diagonalization procedure commonly used in quantum‐chemical calculations. The second is a recently developed solver, called the “Generalized Preconditioned Locally Harmonic Residual (GPLHR) method.” We also present a modification of the Davidson procedure that allows one to solve for a specific transition. The details of the algorithms, their computational scaling, and memory requirements are described. The new algorithms are implemented within the EOM‐CC suite of methods in the Q‐Chem electronic structure program. © 2014 Wiley Periodicals, Inc.     

Connection between the regular approximation and the normalized elimination of the small component in relativistic quantum theory

The regular approximation to the normalized elimination of the small component (NESC) in the modified Dirac equation has been developed and presented in matrix form. The matrix form of the infinite-order regular approximation (IORA) expressions, obtained in [Filatov and Cremer, J. Chem. Phys. 118, 6741 (2003)] using the resolution of the identity, is the exact matrix representation and corresponds to the zeroth-order regular approximation to NESC (NESC-ZORA). Because IORA (=NESC-ZORA) is a variationally stable method, it was used as a suitable starting point for the development of the second-order regular approximation to NESC (NESC-SORA). As shown for hydrogenlike ions, NESC-SORA energies are closer to the exact Dirac energies than the energies from the fifth-order Douglas-Kroll approximation, which is much more computationally demanding than NESC-SORA. For the application of IORA (=NESC-ZORA) and NESC-SORA to many-electron systems, the number of the two-electron integrals that need to be evaluated (identical to the number of the two-electron integrals of a full Dirac-Hartree-Fock calculation) was drastically reduced by using the resolution of the identity technique. An approximation was derived, which requires only the two-electron integrals of a nonrelativistic calculation. The accuracy of this approach was demonstrated for heliumlike ions. The total energy based on the approximate integrals deviates from the energy calculated with the exact integrals by less than 5 x 10(-9) hartree units. NESC-ZORA and NESC-SORA can easily be implemented in any nonrelativistic quantum chemical program. Their application is comparable in cost with that of nonrelativistic methods. The methods can be run with density functional theory and any wave function method. NESC-SORA has the advantage that it does not imply a picture change.     

Gradient-based direct normal-mode analysis

A formulation of a direct, iterative method for obtaining the lowest eigenvalues and eigenvectors of a Hessian matrix is presented. Similar to the iterative schemes in electronic structure configuration interaction calculations (methods due to Lanczos, Davidson, and others), the mass-weighted Hessian matrix K is not constructed explicitly; instead, its operation on a basis vector (a direction in the 3N Cartesian configuration space of the atoms) is computed based on the principles of dynamical equations of motion. By noting that the time derivative of the gradient vector in the harmonic force field is related to the particles' momenta via dg/dt = Kp, a Hessian-vector product can be computed on the fly by finite differentiation of the gradient along the direction specified by the p vector. Thus, only two evaluations of the gradient are required per Davidson-like iteration per root, which leads to a linear scaling behavior of the computational effort with the number of atoms (provided that the evaluation of the gradient scales linearly). Preliminary results are presented for a 27,000-atom 4He nanodroplet.     

Approximate treatment of higher excitations in coupled-cluster theory

The possibilities for the approximate treatment of higher excitations in coupled-cluster (CC) theory are discussed. Potential routes for the generalization of corresponding approximations to lower-level CC methods are analyzed for higher excitations. A general string-based algorithm is presented for the evaluation of the special contractions appearing in the equations specific to those approximate CC models. It is demonstrated that several iterative and noniterative approximations to higher excitations can be efficiently implemented with the aid of our algorithm and that the coding effort is mostly reduced to the generation of the corresponding formulas. The performance of the proposed and implemented methods for total energies is assessed with special regard to quadruple and pentuple excitations. The applicability of our approach is illustrated by benchmark calculations for the butadiene molecule. Our results demonstrate that the proposed algorithm enables us to consider the effect of quadruple excitations for molecular systems consisting of up to 10-12 atoms.