Lagrangian formulation with dissipation of Born-Oppenheimer molecular dynamics using the density-functional tight-binding method

An important element determining the time requirements of Born-Oppenheimer molecular dynamics (BOMD) is the convergence rate of the self-consistent solution of Roothaan equations (SCF). We show here that improved convergence and dynamics stability can be achieved by use of a Lagrangian formalism of BOMD with dissipation (DXL-BOMD). In the DXL-BOMD algorithm, an auxiliary electronic variable (e.g., the electron density or Fock matrix) is propagated and a dissipative force is added in the propagation to maintain the stability of the dynamics. Implementation of the approach in the self-consistent charge density functional tight-binding method makes possible simulations that are several hundred picoseconds in lengths, in contrast to earlier DFT-based BOMD calculations, which have been limited to tens of picoseconds or less. The increase in the simulation time results in a more meaningful evaluation of the DXL-BOMD method. A comparison is made of the number of iterations (and time) required for convergence of the SCF with DXL-BOMD and a standard method (starting with a zero charge guess for all atoms at each step), which gives accurate propagation with reasonable SCF convergence criteria. From tests using NVE simulations of C(2)F(4) and 20 neutral amino acid molecules in the gas phase, it is found that DXL-BOMD can improve SCF convergence by up to a factor of two over the standard method. Corresponding results are obtained in simulations of 32 water molecules in a periodic box. Linear response theory is used to analyze the relationship between the energy drift and the correlation of geometry propagation errors.     

A note on density matrix extrapolation and multiple solutions of the unrestricted Hartree–Fock equations

Unrestricted Hartree–Fock (UHF) SCF–MO calculations on the doublet reaction surface for the addition of methylidyne (CH) to ethylene (C₂H₄) using the standard extrapolation techniques of the GAUSSIAN 70 program show erratic behavior. On the other hand, the potential energy surface calculated without extrapolation of the density matrix and by using the final density matrix of a neighboring point as the initial guess for the density matrix for the new point gave a smooth potential curve without any kinks or erratic pattern. Therefore, the density extrapolation technique should be used with particular caution in UHF calculations.     

Accelerated convergence in SCF calculations and the level shifting technique

The possibility of improving the convergence rate in SCF calculations by exploiting the essential arbitrariness of the diagonal matrix elements of the Fock operator formed in the MO basis is studied. It is shown that it is possible to accelerate convergence in many cases by adopting a different form of level shifting technique (called the reverse level shifting technique by us).     

Geometric integration in Born-Oppenheimer molecular dynamics

Geometric integration schemes for extended Lagrangian self-consistent Born-Oppenheimer molecular dynamics, including a weak dissipation to remove numerical noise, are developed and analyzed. The extended Lagrangian framework enables the geometric integration of both the nuclear and electronic degrees of freedom. This provides highly efficient simulations that are stable and energy conserving even under incomplete and approximate self-consistent field (SCF) convergence. We investigate three different geometric integration schemes: (1) regular time reversible Verlet, (2) second order optimal symplectic, and (3) third order optimal symplectic. We look at energy conservation, accuracy, and stability as a function of dissipation, integration time step, and SCF convergence. We find that the inclusion of dissipation in the symplectic integration methods gives an efficient damping of numerical noise or perturbations that otherwise may accumulate from finite arithmetics in a perfect reversible dynamics.     

An unconventional scf method for calculations on large molecules

An unconventional SCF method for calculations on large molecules with more than 100 basis functions is described. Storage problems which arise in conventional SCF schemes when storing more than 10⁷ integrals are avoided by repeated calculation of integrals. The resulting increase in computational times is kept at a reasonable level by (a) improving the initial guess, (b) accelerating convergence, (c) employing a recursive construction of the Fock matrix, and (d) eliminating insignificant integrals from the calculation by a density-weighted cutoff criterion. Sample calculations show that, compared with conventional SCF calculations, computational times increase by 25%–75% depending on the basis set and the shape of the molecule.     

An efficient self-consistent field method for large systems of weakly interacting components

An efficient method for removing the self-consistent field (SCF) diagonalization bottleneck is proposed for systems of weakly interacting components. The method is based on the equations of the locally projected SCF for molecular interactions (SCF MI) which utilize absolutely localized nonorthogonal molecular orbitals expanded in local subsets of the atomic basis set. A generalization of direct inversion in the iterative subspace for nonorthogonal molecular orbitals is formulated to increase the rate of convergence of the SCF MI equations. Single Roothaan step perturbative corrections are developed to improve the accuracy of the SCF MI energies. The resulting energies closely reproduce the conventional SCF energy. Extensive test calculations are performed on water clusters up to several hundred molecules. Compared to conventional SCF, speedups of the order of (N/O)2 have been achieved for the diagonalization step, where N is the size of the atomic orbital basis, and O is the number of occupied molecular orbitals.     

Large scale Hartree-Fock calculations with conventional SCF algorithm. Influence of integral and index compression on Fock matrix construction

It is shown that a compression of two-electron integrals and their indices significantly improves efficiency of the conventional self-consistent field (SCF) algorithm for a solution of the Hartree-Fock equation by decrease the Fock matrix calculation time. The improvement is reached not only due to a reduction of the integral file size, but mainly because data compression reduces or even can eliminate a cache conflict in data transfer from the hard drive to the main computer memory. Thus, the conventional SCF algorithm with the data compression becomes very efficient and permits to carry out large-scale Hartree-Fock calculations. The largest Hartree-Fock calculations have been performed for RNA 433D structure from the PDB data bank with 6080 basis functions formed from 6-31G basis on a workstation with 1 GHz Alpha processor.     

Exact-exchange Hartree–Fock calculations for periodic systems. I. Illustration of the method

A scheme is presented for performing linear-combination-of-atomic-orbitals (LCAO ) self-consistent-field (SCF ) ab initio Hartree–Fock calculations of the electronic structure of periodic systems. The main aspects which characterize the present method are (i) a thorough discussion of both translational and local symmetry properties and the derivation of general formulas for the transformation of all the relevant monoelectronic and bielectronic terms under symmetry operators. (ii) The use of general yet powerful criteria for the truncation of infinite sums; in particular, the Coulomb electron–electron interactions are subdivided into terms corresponding to intersecting or nonintersecting charge distributions; the latter are grouped into shell contributions and the interaction is evaluated by multipolar expansions; the exchange interaction may be evaluated with great precision by retaining a relatively small number of two-electron integrals according to a truncation criterion which fully preserves its nonlocal character. (iii) The use of a procedure for performing integrals over k , as needed in the evaluation of the Fermi energy and in the reconstruction of the Fock matrix, which is particularly simple because it employs partially intersecting small spheres as integration subdomains where linear extrapolation is admitted. A comparison is finally made of our fundamental equations in the critical SCF stage with those obtainable by a recent proposal which uses Fourier transforms to express Coulomb and exchange integrals.     

Gaussian overlap approximation in the projected Hartree–Fock method

A method for the approximate calculation of matrix elements with respect to projected Hartree–Fock wave functions is proposed. The method is tested on some calculations in the many-parameter AMO method. It is found that the approximation reduces the amount of work, involved in the evaluation of the energy, by a factor of five and that it reproduces the exact values to within a few per cent.     

Implementation of divide-and-conquer method including Hartree-Fock exchange interaction

The divide-and-conquer (DC) method, which is one of the linear-scaling methods avoiding explicit diagonalization of the Fock matrix, has been applied mainly to pure density functional theory (DFT) or semiempirical molecular orbital calculations so far. The present study applies the DC method to such calculations including the Hartree-Fock (HF) exchange terms as the HF and hybrid HF/DFT. Reliability of the DC-HF and DC-hybrid HF/DFT is found to be strongly dependent on the cut-off radius, which defines the localization region in the DC formalism. This dependence on the cut-off radius is assessed from various points of view: that is, total energy, energy components, local energies, and density of states. Additionally, to accelerate the self-consistent field convergence in DC calculations, a new convergence technique is proposed.     

N2-time-dependent SCF scheme

A new procedure for the Fock matrix operator construction is proposed. Its application for RHF calculations on diatomic molecules using Slater orbital basis sets shows that the computation time for the new SCF procedure is proportional to the square of the basis set size.     

LFMO和NBO的定域性和作用能解析

为分析苯并分子C12H6的垂直共振能(VRE), 建立了定域片断分子轨道(LFMO)和自然键轨道(NBO)两种基组, 并在两种基组之上进行NBO能量分析和Morokuma能量分解. 在NBO能量分析中, 两种基组的VRE都是稳定的; 而在Morokuma能量分解中, VRE的稳定性取决于基组. 在NBO能量分析中, Fock矩阵的一次性对角化忽视了σ体系和π体系之间的电子耦合作用. 故NBO基组和NBO能量分析方法在计算VRE时似乎都不合理.     

The trust-region self-consistent field method: towards a black-box optimization in Hartree-Fock and Kohn-Sham theories

The trust-region self-consistent field (TRSCF) method is presented for optimizing the total energy E(SCF) of Hartree-Fock theory and Kohn-Sham density-functional theory. In the TRSCF method, both the Fock/Kohn-Sham matrix diagonalization step to obtain a new density matrix and the step to determine the optimal density matrix in the subspace of the density matrices of the preceding diagonalization steps have been improved. The improvements follow from the recognition that local models to E(SCF) may be introduced by carrying out a Taylor expansion of the energy about the current density matrix. At the point of expansion, the local models have the same gradient as E(SCF) but only an approximate Hessian. The local models are therefore valid only in a restricted region-the trust region-and steps can only be taken with confidence within this region. By restricting the steps of the TRSCF model to be inside the trust region, a monotonic and significant reduction of the total energy is ensured in each iteration of the TRSCF method. Examples are given where the TRSCF method converges monotonically and smoothly, but where the standard DIIS method diverges.     

Structure, stability and vibrational spectrum of the hydrogen-bonded complex between HNO3 and H2O. Ab initio and DFT studies

The structure, stability and vibrational spectrum of the binary complex between HONO2 and H2O have been investigated using ab initio calculations at SCF and MP2 levels with different basis sets and B3LYP/6-31G(d,p) calculations. Full geometry optimization was made for the complex studied. It was established that the hydrogen-bonded H2O...HONO2 complex has a planar structure. The corrected values of the dissociation energy at the SCF and MP2 levels and B3LYP calculations are indicative of relatively strong OH...O hydrogen-bonded interaction. The changes in the vibrational characteristics (vibrational frequencies and infrared intensities) arising from the hydrogen bonding between HONO2 and H2O have been estimated by using the ab initio calculations at SCF and MP2 levels and B3LYP/6-31G(d,p) calculations. It was established that the most sensitive to the complexation is the stretching O-H vibration from HONO2. In agreement with the experiment, its vibrational frequency in the complex is shifted to lower wavenumbers. The predicted frequency shift with the B3LYP/6-31G(d,p) calculations (-439 cm(-1)) is in the best agreement with the experimentally measured (-498 cm(-1)). The intensity of this vibration increases dramatically upon hydrogen bonding. The ab initio calculations at the SCF level predict an increase up to five times; at the MP2 level up to 10 times and the B3LYP/6-31G(d,p) predicted increase is up to 17 times. The good agreement between the predicted values of the frequency shifts and those experimentally observed show that the structure of the hydrogen-bonded complex H2O...HONO2 is reliable.     

Efficient first-principles electronic dynamics

An efficient first-principles electronic dynamics method is introduced in this article. The approach we put forth relies on incrementally constructing a time-dependent Fock∕Kohn-Sham matrix using active space density screening method that reduces the cost of computing two-electron repulsion integrals. An adaptive stepsize control algorithm is developed to optimize the efficiency of the electronic dynamics while maintaining good energy conservation. A selected set of model dipolar push-pull chromophore molecules are tested and compared with the conventional method of direct formation of the Fock∕Kohn-Sham matrix. While both methods considered herein take on identical dynamical simulation pathways for the molecules tested, the active space density screening algorithm becomes much more computationally efficient. The adaptive stepsize control algorithm, when used in conjunction with the dynamically active space method, yields a factor of ～3 speed-up in computational cost as observed in electronic dynamics using the time dependent density functional theory. The total computational cost scales nearly linear with increasing size of the molecular system.     

SCF theory of intermolecular interactions without basis set superposition error

Special SCF LCAO MO type equations are derived, permitting “supermolecule” calculations for intermolecular interactions, excluding basis set superposition error (BSSE) from the beginning on the basis of the “chemical Hamiltonian approach”. (No additional “monomer” calculations are necessary to correct for BSSE.) The formalism excluding the BSSE results in a non-Hermitean Fock matrix; an algorithm is proposed to obtain the required molecular orbitals, in which no integral transformation is needed.     

Accurate potential energy surface and quantum reaction rate calculations for the H+CH4-->H2+CH3 reaction

Calculations for the cumulative reaction probability N(E) (for J=0) and the thermal rate constant k(T) of the H+CH(4)-->H(2)+CH(3) reaction are presented. Accurate electronic structure calculations and a converged Shepard-interpolation approach are used to construct a potential energy surface which is specifically designed to allow the precise calculation of k(T) and N(E). Accurate quantum dynamics calculations employing flux correlation functions and multiconfigurational time-dependent Hartree wave packet propagation compute N(E) and k(T) based on this potential energy surface. The present work describes in detail the various convergence test performed to investigate the accuracy of the calculations at each step. These tests demonstrate the predictive power of the present calculations. In addition, approximate approaches for reaction rate calculations are discussed. A quite accurate approximation can be obtained from a potential energy surface which includes only interpolation points on the minimum energy path.     

Development of a restricted Hartree—Fock program INDMOL on PARAM: A highly parallel computer

Parallelization of the SCF method for closed-shell molecules on the highly parallel transputer-based system PARAM is described. The parallelization has been implemented on three different hardware and software environments: (1) a network of bare 64 transputers; (2) configuration 1 plus a back-end file system (BFS); and (3) configuration 2 with one INTEL i860 processor. The evaluation of electron repulsion integrals (ERIs) and setting up of the Fock matrix is carried out in parallel on 64 nodes using minimal communication strategies. A good load balance is achieved for ERI evaluation with the help of bounds, local symmetry features, and the shell concept, as well as a data randomization technique, resulting into almost linear speedup (for ERI evaluation). In configurations 2 and 3, BFS is used for parallel storage and retrieval of ERIs. Further, in 3 matrix operations are implemented as remote procedure calls on the i860 processor. Routine techniques of level shifting and extrapolation are used for accelerating SCF convergence. The resulting package, INDMOL, has been tested for some randomly selected molecules having up to 255 contractions. Using configuration 3, a factor of 2 to 5 in computation time is obtained over 1, for the systems for which the ERIs cannot be stored in the distributed core memory. In summary, a heterogeneous system, as in configuration 3, can indeed be optimally exploited for programming peculiar diverse requirements of the SCF procedure. © 1993 John Wiley & Sons, Inc.     

Is large-scale ab initio Hartree-Fock calculation chemically accurate? Toward improved calculation of biological molecule properties

Numerical errors in total energy values in large-scale Hartree–Fock calculations are discussed. To obtain total energy values within chemical accuracy, 0.01 kcal/mol, stricter numerical accuracy is required as basis size increases. In molecules with 10,000 basis sizes, such as proteins, numerical accuracy for total energy values must be retained to at least 11 digits (i.e., to the order of 1.0D-10) to keep accumulation of numerical errors less than the chemical accuracy (0.01 kcal/mol). With this criterion, we examined the sensitivity analysis of numerical accuracy in Hartree–Fock calculation by uniformly replacing the last bit of the mantissa part of a double-precision real number by zero in the Fock matrix construction step, the total energy calculation step, and the Fock matrix diagonalization step. Using a partial summation technique in the Fock matrix generation step, the numerical error for total energy value of molecules with basis size greater than 10,000 was within chemical accuracy (0.01 kcal/mol), whereas with the conventional method the numerical error with several thousand basis sets was larger than chemical accuracy. Computation of one Fock matrix element with parallel machines can include the partial summation technique automatically, so that parallel calculation yields not only high-performance computing but also more precise numerical solutions than the conventional sequential algorithm. We also found that the numerical error of the Householder-QR diagonalization routine is equal to or less than chemical accuracy, even with a matrix size of 10,000. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 443–454, 1999     

Half-projected Hartree–Fock natural orbitals for defining CAS–SCF active spaces

We have extended the range of systems to which the half-projected Hartree–Fock (HPHF ) method has been applied, including the triplet state of the wave function. In our implementation, DIIS overcomes the convergence difficulties reported in earlier studies. HPHF allows generation of a symmetry-broken wave function in regions of the potential energy surface where the RHF wave function is triplet-stable. The fractionally occupied natural orbitals (FONOS ) of the HPHF wave function are good starting vectors for CAS –SCF calculations. A CAS –SCF in the space defined by the HPHF FONOS should be used instead of the unrestricted natural orbital CAS –SCF method in regions of triplet stability and for small active space problems. We draw extensive comparisons between the results of both the UNO –CAS and HPNO –CAS methods and those of full CAS –SCF calculations. © 1993 John Wiley & Sons, Inc.