Density functional theory for molecular and periodic systems using density fitting and continuous fast multipole method: Analytical gradients

A full implementation of analytical energy gradients for molecular and periodic systems is reported in the TURBOMOLE program package within the framework of Kohn–Sham density functional theory using Gaussian‐type orbitals as basis functions. Its key component is a combination of density fitting (DF) approximation and continuous fast multipole method (CFMM) that allows for an efficient calculation of the Coulomb energy gradient. For exchange‐correlation part the hierarchical numerical integration scheme (Burow and Sierka, Journal of Chemical Theory and Computation 2011, 7, 3097) is extended to energy gradients. Computational efficiency and asymptotic O(N) scaling behavior of the implementation is demonstrated for various molecular and periodic model systems, with the largest unit cell of hematite containing 640 atoms and 19,072 basis functions. The overall computational effort of energy gradient is comparable to that of the Kohn–Sham matrix formation. © 2016 Wiley Periodicals, Inc.     

Real-time time-dependent density functional theory using density fitting and the continuous fast multipole method

An implementation of real-time time-dependent density functional theory (RT-TDDFT) within the TURBOMOLE program package is reported using Gaussian-type orbitals as basis functions, second and fourth order Magnus propagator, and the self-consistent field as well as the predictor–corrector time integration schemes. The Coulomb contribution to the Kohn–Sham matrix is calculated combining density fitting approximation and the continuous fast multipole method. Performance of the implementation is benchmarked for molecular systems with different sizes and dimensionalities. For linear alkane chains, the wall time for density matrix time propagation step is comparable to the Kohn-Sham (KS) matrix construction. However, for larger two- and three-dimensional molecules, with up to about 5,000 basis functions, the computational effort of RT-TDDFT calculations is dominated by the KS matrix evaluation. In addition, the maximum time step is evaluated using a set of small molecules of different polarities. The photoabsorption spectra of several molecular systems calculated using RT-TDDFT are compared to those obtained using linear response time-dependent density functional theory and coupled cluster methods.     

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.     

Use of promolecular ASA density functions as a general algorithm to obtain starting MO in SCF calculations

Atomic shell approximation (ASA) constitutes a way to fit first‐order density functions to a linear combination of spherical functions. The ASA fitting method makes use of positive definite expansion coefficients to ensure appropriate probability distribution features. The ASA electron density is sufficiently accurate for the practical implementation of quantum similarity measures, as was proved in previous published work. Here, a new application of the ASA density formalism is analyzed, and employed to obtain an initial guess of the density matrix for SCF procedures. The number of cycles needed to assess the convergence criterion in electronic energy calculations appears comparable to or less than those obtained by other means. Several molecular structures of different classes, including organic systems and metal complexes, were chosen as representative test cases. In addition, an ASA basis set for atoms Sc‐Kr fitted to an ab initio 6‐311G basis set is also presented. © 2001 John Wiley & Sons, Inc. Int J Quantum Chem, 2001     

Long‐range analysis of density fitting in extended systems

Density fitting scheme is analyzed for the Coulomb problem in extended systems from the correctness of long‐range behavior point of view. We show that for the correct cancellation of divergent long‐range Coulomb terms it is crucial for the density fitting scheme to reproduce the overlap matrix exactly. It is demonstrated that from all possible fitting metric choices the Coulomb metric is the only one which inherently preserves the overlap matrix for infinite systems with translational periodicity. Moreover, we show that by a small additional effort any non‐Coulomb metric fit can be made overlap‐preserving as well. The problem is analyzed for both ordinary and Poisson basis set choices. © 2008 Wiley Periodicals, Inc. Int J Quantum Chem, 2008     

Attractive electron–electron interactions within robust local fitting approximations

An analysis of Dunlap's robust fitting approach reveals that the resulting two‐electron integral matrix is not manifestly positive semidefinite when local fitting domains or non‐Coulomb fitting metrics are used. We present a highly local approximate method for evaluating four‐center two‐electron integrals based on the resolution‐of‐the‐identity (RI) approximation and apply it to the construction of the Coulomb and exchange contributions to the Fock matrix. In this pair‐atomic resolution‐of‐the‐identity (PARI) approach, atomic‐orbital (AO) products are expanded in auxiliary functions centered on the two atoms associated with each product. Numerical tests indicate that in 1% or less of all Hartree–Fock and Kohn–Sham calculations, the indefinite integral matrix causes nonconvergence in the self‐consistent‐field iterations. In these cases, the two‐electron contribution to the total energy becomes negative, meaning that the electronic interaction is effectively attractive, and the total energy is dramatically lower than that obtained with exact integrals. In the vast majority of our test cases, however, the indefiniteness does not interfere with convergence. The total energy accuracy is comparable to that of the standard Coulomb‐metric RI method. The speed‐up compared with conventional algorithms is similar to the RI method for Coulomb contributions; exchange contributions are accelerated by a factor of up to eight with a triple‐zeta quality basis set. A positive semidefinite integral matrix is recovered within PARI by introducing local auxiliary basis functions spanning the full AO product space, as may be achieved by using Cholesky‐decomposition techniques. Local completion, however, slows down the algorithm to a level comparable with or below conventional calculations. © 2013 Wiley Periodicals, Inc.     

Analytical nuclear excited‐state gradients for the Tamm–Dancoff approximation using uncoupled frozen‐density embedding

We report the derivation and implementation of analytical nuclear gradients for excited states using time‐dependent density functional theory using the Tamm–Dancoff approximation combined with uncoupled frozen‐density embedding using density fitting. Explicit equations are presented and discussed. The implementation is able to treat singlet as well as triplet states and functionals using the local density approximation, the generalized gradient approximation, combinations with Hartree–Fock exchange (hybrids), and range‐separated functionals such as CAM‐B3LYP. The new method is benchmarked against supermolecule calculations in two case studies: The solvatochromic shift of the (vertical) fluorescence energy of 4‐aminophthalimide on solvation, and the first local excitation of the benzonitrile dimer. Whereas for the 4‐aminophthalimide–water complex deviations of about 0.2 eV are obtained to supermolecular calculations, for the benzonitrile dimer the maximum error for adiabatic excitation energies is below 0.01 eV due to a weak coupling of the subsystems. © 2017 Wiley Periodicals, Inc.     

Parallel implementation of divide-and-conquer semiempirical quantum chemistry calculations

We have implemented a parallel divide-and-conquer method for semiempirical quantum mechanical calculations. The standard message passing library, the message passing interface (MPI), was used. In this parallel version, the memory needed to store the Fock and density matrix elements is distributed among the processors. This memory distribution solves the problem of demanding requirement of memory for very large molecules. While the parallel calculation for construction of matrix elements is straightforward, the parallel calculation of Fock matrix diagonalization is achieved via the divide-and-conquer method. Geometry optimization is also implemented with parallel gradient calculations. The code has been tested on a Cray T3E parallel computer, and impressive speedup of calculations has been achieved. Our results indicate that the divide-and-conquer method is efficient for parallel implementation. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1101–1109, 1998     

用密度梯度理论研究缔合流体的汽液成核性质

在密度梯度展开的基础上,将影响参数k 表达成温度的函数,建立了一个适用于均相和非均相缔合流体的状态方程。应用流体的蒸汽压和液相密度实验数据关联分子参数。在密度梯度理论的框架下,计算了水,重水,甲醇,乙醇,正丙醇,正丁醇,正戊醇和正己醇的成核速率并与实验数据进行了对比,计算结果令人满意。结果表明,密度梯度理论与密度泛函理论一样,可研究液核的结构和性质,但通过调整影响参数k, 可获得更为准确的成核速率。     

Density functional theory for molecular and periodic systems using density fitting and continuous fast multipole method: Stress tensor

A full implementation of the analytical stress tensor for periodic systems is reported in the TURBOMOLE program package within the framework of Kohn–Sham density functional theory using Gaussian-type orbitals as basis functions. It is the extension of the implementation of analytical energy gradients (Lazarski et al., Journal of Computational Chemistry 2016, 37, 2518–2526) to the stress tensor for the purpose of optimization of lattice vectors. Its key component is the efficient calculation of the Coulomb contribution by combining density fitting approximation and continuous fast multipole method. For the exchange-correlation (XC) part the hierarchical numerical integration scheme (Burow and Sierka, Journal of Chemical Theory and Computation 2011, 7, 3097–3104) is extended to XC weight derivatives and stress tensor. The computational efficiency and favorable scaling behavior of the stress tensor implementation are demonstrated for various model systems. The overall computational effort for energy gradient and stress tensor for the largest systems investigated is shown to be at most two and a half times the computational effort for the Kohn–Sham matrix formation. © 2019 Wiley Periodicals, Inc.     

Fast and accurate Coulomb calculation with Gaussian functions

Coulomb interaction is one of the major time-consuming components in a density functional theory (DFT) calculation. In the last decade, dramatic progresses have been made to improve the efficiency of Coulomb calculation, including continuous fast multipole method (CFMM) and J-engine method, all developed first inside Q-Chem. The most recent development is the advent of Fourier transform Coulomb method developed by Fusti-Molnar and Pulay, and an improved version of the method has been recently implemented in Q-Chem. It replaces the least efficient part of the previous Coulomb methods with an accurate numerical integration scheme that scales in O(N2) instead of O(N4) with the basis size. The result is a much smaller slope in the linear scaling with respect to the molecular size and we will demonstrate through a series of benchmark calculations that it speeds up the calculation of Coulomb energy by several folds over the efficient existing code, i.e., the combination of CFMM and J-engine, without loss of accuracy. Furthermore, we will show that it is complementary to the latter and together the three methods offer the best performance for Coulomb part of DFT calculations, making the DFT calculations affordable for very large systems involving thousands of basis functions.     

Exploiting graphical processing units to enable quantum chemistry calculation of large solvated molecules with conductor-like polarizable continuum models

The conductor-like polarizable continuum model (C-PCM) with switching/Gaussian smooth discretization is a widely used implicit solvation model in quantum chemistry. We have previously implemented C-PCM solvation for Hartree-Fock (HF) and density functional theory (DFT) on graphical processing units (GPUs), enabling the quantum mechanical treatment of large solvated biomolecules. Here, we first propose a GPU-based algorithm for the PCM conjugate gradient linear solver that greatly improves the performance for very large molecules. The overhead for PCM-related evaluations now consumes less than 15% of the total runtime for DFT calculations on large molecules. Second, we demonstrate that our algorithms tailored for ground state C-PCM are transferable to excited state properties. Using a single GPU, our method evaluates the analytic gradient of the linear response PCM time-dependent density functional theory energy up to 80× faster than a conventional central processing unit (CPU)-based implementation. In addition, our C-PCM algorithms are transferable to other methods that require electrostatic potential (ESP) evaluations. For example, we achieve speed-ups of up to 130× for restricted ESP-based atomic charge evaluations, when compared to CPU-based codes. We also summarize and compare the different PCM cavity discretization schemes used in some popular quantum chemistry packages as a reference for both users and developers.     

Quantum crystallography: A perspective

Extraction of the complete quantum mechanics from X‐ray scattering data is the ultimate goal of quantum crystallography. This article delivers a perspective for that possibility. It is desirable to have a method for the conversion of X‐ray diffraction data into an electron density that reflects the antisymmetry of an N‐electron wave function. A formalism for this was developed early on for the determination of a constrained idempotent one‐body density matrix. The formalism ensures pure‐state N‐representability in the single determinant sense. Applications to crystals show that quantum mechanical density matrices of large molecules can be extracted from X‐ray scattering data by implementing a fragmentation method termed the kernel energy method (KEM). It is shown how KEM can be used within the context of quantum crystallography to derive quantum mechanical properties of biological molecules (with low data‐to‐parameters ratio). © 2017 Wiley Periodicals, Inc.     

Hartree-Fock calculations with linearly scaling memory usage

We present an implementation of a set of algorithms for performing Hartree-Fock calculations with resource requirements in terms of both time and memory directly proportional to the system size. In particular, a way of directly computing the Hartree-Fock exchange matrix in sparse form is described which gives only small addressing overhead. Linear scaling in both time and memory is demonstrated in benchmark calculations for system sizes up to 11 650 atoms and 67 204 Gaussian basis functions on a single computer with 32 Gbytes of memory. The sparsity of overlap, Fock, and density matrices as well as band gaps are also shown for a wide range of system sizes, for both linear and three-dimensional systems.     

Generic approach to the method development of intact protein separations using hydrophobic interaction chromatography

We describe a liquid chromatography method development approach for the separation of intact proteins using hydrophobic interaction chromatography. First, protein retention was determined as function of the salt concentration by isocratic measurements and modeled using linear regression. The error between measured and predicted retention factors was studied while varying gradient time (between 15 and 120 min) and gradient starting conditions, and ranged between 2 and 15%. To reduce the time needed to develop optimized gradient methods for hydrophobic interaction chromatography separations, retention‐time estimations were also assessed based on two gradient scouting runs, resulting in significantly improved retention‐time predictions (average error < 2.5%) when varying gradient time. When starting the scouting gradient at lower salt concentrations (stronger eluent), retention time prediction became inaccurate in contrast to predictions based on isocratic runs. Application of three scouting runs and a nonlinear model, incorporating the effects of gradient duration and mobile‐phase composition at the start of the gradient, provides accurate results (improved fitting compared to the linear solvent‐strength model) with an average error of 1.0% and maximum deviation of –8.3%. Finally, gradient scouting runs and retention‐time modeling have been applied for the optimization of a critical‐pair protein isoform separation encountered in a biotechnological sample.     

Accelerating Kohn‐Sham response theory using density fitting and the auxiliary‐density‐matrix method

An extension of the formulation of the atomic‐orbital‐based response theory of Larsen et al., JCP 113, 8909 (2000) is presented. This new framework has been implemented in LSDalton and allows for the use of Kohn‐Sham density‐functional theory with approximate treatment of the Coulomb and Exchange contributions to the response equations via the popular resolution‐of‐the‐identity approximation as well as the auxiliary‐density matrix method (ADMM). We present benchmark calculations of ground‐state energies as well as the linear and quadratic response properties: vertical excitation energies, polarizabilities, and hyperpolarizabilities. The quality of these approximations in a range of basis sets is assessed against reference calculations in a large aug‐pcseg‐4 basis. Our results confirm that density fitting of the Coulomb contribution can be used without hesitation for all the studied properties. The ADMM treatment of exchange is shown to yield high accuracy for ground‐state and excitation energies, whereas for polarizabilities and hyperpolarizabilities the performance gain comes at a cost of accuracy. Excitation energies of a tetrameric model consisting of units of the P700 special pigment of photosystem I have been studied to demonstrate the applicability of the code for a large system.     

Wave‐function frozen‐density embedding: Approximate analytical nuclear ground‐state gradients

We report the derivation of approximate analytical nuclear ground‐state uncoupled frozen density embedding (FDEu) gradients for the resolution of identity (RI) variant of the second‐order approximate coupled cluster singles and doubles (RICC2) as well as density functional theory (DFT), and an efficient implementation thereof in the KOALA program. In order to guarantee a computationally efficient treatment, those gradient terms are neglected which would require the exchange of orbital information. This approach allows for geometry optimizations of single molecules surrounded by numerous molecules with fixed nuclei at RICC2‐in‐RICC2, RICC2‐in‐DFT, and DFT‐in‐DFT FDE level of theory using a dispersion correction, required due to the DFT‐based treatment of the interaction in FDE theory. Accuracy and applicability are assessed by the example of two case studies: (a) the Watson‐Crick pair adenine‐thymine, for which the optimized structures exhibit a maximum error of about 0.08 Å for our best scheme compared to supermolecular reference calculations, (b) carbon monoxide on a magnesium oxide surface model, for which the error amount up to 0.1 Å for our best scheme. Efficiency is demonstrated by successively including environment molecules and comparing to an optimized conventional supermolecular implementation, showing that the method is able to outperform conventional RICC2 schemes already with a rather small number of environment molecules, gaining significant speed up in computation time. © 2016 Wiley Periodicals, Inc.     

The Varigrad and similar apparatus in gradient liquid chromatography: A systematic computational approach. Part 1: Theoretical development

Different methods for producing a pre-assigned gradient in liquid chromatography are presented. A systematic approach to the calculations involved is given. The derived equations are mostly Poisson and Poisson summation distributions which are tabulated in the literature. The incremental method of gradient elution developed by Scott has been modified. In the modified apparatus two mixing chambers are used instead of one. This leads to an appreciable decrease in the number of reservoirs needed for the same precision in fitting a desired gradient. The application of the derived equations, together with other Varigrad modifications, will be given in Part 2 of this paper (to be published shortly in this Journal).     

Theoretical description of unconstrained thermally induced shape‐memory recovery in crosslinked polyethylenes

A new theoretical approach based on the modified three‐element Eyring‐Halsey model was developed for the derivation of an equation describing the thermally induced recovery of predeformed and crystallized crosslinked polymers. The proposed approach takes into account the influence of crystallizable covalent network and of entangled slipped molecular chains. Modeling of thermally induced shape‐memory (SM) recovery strain and SM recovery rate detected at constant heating rate has been successfully performed for nearly linear and two short‐chain branched polyethylenes, which were crosslinked by peroxide. The values of material constants determined by fitting agree with the estimations existing in literature. Fitting results have shown that increase of degree of branching and crosslink density accompanied with reducing crystallinity results in increasing contribution of the entangled slipped chains to the total stored SM strain. The physical sense of main fitting parameters and their dependences on the material constants such as crystallinity are discussed. © 2014 Wiley Periodicals, Inc. J. Polym. Sci., Part B: Polym. Phys. 2014 , 52, 815–822     

Carrier ampholyte‐free isoelectric focusing on a paper‐based analytical device for the fractionation of proteins

Isoelectric focusing plays a critical role in the analysis of complex protein samples. Conventionally, isoelectric focusing is implemented with carrier ampholytes in capillary or immobilized pH gradient gel. In this study, we successfully exhibited a carrier ampholyte‐free isoelectric focusing on paper‐based analytical device. Proof of the concept was visually demonstrated with color model proteins. Experimental results showed that not only a pH gradient was well established along the open paper fluidic channel as confirmed by pH indicator strip, the pH gradient range could also be tuned by the catholyte or anolyte. Furthermore, the isoelectric focusing fractions from the paper channel can be directly cut and recovered into solutions for post analysis with sodium dodecyl sulfate‐polyacrylamide gel electrophoresis and matrix‐assisted laser desorption/ionization‐time‐of‐flight mass spectrometry. This paper‐based isoelectric focusing method is fast, cheap, simple and easy to operate, and could potentially be used as a cost‐effective protein sample clean‐up method for target protein analysis with mass spectrometry.