Car-Parrinello treatment for an approximate density-functional theory method

The authors formulate a Car-Parrinello treatment for the density-functional-based tight-binding method with and without self-consistent charge corrections. This method avoids the numerical solution of the secular equations, the principal drawback for large systems if the linear combination of atomic orbital ansatz is used. The formalism is applicable to finite systems and for supercells using periodic boundary conditions within the Gamma-point approximation. They show that the methodology allows the application of modern computational techniques such as sparse matrix storage and massive parallelization in a straightforward way. All present bottlenecks concerning computer time and consumption of memory and memory bandwidth can be removed. They illustrate the performance of the method by direct comparison with Born-Oppenheimer molecular dynamics calculations. Water molecules, benzene, the C(60) fullerene, and liquid water have been selected as benchmark systems.     

An efficient and robust integration scheme for the equations of motion of the multiconfiguration time-dependent Hartree (MCTDH) method

An efficient and robust integration scheme tailored to the equations of motion of the multiconfiguration time-dependent Hartree (MCTDH) method is presented. An error estimation allows the automatical adjustment of the step size and hence controls the integration error. The integration scheme decouples the MCTDH equations of motion into several disjoined subsystems, of which one determines the time evolution of the MCTDH-coefficients. While the conventional MCTDH equations are non-linear, the working equation for the MCTDH-coefficients becomes linear in the present integration scheme. To investigate the integrator’s performance it is applied to the photodissociation process of methyl iodide. The results of the novel integration scheme are in perfect agreement to those obtained by solving the MCTDH working equations conventionally. The computation time, however, is reduced by a factor of about ten when the new integration scheme is used to propagate large systems.     

A program system for computer integration of multistep reaction rate equations using the gear integration method

A program system is described for the integration of the rate equations resulting from large systems of elementary reactions. The Gear integration method is used for this problem, which frequently may exhibit stiffness instability when other integration methds are employed. No usage of the quasi-steady-state approximation is necessary. Ease in varying the reaction mechanism and simplicity of input structure are coupled with efficient execution and minimal demands for program storage as key features. The input–output structure, method of operation, and implementation are summarized, along with core storage requirements and execution times for trials using an IBM 360/44 computer.     

Integration of chemical stiff ODEs using exponential propagation method

In this paper, we study the numerical long time integration of large stiff systems of differential equations arising from chemical reactions by exponential propagation methods. These methods, which typically converge faster, use matrix-vector products with the exponential or other related function of the Jacobian that can be effectively approximated by Krylov sabspace methods. We equip these methods to an automatic stepsize control technique and apply the method of order 4 for numerical integration of some famous stiff chemical problems such as Belousov-Zhabotinskii reaction, the Chapman atmosphere, Hydrogen chemistry, chemical Akzo-Nobel problem and air pollution problem.     

Vectorization of the general Monte Carlo classical trajectory program VENUS

The general chemical dynamics computer program VENUS is used to perform classical trajectory simulations for large polyatomic systems, with many atoms and complicated potential energy functions. To simulate an ensemble of many trajectories requires a large amount of CPU time. Since each trajectory is independent, it is possible to parallel process a large set of trajectories instead of processing the trajectories by the conventional sequential approach. This enhances the vectorizability of the VENUS program, since the integration of Hamilton's equations of motion and the gradient evaluation, which comprise 97.8% of the CPU, can each be parallel processed. In this article, the vectorization and ensuing optimization of VENUS on the CRAY-YMP and IBM-3090 are presented in terms of both global strategies and technical details. A switching algorithm is designed to enhance the vector performance and to minimize the memory storage. A performance of 140 MFLOPS and a vector/scalar execution rate ratio of 10.6 are observed when this new version of VENUS is used to study the association of CH₃ with the H(Ar)₁₂ cluster on the CRAY-YMP.     

A contemporary linear representation theory for ordinary differential equations: probabilistic evolutions and related approximants for unidimensional autonomous systems

In this paper, we build on our previous research on probabilistic foundations of dynamical systems and introduce a theory of linear representation for ordinary differential equations. The theory is developed for explicit ODEs and can be further extended to cover implicit cases. In this report, we investigate the case of a canonical single unknown autonomous system. First we construct a linear representation to get an infinite linear ODE set with a constant coefficient matrix which can be transformed into an upper triangular form. Then we find its approximate truncated solutions. We describe a number of properties of the theory using this framework. The companion of this paper expands this canonical approach to cover multidimensional cases using the theory of folded arrays which is another line of research established by our research group.     

Application of Coulson's Contour integral method to linear finite polymers

Coulson's contour integration method has been applied to calculate the charge and bond order matrices and total energies of conjugated systems composed of fragments whose Hückel secular equations have been solved. Integral formulas have been derived for the calculation of these quantities for two coupled conjugated systems and for linear polymers. A rapid method for the numerical evaluation of the integral formulas is presented and applied to a linear finite polyene.     

Evaluation of atomic pressure in the multiple time‐step integration algorithm

In molecular dynamics (MD) calculations, reduction in calculation time per MD loop is essential. A multiple time‐step (MTS) integration algorithm, the RESPA (Tuckerman and Berne, J. Chem. Phys. 1992, 97, 1990–2001), enables reductions in calculation time by decreasing the frequency of time‐consuming long‐range interaction calculations. However, the RESPA MTS algorithm involves uncertainties in evaluating the atomic interaction‐based pressure (i.e., atomic pressure) of systems with and without holonomic constraints. It is not clear which intermediate forces and constraint forces in the MTS integration procedure should be used to calculate the atomic pressure. In this article, we propose a series of equations to evaluate the atomic pressure in the RESPA MTS integration procedure on the basis of its equivalence to the Velocity‐Verlet integration procedure with a single time step (STS). The equations guarantee time‐reversibility even for the system with holonomic constrants. Furthermore, we generalize the equations to both (i) arbitrary number of inner time steps and (ii) arbitrary number of force components (RESPA levels). The atomic pressure calculated by our equations with the MTS integration shows excellent agreement with the reference value with the STS, whereas pressures calculated using the conventional ad hoc equations deviated from it. Our equations can be extended straightforwardly to the MTS integration algorithm for the isothermal NVT and isothermal–isobaric NPT ensembles. © 2017 Wiley Periodicals, Inc.     

A comparative study of molecular dynamics in Cartesian and in internal coordinates: dynamical instability in the latter caused by nonlinearity of the equations of motion

The stability of a general molecular dynamics (MD) integration scheme is examined for simulations in generalized (internal plus external) coordinates (GCs). An analytic expression is derived for the local error in energy during each integration time step. This shows that the explicit dependence of the mass-matrix on GCs, which makes the system's Lagrange equations of motion nonlinear, causes MD simulations in GCs to be less stable than those in Cartesian coordinates (CCs). In terms of CCs, the corresponding mass-matrix depends only on atomic masses and thus atomistic motion is subject to the linear Newton equations, which makes the system more stable. Also investigated are two MD methods in GCs that utilize nonzero elements of the vibrational spectroscopic B-matrices. One updates positions and velocities in GCs that are iteratively adjusted so as to conform to the velocity Verlet equivalent in GCs. The other updates positions in GCs and velocities in CCs that are adjusted to satisfy the internal constraints of the new constrained WIGGLE MD scheme. The proposed methods are applied to an isolated n-octane molecule and their performances are compared with those of several CCMD schemes. The simulation results are found to be consistent with the analytic stability analysis. Finally, a method is presented for computing nonzero elements of B-matrices for external rotations without imposing the Casimir-Eckart conditions.     

Molecular dynamics on a distributed-memory multiprocessor

Dynamics simulations of molecular systems are notoriously computationally intensive. Using parallel computers for these simulations is important for reducing their turnaround time. In this article we describe a parallelization of the simulation program CHARMM for the Intel iPSC/860, a distributed memory multiprocessor. In the parallelization, the computational work is partitioned among the processors for core calculations including the calculation of forces, the integration of equations of motion, the correction of atomic coordinates by constraint, and the generation and update of data structures used to compute nonbonded interactions. Processors coordinate their activity using synchronous communication to exchange data values. Key data structures used are partitioned among the processors in nearly equal pieces, reducing the memory requirement per node and making it possible to simulate larger molecular systems. We examine the effectiveness of the parallelization in the context of a case study of a realistic molecular system. While effective speedup was achieved for many of the dynamics calculations, other calculations fared less well due to growing communication costs for exchanging data among processors. The strategies we used are applicable to parallelization of similar molecular mechanics and dynamics programs for distributed memory multiprocessors. © 1992 by John Wiley & Sons, Inc.     

Time-dependent exchange-correlation current density functionals with memory

Most present applications of time-dependent density functional theory use adiabatic functionals, i.e., the effective potential at time t is determined solely by the density at the same time. This paper discusses a method that aims to go beyond this approximation, by incorporating "memory" effects: the derived exchange-correlation potential will depend not only on present densities but also on the past. In order to ensure the potentials are causal, we formulate the action on the Keldysh contour for electrons in electromagnetic fields, from which we derive suitable Kohn-Sham equations. The exchange-correlation action is now a functional of the electron density and velocity field. A specific action functional is constructed which is Galilean invariant and yields a causal exchange-correlation vector potential for the Kohn-Sham equations incorporating memory effects. We show explicitly that the net exchange-correlation Lorentz force is zero. The potential is consistent with known dynamical properties of the homogeneous electron gas (in the linear response limit).     

Linear compartmental systems. I. kinetic analysis and derivation of their optimized symbolic equations

The study of many biological systems requires the application of a compartmental analysis, together with the use of isotopic tracers, parameter identification and methods to evaluate the mean parameters. For all this, the kinetic equations of the compartmental system as a function of its parameters are needed. In this paper, we present some considerations on the diagrams of connectivity of linear compartmental systems and obtain new properties from the matrix corresponding to the ordinary first-order linear differential equation systems which describe their kinetic behaviour. Using these properties, symbolic equations are obtained in a simplified form. These equations provide the instantaneous amount of substance in any compartment of the system when zero input is injected into one or more of the system compartments, solely as a function of those parameters of compartmental systems which really have an influence on the sought expression. This is unlike what happens in the other symbolic equations obtained in a previous contribution that included all the fractional transfer coefficients involved in the compartmental system, regardless of whether or not they had an influence on the instantaneous amount of substance.     

A Program package using stiff,sparse integration methods for the automatic solution of mass action kinetics equations

A portable program package, MACKSIM, for mass action chemical kinetics simulation, is discussed. As these kinetics are readily expressed in explicit mathematical terms, such a package contains two major and distinct modules, the numerical analysis and the user interface. For the first, MACKSIM uses the latest proven developments incorporating sparse matrix techniques in the backward difference predictor corrector methods originated by Gear for the integration of stiff ordinary differential equations, and thus requires minimal computing time to solve large systems of equations. For the second, the program provides a flexible interface which permits simple specification and variation of reactions, requires no special character input, and has no limit on the number of reactions or species involved other than that imposed by the size of the computer. The technology of these components is discussed briefly, the use of the package for standard reactions is illustrated, and current applications are mentioned.     

Analytical second derivatives of the energy in MNDO methods

Analytical second derivatives of the energy are derived and efficiently implemented for semiempirical MNDO-type methods including AM1, PM3, and MNDO/d. A new algorithm for the simultaneous solution of several CPHF equations is proposed in which the amount of memory required is independent of the number of iterations. The analytical approach is faster than the numerical approach typically by a factor of 5 and exhibits a reliable convergence over a wide range of molecules. The asymptotic memory and secondary storage requirements of the reported procedure can be as low as O(N²) without significant degradation of the performance. © 1996 by John Wiley & Sons, Inc.     

Effective solution of linear Diophantine equation systems with an application in chemistry

Systematic methods of the solution of linear Diophantine systems of equations and their application in decomposing complex chemical reactions are presented. The Contejean–Devie algorithm is improved. A new, linear programming based enumerative algorithm is described, which is applicable to large systems with large solutions. Mathematica implementations are tested and compared in important chemical examples.     

Convergence acceleration of iterative sequences. the case of scf iteration

Based on a recent method of Pople et al for the solution of large systems of linear equations, a procedure is given for accelerating the convergence of slowly converging quasi-Newton— Raphson type algorithms. This procedure is particularly advantageous if the number of parameters is so large that the calculation and storage of the hessian is no longer practical. Application to the SCF problem is treated in detail.     

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.     

Comparison of compuational methods for simulating nuclear Overhauser effects in NMR spectroscopy

Various algorithms for solving the Solomon equations describing nuclear Overhauser effects (nOes) in NMR spectroscopy have been compared. The applicability of the eigenvalue/eigenvector and the numerical integration approaches have been investigated. The eigenvalue/eigenvector approach is not a computationally efficient means of simulating nOe experiments in which a saturating radiofrequency field is applied during the time course. For experiments in which nOes develop in the absence of an RF field, this approach should only be used in simulating a full NOESY spectrum. Integration schemes have been found to be more efficient at simulating nOe experiments in which the nOe evolves in the presence of a saturating field, at simulating a partial set of initial perturbation experiments and at simulating a few rows or columns in a NOESY spectrum. Various integration schemes were applied to a two-spin system for which an analytic solution is available and to a model B-DNA oligonucleotide hexamer. The previously unused Taylor series algorithm was found to be superior to the Euler, midpoint, and fourth-order Runge–Kutta methods with regard to integration accuracy/computation time. An adaptive step size control routine for the Taylor series integration scheme was developed. Integration schemes can be speeded up in a simple fashion by introducing a distance cutoff for the dipolar interaction. Using a cutoff of 8 Å the Taylor series algorithm was able to compute the NOESY spectrum more rapidly than the eigenvalue/eigenvector algorithm for large spin systems at short mixing times. At longer mixing times the eigenvalue/eigenvector approach becomes the more efficient scheme.