Volume 134, Issue 5, 07 February 2011

A coarsegrained molecular model, which consists of a spherical particle and an orientation vector, is proposed to simulate lipidmembrane on a large length scale. The solvent is implicitly represented by an effective attractive interaction between particles. A bilayer structure is formed by orientationdependent (tilt and bending) potentials. In this model, the membrane properties (bending rigidity, line tension of membrane edge, area compression modulus, lateral diffusion coefficient, and flipflop rate) can be varied over broad ranges. The stability of the bilayer membrane is investigated via dropletvesicle transition. The rupture of the bilayer and wormlike micelle formation can be induced by an increase in the spontaneous curvature of the monolayermembrane.
 COMMUNICATIONS


Communication: Driven Brownian transport in eccentric septate channels
View Description Hide DescriptionIn eccentric septate channels the pores connecting adjacent compartments are shifted offaxis, either periodically or randomly, so that straight trajectories parallel to the axis are not allowed. Driven transport of a Brownian particle in such a channel is characterized by a strong suppression of the current and its dispersion. For large driving forces, both quantities approach an asymptotic value, which can be analytically approximated in terms of the stationary distribution of the particle exit times out of a single channel compartment.
 Top

 ARTICLES

 Theoretical Methods and Algorithms

Stability, structural, and magnetic phase diagrams of ternary ferromagnetic 3dtransitionmetal clusters with five and six atoms
View Description Hide DescriptionWe report a theoretical investigation of freestanding ternary clusters with x + y + z = 5 and 6. Our study is performed within density functional theory as implemented in the GAUSSIAN 03 set of programs and with the BPW91/SDD level of theory. We analyze the geometries, chemical order, local and total magnetic moments, binding energies, excess energies, and second difference in the energy in the whole range of composition, from which structural, magnetic, and stabilityphase diagrams are predicted for these cluster sizes. We determine the optimal stoichiometries for these clusters as regards the maximum total magnetic moment and stability.

Efficient implementation of the Hellmann–Feynman theorem in a diffusion Monte Carlo calculation
View Description Hide DescriptionKinetic and potential energies of systems of He atoms in the solid phase are computed at T = 0. Results at two densities of the liquid phase are presented as well. Calculations are performed by the multiweight extension to the diffusionMonte Carlo method that allows the application of the Hellmann–Feynman theorem in a robust and efficient way. This is a general method that can be applied in other situations of interest as well.

Timeoptimal control of spin 1/2 particles in the presence of radiation damping and relaxation
View Description Hide DescriptionWe consider the timeoptimal control of an ensemble of uncoupled spin 1/2 particles in the presence of relaxation and radiation damping effects, whose dynamics is governed by nonlinear equations generalizing the standard linear Bloch equations. For a single spin, the optimal control strategy can be fully characterized analytically. However, in order to take into account the inhomogeneity of the static magnetic field, an ensemble of isochromats at different frequencies must be considered. For this case, numerically optimized pulse sequences are computed and the dynamics under the corresponding optimal field is experimentally demonstrated using nuclear magnetic resonance techniques.

Selfassembly in block polyelectrolytes
View Description Hide DescriptionThe selfconsistent field theory (SCFT) complemented with the Poisson–Boltzmann equation is employed to explore selfassembly of polyelectrolytecopolymers composed of charged blocks A and neutral blocks B. We have extended SCFT to dissociating triblock copolymers and demonstrated our approach on three characteristic examples: (1) diblock copolymer (AB) melt, (2) symmetric triblock copolymer (ABA) melt, (3) triblock copolymer (ABA) solution with added electrolyte. For copolymer melts, we varied the composition (that is, the total fraction of Asegments in the system) and the charge density on A blocks and calculated the phase diagram that contains ordered mesophases of lamellar, gyroid, hexagonal, and bcc symmetries, as well as the uniform disordered phase. The phase diagram of charged block copolymer melts in the charge density – system composition coordinates is similar to the classical phase diagram of neutral block copolymer melts, where the composition and the Flory mismatch interaction parameter are used as variables. We found that the transitions between the polyelectrolyte mesophases with the increase of charge density occur in the same sequence, from lamellar to gyroid to hexagonal to bcc to disordered morphologies, as the mesophase transitions for neutral diblocks with the decrease of . In a certain range of compositions, the phase diagram for charged triblock copolymers exhibits unexpected features, allowing for transitions from hexagonal to gyroid to lamellar mesophases as the charge density increases. Triblock polyelectrolyte solutions were studied by varying the charge density and solvent concentration at a fixed copolymer composition. Transitions from lamellar to gyroid and gyroid to hexagonal morphologies were observed at lower polymer concentrations than the respective transitions in the similar neutral copolymer, indicating a substantial influence of the charge density on phase behavior.

Stiffness detection and reduction in discrete stochastic simulation of biochemical systems
View Description Hide DescriptionTypical multiscale biochemical models contain fastscale and slowscale reactions, where “fast” reactions fire much more frequently than “slow” ones. This feature often causes stiffness in discrete stochastic simulation methods such as Gillespie's algorithm and the TauLeaping method leading to inefficient simulation. This paper proposes a new strategy to automatically detect stiffness and identify species that cause stiffness for the TauLeaping method, as well as two stiffness reduction methods. Numerical results on a stiff decaying dimerization model and a heat shock protein regulation model demonstrate the efficiency and accuracy of the proposed methods for multiscale biochemical systems.

Sources of the deficiencies in the popular SPC/E and TIP3P models of water
View Description Hide DescriptionMotivated by the results of Vega et al. [J. Phys. Condens. Matter20, 153101 (2008)] about the phase diagram of water, and by the results of Kiss and Baranyai [J. Chem. Phys.131, 204310 (2009)] about the properties of gasphase clusters, we carried out a comparative study of the structure modeled by SPC/E and TIP3P interactions in ambient liquidwater. The gasphase clusters of SPC/E and TIP3P models show erroneous structures, while TIP4Ptype models, either polarizable or not, provide qualitatively correct results. The trimers of SPC/E and TIP3P are planar in gas phase, contrary to experimental and TIP4Ptype models. The aim of this study was to see whether traces of these false geometriescharacteristic to SPC/E and TIP3P in gas phase can also be found in the liquid phase. For this purpose we selected trimers formed by adjacent neighbors of water molecules in the liquid and calculated their geometrical features. We determined angles formed by the HO bonds of the molecules with OO vectors and with the normal vector of the OOO plane in the selected trimers. Our results showed that, despite high temperature, the SPC/E and TIP3P water contains larger number of planar arrangements than other TIP4Ptype models. Although structural differences presented in this study are small, they are accurately detectable. These results weaken the reliability of studies obtained by the SPC/E or TIP3P models even in the liquid phase.

Escorted free energy simulations
View Description Hide DescriptionWe describe a strategy to improve the efficiency of free energy estimates by reducing dissipation in nonequilibrium Monte Carlo simulations. This strategy generalizes the targeted free energyperturbation approach [C. Jarzynski, Phys. Rev. E 65, 046122 (2002)] to nonequilibrium switching simulations, and involves generating artificial, “escorted” trajectories by coupling the evolution of the system to updates in external work parameter. Our central results are: (1) a generalized fluctuation theorem for the escorted trajectories, and (2) estimators for the free energy difference ΔF in terms of these trajectories. We illustrate the method and its effectiveness on model systems.

A multiple replica approach to simulate reactive trajectories
View Description Hide DescriptionA method to generate reactive trajectories, namely equilibrium trajectories leaving a metastable state and ending in another one is proposed. The algorithm is based on simulating in parallel many copies of the system, and selecting the replicas which have reached the highest values along a chosen onedimensional reaction coordinate. This reaction coordinate does not need to precisely describe all the metastabilities of the system for the method to give reliable results. An extension of the algorithm to compute transition times from one metastable state to another one is also presented. We demonstrate the interest of the method on two simple cases: A onedimensional twowell potential and a twodimensional potential exhibiting two channels to pass from one metastable state to another one.

Ringpolymer instanton method for calculating tunneling splittings
View Description Hide DescriptionThe semiclassical instanton expression for the tunneling splitting between two symmetric wells is rederived, starting from the ringpolymer representation of the quantum partition function. This leads to simpler mathematics by replacing functional determinants with matrix determinants. By exploiting the simple Hückellike structure of the matrices, we derive an expression for the instanton tunneling splitting in terms of a minimum on the potential surface of a linearpolymer. The latter is a section cut out of a ringpolymer, consisting of an infinite number of beads, which describes a periodic orbit on the inverted potential surface. The approach is straightforward to generalize to multiple dimensions, and we demonstrate that it is computationally practical by carrying out instanton calculations of tunneling splittings in and malonaldehyde in full dimensionality.

Cartesian coupled coherent states simulations: Ne_{ n }Br_{2} dissociation as a test case
View Description Hide DescriptionIn this article, we describe coupled coherent states (CCS) simulations of vibrational predissociation of weakly bounded complexes. The CCS method is implemented in the Cartesian frame in a manner that is similar to classical molecular dynamics. The calculated lifetimes of the vibrationally excited NeBr_{2}(ν) complexes agree with experiment and previous calculations. Although the CCS method is, in principle, a fully quantum approach, in practice it typically becomes a semiclassical technique at long times. This is especially true following dissociation events. Consequently, it is very difficult to converge the quantum calculations of the final Br_{2} vibrational distributions after predissociation and of the autocorrelation functions. However, the main advantage of the method is that it can be applied with relative ease to determine the lifetimes of larger complexes and, in order to demonstrate this, preliminary results for tetra and pentaatomic clusters are reported.

Implementation of the analytic energy gradient for the combined timedependent density functional theory/effective fragment potential method: Application to excitedstate molecular dynamics simulations
View Description Hide DescriptionExcitedstate quantum mechanics/molecular mechanics molecular dynamics simulations are performed, to examine the solvent effects on the fluorescence spectra of aqueous formaldehyde. For that purpose, the analytical energy gradient has been derived and implemented for the linearresponse timedependent density functional theory (TDDFT) combined with the effective fragment potential (EFP) method. The EFP method is an efficient ab initio based polarizable model that describes the explicit solvent effects on electronic excitations, in the present work within a hybrid TDDFT/EFP scheme. The new method is applied to the excitedstate MD of aqueous formaldehyde in the nπ* state. The calculated π*→n transition energy and solvatochromic shift are in good agreement with other theoretical results.

Maximum likelihoodbased analysis of singlemolecule photon arrival trajectories
View Description Hide DescriptionIn this work we explore the statistical properties of the maximum likelihoodbased analysis of onecolor photon arrival trajectories. This approach does not involve binning and, therefore, all of the information contained in an observed photon strajectory is used. We study the accuracy and precision of parameter estimates and the efficiency of the Akaike information criterion and the Bayesian information criterion (BIC) in selecting the true kinetic model. We focus on the low excitation regime where photon trajectories can be modeled as realizations of Markov modulated Poisson processes. The number of observed photons is the key parameter in determining model selection and parameter estimation. For example, the BIC can select the true threestate model from competing two, three, and fourstate kinetic models even for relatively short trajectories made up of 2 × 10^{3}photons. When the intensity levels are wellseparated and 10^{4}photons are observed, the twostate model parameters can be estimated with about 10% precision and those for a threestate model with about 20% precision.

Explicitly timedependent coupled cluster singles doubles calculations of laserdriven manyelectron dynamics
View Description Hide DescriptionWe report explicitly timedependent coupled cluster singles doubles (TDCCSD) calculations, which simulate the laserdriven correlated manyelectron dynamics in molecular systems. Small molecules, i.e., HF, H_{2}O, NH_{3}, and CH_{4}, are treated mostly with polarized valence double zeta basis sets. We determine the coupled clusterground states by imaginary time propagation for these molecules. Excited stateenergies are obtained from the Fourier transform of the timedependent dipole moment after an ultrashort, broadband laser excitation. The timedependent expectation values are calculated from the complex cluster amplitudes using the corresponding configuration interaction singles doubles wave functions. Also resonant laser excitations of these excited states are simulated, in order to explore the limits for the numerical stability of our current TDCCSD implementation, which uses timeindependent molecular orbitals to form excited configurations.

A waterswap reaction coordinate for the calculation of absolute protein–ligand binding free energies
View Description Hide DescriptionThe accurate prediction of absolute protein–ligand binding free energies is one of the grand challenge problems of computational science. Binding free energy measures the strength of binding between a ligand and a protein, and an algorithm that would allow its accurate prediction would be a powerful tool for rational drug design. Here we present the development of a new method that allows for the absolute binding free energy of a protein–ligand complex to be calculated from first principles, using a single simulation. Our method involves the use of a novel reaction coordinate that swaps a ligand bound to a protein with an equivalent volume of bulk water. This waterswap reaction coordinate is built using an identity constraint, which identifies a cluster of water molecules from bulk water that occupies the same volume as the ligand in the protein active site. A dual topology algorithm is then used to swap the ligand from the active site with the identified water cluster from bulk water. The free energy is then calculated using replica exchange thermodynamic integration. This returns the free energy change of simultaneously transferring the ligand to bulk water, as an equivalent volume of bulk water is transferred back to the protein active site. This, directly, is the absolute binding free energy. It should be noted that while this reaction coordinate models the binding process directly, an accurate force field and sufficient sampling are still required to allow for the binding free energy to be predicted correctly. In this paper we present the details and development of this method, and demonstrate how the potential of mean force along the waterswap coordinate can be improved by calibrating the softcore Coulomb and LennardJones parameters used for the dual topology calculation. The optimal parameters were applied to calculations of protein–ligand binding free energies of a neuraminidase inhibitor (oseltamivir), with these results compared to experiment. These results demonstrate that the waterswap coordinate provides a viable and potentially powerful new route for the prediction of protein–ligand binding free energies.

Variational second order density matrix study of : Importance of subspace constraints for sizeconsistency
View Description Hide DescriptionVariational second order density matrix theory under “twopositivity” constraints tends to dissociate molecules into unphysical fractionally charged products with too low energies. We aim to construct a qualitatively correct potential energy surface for by applying subspaceenergy constraints on mono and diatomic subspaces of the molecular basis space. Monoatomic subspace constraints do not guarantee correct dissociation: the constraints are thus geometry dependent. Furthermore, the number of subspace constraints needed for correct dissociation does not grow linearly with the number of atoms. The subspace constraints do impose correct chemical properties in the dissociation limit and sizeconsistency, but the structure of the resulting second order density matrix method does not exactly correspond to a system of noninteracting units.

Efficient timedependent density functional theory approximations for hybrid density functionals: Analytical gradients and parallelization
View Description Hide DescriptionIn this paper, we present the implementation of efficient approximations to timedependent 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 (RIJ) 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 groundstate 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 AgTB2helicate (156 atoms, 2732 basis functions), it is demonstrated that the COSX algorithm parallelizes almost perfectly (speedup ∼26–29 for 30 processors). The exchangecorrelation terms also parallelize well (speedup ∼27–29 for 30 processors). The solution of the Zvector 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 Becke3LeeYangParr energy and gradient calculation on the AgTB2helicate 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 frozencore approximation. The algorithms described in this work are implemented into the ORCA electronic structure system.

A quantum propagator for pathintegral simulations of rigid molecules
View Description Hide DescriptionThe expression for the quantum propagator for rigid tops, proposed by Müser and Berne [Phys. Rev. Lett. 77, 2638 (1996)], has been extended to asymmetric tops. Pathintegral Monte Carlo simulations are provided that show that the quantum propagator proposed in this work exactly reproduces the rotational energy of free asymmetric tops as evaluated from the partition function. This propagator can subsequently be used in pathintegral simulations of condensed phases if a rigid molecular model is used.

Tensor decomposition in postHartree–Fock methods. I. Twoelectron integrals and MP2
View Description Hide DescriptionA new approximation for postHartree–Fock (HF) methods is presented applying tensor decomposition techniques in the canonical product tensor format. In this ansatz, multidimensional tensors like integrals or wavefunction parameters are processed as an expansion in onedimensional representing vectors. This approach has the potential to decrease the computational effort and the storage requirements of conventional algorithms drastically while allowing for rigorous truncation and error estimation. For postHF ab initio methods, for example, storage is reduced to with d being the number of dimensions of the full tensor,R being the expansion length (rank) of the tensor decomposition, and n being the number of entries in each dimension (i.e., the orbital index). If all tensors are expressed in the canonical format, the computational effort for any subsequent tensor contraction can be reduced to . We discuss details of the implementation, especially the decomposition of the twoelectron integrals, the AO–MO transformation, the Møller–Plesset perturbation theory (MP2) energy expression and the perspective for coupled cluster methods. An algorithm for rank reduction is presented that parallelizes trivially. For a set of representative examples, the scaling of the decomposition rank with system and basis set size is found to be for the AO integrals, for the MO integrals, and for the MP2 t _{2}amplitudes (N denotes a measure of system size) if the upper bound of the error in the ℓ^{2}norm is chosen as ε = 10^{−2}. This leads to an error in the MP2 energy in the order of mHartree.

Vibrational coupled cluster response theory: A general implementation
View Description Hide DescriptionThe calculation of vibrational contributions to molecular properties using vibrational coupled cluster (VCC) response theory is discussed. General expressions are given for expectation values, linear response functions, and transition moments. It is shown how these expressions can be evaluated for arbitrary levels of excitation in the wave function parameterization as well as for arbitrary coupling levels in the potential and property surfaces. The convergence of the method is assessed by benchmark calculations on formaldehyde. Furthermore, excitation energies and infrared intensities are calculated for the fundamental vibrations of furan using VCC limited to up to twomode and up to threemode excitations, VCC[2] and VCC[3], as well as VCC with full twomode and approximate threemode couplings, VCC[2pt3].