Search

In this paper, we present the implementation of efficient approximations to time-dependent density functional theory (TDDFT) within the Tamm-Dancoff approximation (TDA) for hybrid density functionals. For the calculation of the TDDFT/TDA excitation energies and analytical gradients, we combine the resolution of identity (RI-J) algorithm for the computation of the Coulomb terms and the recently introduced "chain of spheres exchange" (COSX) algorithm for the calculation of the exchange terms. It is shown that for extended basis sets, the RIJCOSX approximation leads to speedups of up to 2 orders of magnitude compared to traditional methods, as demonstrated for hydrocarbon chains. The accuracy of the adiabatic transition energies, excited state structures, and vibrational frequencies is assessed on a set of 27 excited states for 25 molecules with the configuration interaction singles and hybrid TDDFT/TDA methods using various basis sets. Compared to the canonical values, the typical error in transition energies is of the order of 0.01 eV. Similar to the ground-state results, excited state equilibrium geometries differ by less than 0.3 pm in the bond distances and 0.5° in the bond angles from the canonical values. The typical error in the calculated excited state normal coordinate displacements is of the order of 0.01, and relative error in the calculated excited state vibrational frequencies is less than 1%. The errors introduced by the RIJCOSX approximation are, thus, insignificant compared to the errors related to the approximate nature of the TDDFT methods and basis set truncation. For TDDFT/TDA energy and gradient calculations on Ag-TB2-helicate (156 atoms, 2732 basis functions), it is demonstrated that the COSX algorithm parallelizes almost perfectly (speedup ~26-29 for 30 processors). The exchange-correlation terms also parallelize well (speedup ~27-29 for 30 processors). The solution of the Z-vector equations shows a speedup of ~24 on 30 processors. The parallelization efficiency for the Coulomb terms can be somewhat smaller (speedup ~15-25 for 30 processors), but their contribution to the total calculation time is small. Thus, the parallel program completes a Becke3-Lee-Yang-Parr energy and gradient calculation on the Ag-TB2-helicate in less than 4 h on 30 processors. We also present the necessary extension of the Lagrangian formalism, which enables the calculation of the TDDFT excited state properties in the frozen-core approximation. The algorithms described in this work are implemented into the ORCA electronic structure system.

The weighted stochastic simulation algorithm (wSSA) was developed by Kuwahara and Mura [J. Chem. Phys. 129, 165101 (2008)] to efficiently estimate the probabilities of rare events in discrete stochastic systems. The wSSA uses importance sampling to enhance the statistical accuracy in the estimation of the probability of the rare event. The original algorithm biases the reaction selection step with a fixed importance sampling parameter. In this paper, we introduce a novel method where the biasing parameter is state-dependent. The new method features improved accuracy, efficiency, and robustness.

The dynamics of many biological processes of interest, such as the folding of a protein, are slow and complicated enough that a single molecular dynamics simulation trajectory of the entire process is difficult to obtain in any reasonable amount of time. Moreover, one such simulation may not be sufficient to develop an understanding of the mechanism of the process, and multiple simulations may be necessary. One approach to circumvent this computational barrier is the use of Markov state models. These models are useful because they can be constructed using data from a large number of shorter simulations instead of a single long simulation. This paper presents a new Bayesian method for the construction of Markov models from simulation data. A Markov model is specified by (τ,P,T), where τ is the mesoscopic time step, P is a partition of configuration space into mesostates, and T is an N(P)×N(P) transition rate matrix for transitions between the mesostates in one mesoscopic time step, where N(P) is the number of mesostates in P. The method presented here is different from previous Bayesian methods in several ways. (1) The method uses Bayesian analysis to determine the partition as well as the transition probabilities. (2) The method allows the construction of a Markov model for any chosen mesoscopic time-scale τ. (3) It constructs Markov models for which the diagonal elements of T are all equal to or greater than 0.5. Such a model will be called a "consistent mesoscopic Markov model" (CMMM). Such models have important advantages for providing an understanding of the dynamics on a mesoscopic time-scale. The Bayesian method uses simulation data to find a posterior probability distribution for (P,T) for any chosen τ. This distribution can be regarded as the Bayesian probability that the kinetics observed in the atomistic simulation data on the mesoscopic time-scale τ was generated by the CMMM specified by (P,T). An optimization algorithm is used to find the most probable CMMM for the chosen mesoscopic time step. We applied this method of Markov model construction to several toy systems (random walks in one and two dimensions) as well as the dynamics of alanine dipeptide in water. The resulting Markov state models were indeed successful in capturing the dynamics of our test systems on a variety of mesoscopic time-scales.

Gillespie's direct method (DM) [D. Gillespie, J. Chem. Phys. 81, 2340 (1977)] for exact stochastic simulation of chemical reaction systems has been widely adopted. It is easy to implement but requires large computation for relatively large systems. Recently, two more efficient methods, next reaction method (NRM) [M. A. Gibson and J. Bruck, J. Phys. Chem. A 105, 1876 (2000)] and optimized DM (ODM) [Y. Cao et al., J. Chem. Phys. 121, 4059 (2004)], have been developed to improve simulation speed. It has been demonstrated that the ODM is the state-of-the-art most efficient method for exact stochastic simulation of most practical reaction systems. In this paper, we first develop an exact stochastic simulation algorithm named ODMK that is more efficient than the ODM. We then develop an approximate method named K-skip method to further accelerate simulation. Using two chemical reaction systems, we demonstrate that our ODMK and K-skip method can save 20%-30% and 70%-80% simulation time, respectively, comparing to the ODM. We also show that our ODMK and K-skip method provide almost the same simulation accuracy as the ODM.

In this paper a Clifford algebra-based method is applied to calculate polymer chain conformations. The approach enables the calculation of the position of an atom in space with the knowledge of the bond length (l), valence angle (theta), and rotation angle (phi) of each of the preceding bonds in the chain. Hence, the set of geometrical parameters {l(i),theta(i),phi(i)} yields all the position coordinates p(i) of the main chain atoms. Moreover, the method allows the calculation of side chain conformations and the computation of rotations of chain segments. With these features it is, in principle, possible to generate conformations of any type of chemical structure. This method is proposed as an alternative for the classical approach by matrix algebra. It is more straightforward and its final symbolic representation considerably simpler than that of matrix algebra. Approaches for realistic modeling by means of incorporation of energetic considerations can be combined with it. This article, however, is entirely focused at showing the suitable mathematical framework on which further developments and applications can be built.

The molecular dynamics expression of heat flux, originally derived by Irving and Kirkwood [J. Chem. Phys. 18, 817 (1950)] for pairwise potentials, is generalized in this paper for systems with many-body potentials. The original formula consists of a kinetic part and a potential part, and the latter term is found in the present study to be expressible as a summation of contributions from all the many-body potentials defined in the system. The energy transfer among a set of sites for which a many-body potential is defined is discussed and evaluated by the rate of increase in the kinetic energy of each site due to the potential, and its accumulation over all the potentials in the system is shown to make up the potential part of the generalized expression. A molecular dynamics simulation for liquid n-octane was performed to demonstrate the applicability of the new expression obtained in this study to measure the heat flux and to elucidate the contributions of inter- and intramolecular potentials to heat conduction.

The evaluation of interactions between nearby particles constitutes the majority of the computational workload involved in classical molecular dynamics (MD) simulations. In this paper, we introduce a new method for the parallelization of range-limited particle interactions that proves particularly suitable to MD applications. Because it applies not only to pairwise interactions but also to interactions involving three or more particles, the method can be used for evaluation of both nonbonded and bonded forces in a MD simulation. It requires less interprocessor data transfer than traditional spatial decomposition methods at all but the lowest levels of parallelism. It gains an additional practical advantage in certain commonly used interprocessor communication networks by distributing the communication burden more evenly across network links and by decreasing the associated latency. When used to parallelize MD, it further reduces communication requirements by allowing the computations associated with short-range nonbonded interactions, long-range electrostatics, bonded interactions, and particle migration to use much of the same communicated data. We also introduce certain variants of this method that can significantly improve the balance of computational load across processors.

In this paper, we present simulations of the decay of quantum coherence between vibrational states of I(2) in its ground (X) electronic state embedded in a cryogenic Kr matrix. We employ a numerical method based on the semiclassical limit of the quantum Liouville equation, which allows the simulation of the evolution and decay of quantum vibrational coherence using classical trajectories and ensemble averaging. The vibrational level-dependent interaction of the I(2)(X) oscillator with the rare-gas environment is modeled using a recently developed method for constructing state-dependent many-body potentials for quantum vibrations in a many-body classical environment [J. M. Riga, E. Fredj, and C. C. Martens, J. Chem. Phys. 122, 174107 (2005)]. The vibrational dephasing rates gamma(0n) for coherences prepared between the ground vibrational state mid R:0 and excited vibrational state mid R:n are calculated as a function of n and lattice temperature T. Excellent agreement with recent experiments performed by Karavitis et al. [Phys. Chem. Chem. Phys. 7, 791 (2005)] is obtained.

In this paper we present the results of a stochastic model of reversible biochemical reaction networks that are being driven through an open boundary, such that the system is interacting with its surrounding environment with explicit material exchange. The stochastic model is based on the master equation approach and is intimately related to the grand canonical ensemble of statistical mechanics. We show that it is possible to analytically calculate the joint probability function of the random variables describing the number of molecules in each state of the system for general linear networks. Definitions of reaction chemical potentials and conductances follow from inherent properties of this model, providing a description of energy dissipation in the system. We are also able to suggest novel methods for experimentally determining reaction fluxes and biochemical affinities at nonequilibrium steady state as well as the overall network connectivity.

The finite-temperature string method proposed by E, et al. [W. E, W. Ren, and E. Vanden-Eijnden, Phys. Rev. B 66, 052301 (2002)] is a very effective way of identifying transition mechanisms and transition rates between metastable states in systems with complex energy landscapes. In this paper, we discuss the theoretical background and algorithmic details of the finite-temperature string method, as well as the application to the study of isomerization reaction of the alanine dipeptide, both in vacuum and in explicit solvent. We demonstrate that the method allows us to identify directly the isocommittor surfaces, which are approximated by hyperplanes, in the region of configuration space where the most probable transition trajectories are concentrated. These results are verified subsequently by computing directly the committor distribution on the hyperplanes that define the transition state region.

We describe the dynamic behavior of a 1-stearoyl-2-oleoyl-phosphatidylethanolamine (SOPE) bilayer from a 20 ns molecular dynamics simulation. The dynamics of individual molecules are characterized in terms of (2)H spin-lattice relaxation rates, nuclear overhauser enhancement spectroscopy (NOESY) cross-relaxation rates, and lateral diffusion coefficients. Additionally, we describe the dynamics of hydrogen bonding through an analysis of hydrogen bond lifetimes and the time evolution of clusters of hydrogen bonded lipids. The simulated trajectory is shown to be consistent with experimental measures of internal, intermolecular, and diffusive motion. Consistent with our analysis of SOPE structure in the companion paper, we see hydrogen bonding dominating the dynamics of the interface region. Comparison of (2)H T(1) relaxation rates for chain methylene segments in phosphatidylcholine and phosphatidylethanolamine bilayers indicates that slower motion resulting from hydrogen bonding extends at least three carbons into the hydrophobic core. NOESY cross-relaxation rates compare well with experimental values, indicating the observed hydrogen bonding dynamics are realistic. Calculated lateral diffusion rates (4 +/ -1 x 10(-8) cm(2)s) are comparable, though somewhat lower than, those determined by pulsed field gradient NMR methods.

The dynamical solution of a well-mixed, nonlinear stochastic chemical kinetic system, described by the Master equation, may be exactly computed using the stochastic simulation algorithm. However, because the computational cost scales with the number of reaction occurrences, systems with one or more "fast" reactions become costly to simulate. This paper describes a hybrid stochastic method that partitions the system into subsets of fast and slow reactions, approximates the fast reactions as a continuous Markov process, using a chemical Langevin equation, and accurately describes the slow dynamics using the integral form of the "Next Reaction" variant of the stochastic simulation algorithm. The key innovation of this method is its mechanism of efficiently monitoring the occurrences of slow, discrete events while simultaneously simulating the dynamics of a continuous, stochastic or deterministic process. In addition, by introducing an approximation in which multiple slow reactions may occur within a time step of the numerical integration of the chemical Langevin equation, the hybrid stochastic method performs much faster with only a marginal decrease in accuracy. Multiple examples, including a biological pulse generator and a large-scale system benchmark, are simulated using the exact and proposed hybrid methods as well as, for comparison, a previous hybrid stochastic method. Probability distributions of the solutions are compared and the weak errors of the first two moments are computed. In general, these hybrid methods may be applied to the simulation of the dynamics of a system described by stochastic differential, ordinary differential, and Master equations.

Recently, we have proposed an efficient scheme for Monte Carlo simulations, the multiple "time step" Monte Carlo (MTS-MC) [J. Chem. Phys. 117, 8203 (2002)] based on the separation of the potential interactions into two additive parts. In this paper, the structural and thermodynamic properties of the simple point charge water model combined with the Ewald sum are compared for the MTS-MC real-/reciprocal-space split of the Ewald summation and the common Metropolis Monte Carlo method. We report a number of observables as a function of CPU time calculated using MC and MTS-MC. The correlation functions indicate that speedups on the order of 4.5-7.5 can be obtained for systems of 108-500 waters for n=10 splitting parameter., ((c) 2004 American Institute of Physics.)

The accurate simulation of biologically active macromolecules faces serious limitations that originate in the treatment of electrostatics in the empirical force fields. The current use of "partial charges" is a significant source of errors, since these vary widely with different conformations. By contrast, the molecular electrostatic potential (MEP) obtained through the use of a distributed multipole moment description, has been shown to converge to the quantum MEP outside the van der Waals surface, when higher order multipoles are used. However, in spite of the considerable improvement to the representation of the electronic cloud, higher order multipoles are not part of current classical biomolecular force fields due to the excessive computational cost. In this paper we present an efficient formalism for the treatment of higher order multipoles in Cartesian tensor formalism. The Ewald "direct sum" is evaluated through a McMurchie-Davidson formalism [L. McMurchie and E. Davidson, J. Comput. Phys. 26, 218 (1978)]. The "reciprocal sum" has been implemented in three different ways: using an Ewald scheme, a particle mesh Ewald (PME) method, and a multigrid-based approach. We find that even though the use of the McMurchie-Davidson formalism considerably reduces the cost of the calculation with respect to the standard matrix implementation of multipole interactions, the calculation in direct space remains expensive. When most of the calculation is moved to reciprocal space via the PME method, the cost of a calculation where all multipolar interactions (up to hexadecapole-hexadecapole) are included is only about 8.5 times more expensive than a regular AMBER 7 [D. A. Pearlman et al., Comput. Phys. Commun. 91, 1 (1995)] implementation with only charge-charge interactions. The multigrid implementation is slower but shows very promising results for parallelization. It provides a natural way to interface with continuous, Gaussian-based electrostatics in the future. It is hoped that this new formalism will facilitate the systematic implementation of higher order multipoles in classical biomolecular force fields., ((c) 2004 American Institute of Physics)

In this paper, we describe a new and fully coherent stochastic surface hopping method for simulating mixed quantum-classical systems. We illustrate the approach on the simple but unforgiving problem of quantum evolution of a two-state quantum system in the limit of unperturbed pure state dynamics and for dissipative evolution in the presence of both stationary and nonstationary random environments. We formulate our approach in the Liouville representation and describe the density matrix elements by ensembles of trajectories. Population dynamics are represented by stochastic surface hops for trajectories representing diagonal density matrix elements. These are combined with an unconventional coherent stochastic hopping algorithm for trajectories representing off-diagonal quantum coherences. The latter generalizes the binary (0,1) "probability" of a trajectory to be associated with a given state to allow integers that can be negative or greater than unity in magnitude. Unlike existing surface hopping methods, the dynamics of the ensembles are fully entangled, correctly capturing the coherent and nonlocal structure of quantum mechanics. [ABSTRACT FROM AUTHOR]

Stochastic simulation for in silico studies of large biochemical networks requires a great amount of computational time. We recently proposed a new exact simulation algorithm, called the rejectionbased stochastic simulation algorithm (RSSA) [Thanh et al., J. Chem. Phys. 141(13), 134116 (2014)], to improve simulation performance by postponing and collapsing as much as possible the propensity updates. In this paper, we analyze the performance of this algorithm in detail, and improve it for simulating large-scale biochemical reaction networks. We also present a new algorithm, called simultaneous RSSA (SRSSA), which generates many independent trajectories simultaneously for the analysis of the biochemical behavior. SRSSA improves simulation performance by utilizing a single data structure across simulations to select reaction firings and forming trajectories. The memory requirement for building and storing the data structure is thus independent of the number of trajectories. The updating of the data structure when needed is performed collectively in a single operation across the simulations. The trajectories generated by SRSSA are exact and independent of each other by exploiting the rejection-based mechanism. We test our new improvement on real biological systems with a wide range of reaction networks to demonstrate its applicability and efficiency. [ABSTRACT FROM AUTHOR]

The thermodynamic stability of a folded RNA is intricately tied to the counterions and the free energy of this interaction must be accounted for in any realistic RNA simulations. Extending a tight-binding model published previously, in this paper we investigate the fundamental structure of charges arising from the interaction between small functional RNA molecules and divalent ions such as Mg2+ that are especially conducive to stabilizing folded conformations. The characteristic nature of these charges is utilized to construct a discretely connected energy landscape that is then traversed via a novel application of a deterministic graph search technique. This search method can be incorporated into larger simulations of small RNA molecules and provides a fast and accurate way to calculate the free energy arising from the interactions between an RNA and divalent counterions. The utility of this algorithm is demonstrated within a fully atomistic Monte Carlo simulation of the P4-P6 domain of the Tetrahymena group I intron, in which it is shown that the counterionmediated free energy conclusively directs folding into a compact structure. [ABSTRACT FROM AUTHOR]

In this paper, a previous coarse-grain model [J. T. Padding and W. J. Briels, J. Chem. Phys. 117, 925 (2002)] to simulate melts of linear polymers has been adapted to simulate polymers with more complex hierarchies. Bond crossings between highly coarse-grained soft particles are prevented by applying an entanglement algorithm. We first test our method on a virtual branch point inside a linear chain to make sure it works effectively when linking two linear arms. Next, we apply our method to study the diffusive and rheological behaviors of a melt of three-armed stars. We find that the diffusive behavior of the three-armed star is very close to that of a linear polymer with the same molecular weight, while its rheological properties are close to those of a linear chain with molecular mass equal to that of the longest linear sub-chain in the star. [ABSTRACT FROM AUTHOR]

In this paper, we present time-reversible simulation algorithms for rigid bodies in the quaternion representation. By advancing a time-reversible algorithm [Y. Kajima, M. Hiyama, S. Ogata, and T. Tamura, J. Phys. Soc. Jpn. 80, 114002 (2011)] that requires iterations in calculating the angular velocity at each time step, we propose two kinds of iteration-free fast time-reversible algorithms. They are easily implemented in codes. The codes are compared with that of existing algorithms through demonstrative simulation of a nanometer-sized water droplet to find their stability of the total energy and computation speeds. [ABSTRACT FROM AUTHOR]

We propose a parameter-free algorithm for the identification of nearest neighbors. The algorithm is very easy to use and has a number of advantages over existing algorithms to identify nearest-neighbors. This solid-angle based nearest-neighbor algorithm (SANN) attributes to each possible neighbor a solid angle and determines the cutoff radius by the requirement that the sum of the solid angles is 4π. The algorithm can be used to analyze 3D images, both from experiments as well as theory, and as the algorithm has a low computational cost, it can also be used 'on the fly' in simulations. In this paper, we describe the SANN algorithm, discuss its properties, and compare it to both a fixed-distance cutoff algorithm and to a Voronoi construction by analyzing its behavior in bulk phases of systems of carbon atoms, Lennard-Jones particles and hard spheres as well as in Lennard-Jones systems with liquid-crystal and liquid-vapor interfaces. [ABSTRACT FROM AUTHOR]

Using computer simulations to model the folding of proteins into their native states is computationally expensive due to the extraordinarily low degeneracy of the ground state. In this paper, we develop an efficient way to sample these folded conformations using Wang Landau sampling coupled with the configurational bias method (which uses an unphysical 'temperature' that lies between the collapse and folding transition temperatures of the protein). This method speeds up the folding process by roughly an order of magnitude over existing algorithms for the sequences studied. We apply this method to study the adsorption of intrinsically disordered hydrophobic polar protein fragments on a hydrophobic surface. We find that these fragments, which are unstructured in the bulk, acquire secondary structure upon adsorption onto a strong hydrophobic surface. Apparently, the presence of a hydrophobic surface allows these random coil fragments to fold by providing hydrophobic contacts that were lost in protein fragmentation. [ABSTRACT FROM AUTHOR]

A controlled quantum system possesses a search landscape defined by the target physical objective as a function of the controls. This paper focuses on the landscape for the transition probability Pi → f between the states of a finite level quantum system. Traditionally, the controls are applied fields; here, we extend the notion of control to also include the Hamiltonian structure, in the form of time independent matrix elements. Level sets of controls that produce the same transition probability value are shown to exist at the bottom Pi → f = 0.0 and top Pi → f = 1.0 of the landscape with the field and/or Hamiltonian structure as controls. We present an algorithm to continuously explore these level sets starting from an initial point residing at either extreme value of Pi → f. The technique can also identify control solutions that exhibit the desirable properties of (a) robustness at the top and (b) the ability to rapidly rise towards an optimal control from the bottom. Numerical simulations are presented to illustrate the varied control behavior at the top and bottom of the landscape for several simple model systems. [ABSTRACT FROM AUTHOR]

Density functional theory as proposed by Rosenfeld [Phys. Rev. Lett. 63, 980 (1989)] is used to study hard sphere mixture exposed by cylindrically symmetric external field. Exploiting the symmetry of the system, explicit formulas for the weighted densities are derived. The resulting density profiles are compared with new grand canonical Monte Carlo simulations. The comparison reveals very good agreement between the predicted and simulated results even at high densities and very narrow pores. Finally, simple algorithms for computing complete elliptic functions of the first and second kinds that occur in the derived formulae are presented to make the paper self-contained. [ABSTRACT FROM AUTHOR]

A novel theory of an equation of state based on excluded volume and formulated in two preceding papers for gases and gaseous mixtures is extended to the entire density range by considering higher (beginning from the third) approximations of the theory. The algorithm of constructing higher approximations is elaborated. Equations of state are deduced using the requirement of maximum simplicity and contain a single free parameter to be chosen by reason of convenience or simplicity or to be used as a fitting parameter with respect to the computer simulation database. In this way, precise equations of state are derived for the hard-sphere fluid in the entire density range. On the side, the theory reproduces most known earlier equations of state for hard spheres and determines their place in the hierarchy of approximations. Equations of state for van der Waals fluids are also presented, and their critical parameters are estimated. © 2004 American Institute of Physics. [ABSTRACT FROM AUTHOR]

In this paper, we present computer simulation results concerning interdiffusion of fully compatible components in symmetric binary (AB) polymer mixtures in solutions. The simulation is performed in two dimensions using the algorithm based on the dynamic lattice liquid model. The solvent molecules are taken into account explicitly. The evolution of the concentration profiles in time at an interface is studied for chain lengths N=2,4,8,16 for three polymer concentrations [lowercase_phi_synonym]=0.1,0.5,0.9. The tracer diffusion coefficients for polymer chains and for the solvent are obtained by monitoring the mean square displacements of their center of mass. The relationships between coefficients of interdiffusion and self-diffusion are tested. © 2004 American Institute of Physics. [ABSTRACT FROM AUTHOR]

A lattice computer simulation of a symmetric A–B–A triblock copolymer melt is reported. This melt is quenched, in simulation, from an athermal state to 39 different temperatures using cooperative motion algorithm. Energy, specific heat, copolymer end-to-end distance, bridging fraction, lamellar spacing, concentration profiles, and microstructure visualizations are reported. The quenching simulation results are compared with those obtained by alternative thermal treatments, that is by slow heating and slow cooling. Quenches yield data consistent with theory and experiment, whereas slow cooling and slow heating results do not capture the expected behavior for the lamellar spacing and the bridging fraction. Finally, at very low temperatures, below the conventional order–disorder transition temperature, an additional ordering is recorded, from a conventional lamellar phase to a lamellar structure showing copolymer junction points condensed into a two-dimensional plane. © 2003 American Institute of Physics. [ABSTRACT FROM AUTHOR]

The reversible reference system propagator algorithm (r-RESPA), based on a Trotter factorization of the classical propagator, is tested in a molecular dynamics simulation of solid C60. We show how, with an appropriate subdivision of the interaction potential and with a careful balancing of the integration parameters, one can adopt large time steps and impressively efficient r-RESPA integrators which yield the same dynamics obtained by means of the small time-step Verlet algorithm. The results presented here show that the use of r-RESPA integrators speeds up the simulation by a factor of between 20–40 with respect to the standard Verlet algorithms. [ABSTRACT FROM AUTHOR]

The Hartree–Fock problem provides the conceptual and mathematical underpinning of a large portion of quantum chemistry. As efforts in quantum technology aim to enhance computational chemistry algorithms, the Hartree–Fock method, central to many other numerical approaches, is a natural target for quantum enhanced algorithms. While quantum computers and quantum simulation offer many prospects for the future of modern chemistry, the non-deterministic polynomial-complete Hartree–Fock problem is not a likely candidate. We highlight this fact from a number of perspectives including computational complexity, practical examples, and the full characterization of energy landscapes for simple systems. [ABSTRACT FROM AUTHOR]

The Parrinello-Rahman algorithm for imposing a general state of stress in periodic molecular dynamics simulations is widely used in the literature and has been implemented in many readily available molecular dynamics codes. However, what is often overlooked is that this algorithm controls the second Piola-Kirchhoff stress as opposed to the true (Cauchy) stress. This can lead to misinterpretation of simulation results because (1) the true stress that is imposed during the simulation depends on the deformation of the periodic cell, (2) the true stress is potentially very different from the imposed second Piola-Kirchhoff stress, and (3) the true stress can vary significantly during the simulation even if the imposed second Piola-Kirchhoff is constant. We propose a simple modification to the algorithm that allows the true Cauchy stress to be controlled directly. We then demonstrate the efficacy of the new algorithm with the example of martensitic phase transformations under applied stress. [ABSTRACT FROM AUTHOR]

The weighted histogram analysis method (WHAM) including its binless extension has been developed independently in several different contexts, and widely used in chemistry, physics, and statistics, for computing free energies and expectations from multiple ensembles. However, this method, while statistically efficient, is computationally costly or even infeasible when a large number, hundreds or more, of distributions are studied.We develop a locally WHAM (local WHAM) from the perspective of simulations of simulations (SOS), using generalized serial tempering (GST) to resample simulated data from multiple ensembles. The local WHAM equations based on one jump attempt per GST cycle can be solved by optimization algorithms orders of magnitude faster than standard implementations of global WHAM, but yield similarly accurate estimates of free energies to global WHAM estimates. Moreover, we propose an adaptive SOS procedure for solving local WHAM equations stochastically when multiple jump attempts are performed per GST cycle. Such a stochastic procedure can lead to more accurate estimates of equilibrium distributions than local WHAM with one jump attempt per cycle. The proposed methods are broadly applicable when the original data to be "WHAMMED" are obtained properly by any sampling algorithm including serial tempering and parallel tempering (replica exchange). To illustrate the methods, we estimated absolute binding free energies and binding energy distributions using the binding energy distribution analysis method from one and two dimensional replica exchange molecular dynamics simulations for the beta-cyclodextrin-heptanoate host-guest system. In addition to the computational advantage of handling large datasets, our two dimensional WHAM analysis also demonstrates that accurate results similar to those from well-converged data can be obtained from simulations for which sampling is limited and not fully equilibrated. [ABSTRACT FROM AUTHOR]

To respect the nature of discrete parts in a system, stochastic simulation algorithms (SSAs) must update for each action (i) all part counts and (ii) each action's probability of occurring next and its timing. This makes it expensive to simulate biological networks with well-connected "hubs" such as ATP that affect many actions. Temperature and volume also affect many actions and may be changed significantly in small steps by the network itself during fever and cell growth, respectively. Such trends matter for evolutionary questions, as cell volume determines doubling times and fever may affect survival, both key traits for biological evolution. Yet simulations often ignore such trends and assume constant environments to avoid many costly probability updates. Such computational con-venience precludes analyses of important aspects of evolution. Here we present "Lazy Updating," an add-on for SSAs designed to reduce the cost of simulating hubs. When a hub changes, Lazy Updating postpones all probability updates for reactions depending on this hub, until a threshold is crossed. Speedup is substantial if most computing time is spent on such updates. We implemented Lazy Updating for the Sorting Direct Method and it is easily integrated into other SSAs such as Gillespie's Direct Method or the Next Reaction Method. Testing on several toy models and a cel-lular metabolism model showed >10× faster simulations for its use-cases--with a small loss of accuracy. Thus we see Lazy Updating as a valuable tool for some special but important simulation problems that are difficult to address efficiently otherwise. [ABSTRACT FROM AUTHOR]

Calculation of the effective diffusion coefficient D(x), depending on the longitudinal coordinate x in 2D channels with periodically corrugated walls, is revisited. Instead of scaling the transverse lengths and applying the standard homogenization techniques, we propose an algorithm based on formulation of the problem in the complex plane. A simple model is solved to explain the behavior of D(x) in the channels with short periods L, observed by Brownian simulations of Dagdug et al. [J. Chem. Phys. 133, 034707 (2010)]. [ABSTRACT FROM AUTHOR]

Numerical simulation of the time evolution of a spatially homogeneous chemical system is always of great interest. Gillespie first developed the exact stochastic simulation algorithm (SSA), which is accurate but time-consuming. Recently, many approximate schemes of the SSA are proposed to speed up simulation. Presented here is the N-leap method, which guarantees the validity of the leap condition and at the same time keeps the efficiency. In many cases, N-leap has better performance than the widely-used τ-leap method. The details of the N-leap method are described and several examples are presented to show its validity. [ABSTRACT FROM AUTHOR]

We demonstrate theoretically that it may be possible to encode states of a multi-qubit system into the progression of quantized motional/vibrational levels of an ion trapped in a weakly anharmonic potential. Control over such register of quantum information is achieved by applying oscillatory radio-frequency fields shaped optimally for excitation of the desired state-to-state transitions. Anharmonicity of the vibrational spectrum plays a key role in this approach to the control and quantum computation, since it allows resolving different state-to-state transitions and addressing them selectively. Optimal control theory is used to derive pulses for implementing the four-qubit version of Shor's algorithm in a single step. Accuracy of the qubit-state transformations, reached in the numerical simulations, is around 0.999. Very detailed insight is obtained by analysis of the time-evolution of state populations and by spectral analysis of the optimized pulse. [ABSTRACT FROM AUTHOR]

Within the framework of transition path sampling (TPS), activation energies can be computed as path ensemble averages without a priori information about the reaction mechanism [C. Dellago and P. G. Bolhuis, Mol. Simul. 30, 795 (2004)]. Activation energies computed for different conditions can then be used to determine by numerical integration the rate constant for a system of interest from the rate constant known for a reference system. However, in systems with complex potential energy surfaces, multiple reaction pathways may exist making ergodic sampling of trajectory space difficult. Here, we present a combination of TPS with the Wang-Landau (WL) flat-histogram algorithm for an efficient sampling of the transition path ensemble. This method, denoted by WL-TPS, has the advantage that from one single simulation, activation energies at different temperatures can be determined even for systems with multiple reaction mechanisms. The proposed methodology for rate constant calculations does not require the knowledge of the reaction coordinate and is generally applicable to Arrhenius and non-Arrhenius processes. We illustrate the applicability of this technique by studying a two-dimensional toy system consisting of a triatomic molecule immersed in a fluid of repulsive soft disks. We also provide an expression for the calculation of activation volumes from path averages such that the pressure dependence of the rate constant can be obtained by numerical integration. [ABSTRACT FROM AUTHOR]

Parallel tempering (PT) is a set of techniques for accelerating thermal-equilibrium sampling in systems where the exploration of configuration space is hindered by energy barriers. With standard PT algorithms, the computational effort scales unfavorably with system size, so that it is difficult to apply them to large systems. We propose local PT algorithms, for which the computational effort is proportional to the number of degrees of freedom. We demonstrate the effectiveness of the new algorithms on two one-dimensional model systems, showing that results for selected observables are correctly reproduced, and that practical linear scaling is achieved. We show also that the algorithms are readily applied to systems in higher dimensions. We note the prospects for studying large extended systems, including surfaces and interfaces. [ABSTRACT FROM AUTHOR]

We present the partial-propensity stochastic simulation algorithm with composition-rejection sampling (PSSA-CR). It is an exact formulation of the stochastic simulation algorithm (SSA) for well-stirred systems of coupled chemical reactions. The new formulation is a partial-propensity variant [R. Ramaswamy, N. González-Segredo, and I. F. Sbalzarini, J. Chem. Phys. 130, 244104 (2009)] of the composition- rejection SSA [A. Slepoy, A. P. Thompson, and S. J. Plimpton, J. Chem. Phys. 128, 205101 (2008)]. The computational cost of this new formulation is bounded by a constant for weakly coupled reaction networks, and it increases at most linearly with the number of chemical species for strongly coupled reaction networks. PSSA-CR thus combines the advantages of partial-propensity methods and the composition-rejection SSA, providing favorable scaling of the computational cost for all classes of reaction networks. [ABSTRACT FROM AUTHOR]

Distributed replica sampling [T. Rodinger et al., J. Chem. Theory Comput. 2, 725 (2006)] is a simple and general scheme for Boltzmann sampling of conformational space by computer simulation in which multiple replicas of the system undergo a random walk in reaction coordinate or temperature space. Individual replicas are linked through a generalized Hamiltonian containing an extra potential energy term or bias which depends on the distribution of all replicas, thus enforcing the desired sampling distribution along the coordinate or parameter of interest regardless of free energy barriers. In contrast to replica exchange methods, efficient implementation of the algorithm does not require synchronicity of the individual simulations. The algorithm is inherently suited for large-scale simulations using shared or heterogeneous computing platforms such as a distributed network. In this work, we build on our original algorithm by introducing Boltzmann-weighted jumping, which allows moves of a larger magnitude and thus enhances sampling efficiency along the reaction coordinate. The approach is demonstrated using a realistic and biologically relevant application; we calculate the standard binding free energy of benzene to the L99A mutant of T4 lysozyme. Distributed replica sampling is used in conjunction with thermodynamic integration to compute the potential of mean force for extracting the ligand from protein and solvent along a nonphysical spatial coordinate. Dynamic treatment of the reaction coordinate leads to faster statistical convergence of the potential of mean force than a conventional static coordinate, which suffers from slow transitions on a rugged potential energy surface. [ABSTRACT FROM AUTHOR]

We investigated the generation dependent shape and internal structure of star-burst dendrimers under good solvent conditions using small angle x-ray scattering and molecular modeling. Measurements have been performed on poly(amidoamine) dendrimers with generations ranging from g=0 up to g=8 at low concentrations in methanol. We described the static form factor P(q) by a model taking into account the compact, globular shape as well as the loose, polymeric character of dendrimers. Monomer distributions within dendrimers are of special interest for potential applications and have been characterized by the pair correlation function γ(r), as well as by the monomer and end-group density profile, ρ(r) and ρe(r), respectively. Monomer density profiles and γ(r) can be derived from P(q) by modeling and via a model independent approach using the inverse Fourier transformation algorithm first introduced by Glatter. Experimental results are compared with computer simulations performed for single dendrimers of various generations using the cooperative motion algorithm. The simulation gives direct access to γ(r) and ρ(r), allows an independent determination of P(q), and yields in addition to the scattering experiment information about the distribution of the end groups. Excellent qualitative agreement between experiment and simulation has been found. © 2004 American Institute of Physics. [ABSTRACT FROM AUTHOR]

The shape of isolated two-dimensional linear and star polymers, with and without excluded volume, have been determined by Monte Carlo simulations using a pivot algorithm. Stars with 3, 4, 5, and 6 arms are studied. The high efficiency of the pivot algorithm has allowed us to obtain accurate data for long chains. Good agreement is found with existing analytical expressions for nonexcluded volume chains and with other computer simulations for excluded volume chains. [ABSTRACT FROM AUTHOR]

The stochastic simulation algorithm (SSA) has been widely used for simulating biochemical reaction networks. SSA is able to capture the inherently intrinsic noise of the biological system, which is due to the discreteness of species population and to the randomness of their reciprocal interactions. However, SSA does not consider other sources of heterogeneity in biochemical reaction systems, which are referred to as extrinsic noise. Here, we extend two simulation approaches, namely, the integration-based method and the rejection-based method, to take extrinsic noise into account by allowing the reaction propensities to vary in time and state dependent manner. For both methods, new efficient implementations are introduced and their efficiency and applicability to biological models are investigated. Our numerical results suggest that the rejection-based method performs better than the integration-based method when the extrinsic noise is considered. [ABSTRACT FROM AUTHOR]

Numerical approaches can provide useful information about the microscopic processes underlying photocurrent generation in organic solar cells (OSCs). Among them, the Kinetic Monte Carlo (KMC) method is conceptually the simplest, but computationally the most intensive. A less demanding alternative is potentially represented by so-called Master Equation (ME) approaches, where the equations describing particle dynamics rely on the mean-field approximation and their solution is attained numerically, rather than stochastically. The description of charge separation dynamics, the treatment of electrostatic interactions and numerical stability are some of the key issues which have prevented the application of these methods to OSC modelling, despite of their successes in the study of charge transport in disordered system. Here we describe a three-dimensional ME approach to photocurrent generation in OSCs which attempts to deal with these issues. The reliability of the proposed method is tested against reference KMC simulations on bilayer heterojunction solar cells. Comparison of the current-voltage curves shows that the model well approximates the exact result for most devices. The largest deviations in current densities are mainly due to the adoption of the mean-field approximation for electrostatic interactions. The presence of deep traps, in devices characterized by strong energy disorder, may also affect result quality. Comparison of the simulation times reveals that the ME algorithm runs, on the average, one order of magnitude faster than KMC. [ABSTRACT FROM AUTHOR]

Explicit integrators for real-time propagation of time-dependent Kohn-Sham equations are compared regarding their suitability for performing large-scale simulations. Four algorithms are implemented and assessed for both stability and accuracy within a plane-wave pseudopotential framework, employing the adiabatic approximation to the exchange-correlation functional. Simulation results for a single sodium atom and a sodium atom embedded in bulk magnesium oxide are discussed. While the first-order Euler scheme and the second-order finite-difference scheme are unstable, the fourth-order Runge-Kutta scheme is found to be conditionally stable and accurate within this framework. Excellent parallel scalability of the algorithm up to more than a thousand processors is demonstrated for a system containing hundreds of electrons, evidencing the suitability for large-scale simulations based on real-time propagation of time-dependent Kohn-Sham equations. [ABSTRACT FROM AUTHOR]

The opportunities offered by utilizing time-independent Hamiltonian structure as controls are explored for manipulating quantum dynamics. Two scenarios are investigated using different manifestations of Hamiltonian structure to illustrate the generality of the concept. In scenario I, optimally shaped electrostatic potentials are generated to flexibly control electron scattering in a two-dimensional subsurface plane of a semiconductor. A simulation is performed showing the utility of optimally setting the individual voltages applied to a multi-pixel surface gate array in order to produce a spatially inhomogeneous potential within the subsurface scattering plane. The coherent constructive and destructive electron wave interferences are manipulated by optimally adjusting the potential shapes to alter the scattering patterns. In scenario II, molecular vibrational wave packets are controlled by means of optimally selecting the Hamiltonian structure in cooperation with an applied field. As an illustration of the concept, a collection (i.e., a level set) of dipole functions is identified where each member serves with the same applied electric field to produce the desired final transition probability. The level set algorithm additionally found Hamiltonian structure controls exhibiting desirable physical properties. The prospects of utilizing the applied field and Hamiltonian structure simultaneously as controls is also explored. The control scenarios I and II indicate the gains offered by algorithmically guided molecular or material discovery for manipulating quantum dynamics phenomenon. [ABSTRACT FROM AUTHOR]

We develop and test a new elementary Monte Carlo move for use in the off-lattice simulation of polymer systems. This novel Parallel-Rotation algorithm (ParRot) permits moving very efficiently torsion angles that are deeply inside long chains in melts. The parallel-rotation move is extremely simple and is also demonstrated to be computationally efficient and appropriate for Monte Carlo simulation. The ParRot move does not affect the orientation of those parts of the chain outside the moving unit. The move consists of a concerted rotation around four adjacent skeletal bonds. No assumption is made concerning the backbone geometry other than that bond lengths and bond angles are held constant during the elementary move. Properly weighted sampling techniques are needed for ensuring detailed balance because the new move involves a correlated change in four degrees of freedom along the chain backbone. The ParRot move is supplemented with the classical Metropolis Monte Carlo, the Continuum-Configurational-Bias, and Reptation techniques in an isothermal-isobaric Monte Carlo simulation of melts of short and long chains. Comparisons are made with the capabilities of other Monte Carlo techniques to move the torsion angles in the middle of the chains. We demonstrate that ParRot constitutes a highly promising Monte Carlo move for the treatment of long polymer chains in the off-lattice simulation of realistic models of dense polymer systems. © 2001 American Institute of Physics. [ABSTRACT FROM AUTHOR]