Thermodynamics and equilibrium structure of Ne38 cluster: quantum mechanics versus classical

The equilibrium properties of classical Lennard-Jones (LJ38) versus quantum Ne38 Lennard-Jones clusters are investigated. The quantum simulations use both the path-integral Monte Carlo (PIMC) and the recently developed variational-Gaussian wave packet Monte Carlo (VGW-MC) methods. The PIMC and the classical MC simulations are implemented in the parallel tempering framework. The classical heat capacity Cv(T) curve agrees well with that of Neirotti et al. [J. Chem. Phys. 112, 10340 (2000)], although a much larger confining sphere is used in the present work. The classical Cv(T) shows a peak at about 6 K, interpreted as a solid-liquid transition, and a shoulder at approximately 4 K, attributed to a solid-solid transition involving structures from the global octahedral (Oh) minimum and the main icosahedral (C5v) minimum. The VGW method is used to locate and characterize the low energy states of Ne38, which are then further refined by PIMC calculations. Unlike the classical case, the ground state of Ne38 is a liquidlike structure. Among the several liquidlike states with energies below the two symmetric states (Oh and C5v), the lowest two exhibit strong delocalization over basins associated with at least two classical local minima. Because the symmetric structures do not play an essential role in the thermodynamics of Ne38, the quantum heat capacity is a featureless curve indicative of the absence of any structural transformations. Good agreement between the two methods, VGW and PIMC, is obtained. The present results are also consistent with the predictions by Calvo et al. [J. Chem. Phys. 114, 7312 (2001)] based on the quantum superposition method within the harmonic approximation. However, because of its approximate nature, the latter method leads to an incorrect assignment of the Ne38 ground state as well as to a significant underestimation of the heat capacity.     

Equilibrium density of states and thermodynamic properties of a model glass former

This paper presents an analysis of the thermodynamics of a model glass former. We have performed equilibrium sampling of a popular binary Lennard-Jones model, employing parallel tempering Monte Carlo to cover the crystalline, amorphous, and liquid regions of configuration space. Disconnectivity graphs are used to visualize the potential energy landscape in the vicinity of a crystalline geometry and in an amorphous region of configuration space. The crystalline global minimum is separated from the bulk of the minima by a large potential energy gap, leading to broken ergodicity in conventional simulations. Our sampling reveals crystalline global minima that are lower in potential energy than some of the previous candidates. We present equilibrium thermodynamic properties based on parallel tempering simulations, including heat capacities and free energy profiles, which depend explicitly on the crystal structure. We also report equilibrium melting temperatures.     

Relative solvation free energies calculated using an ab initio QM/MM-based free energy perturbation method: dependence of results on simulation length

Molecular dynamics (MD) simulations in conjunction with thermodynamic perturbation approach was used to calculate relative solvation free energies of five pairs of small molecules, namely; (1) methanol to ethane, (2) acetone to acetamide, (3) phenol to benzene, (4) 1,1,1 trichloroethane to ethane, and (5) phenylalanine to isoleucine. Two studies were performed to evaluate the dependence of the convergence of these calculations on MD simulation length and starting configuration. In the first study, each transformation started from the same well-equilibrated configuration and the simulation length was varied from 230 to 2,540 ps. The results indicated that for transformations involving small structural changes, a simulation length of 860 ps is sufficient to obtain satisfactory convergence. In contrast, transformations involving relatively large structural changes, such as phenylalanine to isoleucine, require a significantly longer simulation length (>2,540 ps) to obtain satisfactory convergence. In the second study, the transformation was completed starting from three different configurations and using in each case 860 ps of MD simulation. The results from this study suggest that performing one long simulation may be better than averaging results from three different simulations using a shorter simulation length and three different starting configurations.     

Optimized parallel tempering simulations of proteins

We apply a recently developed adaptive algorithm that systematically improves the efficiency of parallel tempering or replica exchange methods in the numerical simulation of small proteins. Feedback iterations allow us to identify an optimal set of temperatures/replicas which are found to concentrate at the bottlenecks of the simulations. A measure of convergence for the equilibration of the parallel tempering algorithm is discussed. We test our algorithm by simulating the 36-residue villin headpiece subdomain HP-36 where we find a lowest-energy configuration with a root-mean-square deviation of less than 4 A to the experimentally determined structure.     

All-exchanges parallel tempering

An alternative exchange strategy for parallel tempering simulations is introduced. Instead of attempting to swap configurations between two randomly chosen but adjacent replicas, the acceptance probabilities of all possible swap moves are calculated a priori. One specific swap move is then selected according to its probability and enforced. The efficiency of the method is illustrated first on the case of two Lennard-Jones (LJ) clusters containing 13 and 31 atoms, respectively. The convergence of the caloric curve is seen to be at least twice as fast as in conventional parallel tempering simulations, especially for the difficult case of LJ31. Further evidence for an improved efficiency is reported on the ergodic measure introduced by Mountain and Thirumalai [J. Phys. Chem. 93, 6975 (1989)], calculated here for LJ13 close to the melting point. Finally, tests on two simple spin systems indicate that the method should be particularly useful when a limited number of replicas are available.     

Rigid quantum Monte Carlo simulations of condensed molecular matter: water clusters in the n=2-->8 range

The numerical advantage of quantum Monte Carlo simulations of rigid bodies relative to the flexible simulations is investigated for some simple systems. The results show that if high frequency modes in molecular condensed matter are predominantly in the ground state, the convergence of path integral simulations becomes nonuniform. Rigid body quantum parallel tempering simulations are necessary to accurately capture thermodynamic phenomena in the temperature range where the dynamics are influenced by intermolecular degrees of freedom; the stereographic projection path integral adapted for quantum simulations of asymmetric tops is a significantly more efficient strategy compared with Cartesian coordinate simulations for molecular condensed matter under these conditions. The reweighted random series approach for stereographic path integral Monte Carlo is refined and implemented for the quantum simulation of water clusters treated as an assembly of rigid asymmetric tops.     

Low-energy electron scattering by cubane: resonant states and Ramsauer-Townsend features from quantum calculations in the gas phase

Calculations are carried out, using a nonempirical modeling of the interaction potential and solving the quantum scattering coupled channel equations, for low energy electron scattering from cubane (C8H8) molecules in the gas phase. Total integral cross sections are obtained and partial contributions are analyzed for the most important irreducible representations that describe the continuum electron in the Oh molecular symmetry. Several trapping resonances are found and analyzed in terms of the molecular-type features of the resonant electron states associated with them. A Ramsauer-Townsend minimum is also found and its possible behavior related to features of the scattering length as k --> 0.     

Structure of protonated water clusters: low-energy structures and finite temperature behavior

The structure of protonated water clusters H+(H2O)n (n=5-22) are examined by two Monte Carlo methods in conjunction with the OSS2 potential [L. Ojamae, I. Shavitt, and S. J. Singer J. Chem. Phys. 109, 5547 (1998)]. The basin-hopping method is employed to explore the OSS2 potential energy surface and to locate low-energy structures. The topology of the "global minimum," the most stable low-energy structure, changes from single ring to multiple ring to polyhedral cage as the cluster size grows. The temperature dependence of the cluster geometry is examined by carrying out parallel tempering Monte Carlo simulations. Over the temperature range we studied (25-330 K), all water clusters undergo significant structural changes. The trends are treelike structures dominating at high temperature and single-ring structures appearing in slightly lower temperatures. For n> or =7, an additional transition from single ring to multiple rings appears as the temperature decreases. Only for n> or =16 do polyhedral structures dominate the lowest temperature range. Our results indicate very dynamic structural changes at temperature range relevant to atmospheric chemistry and current experiments. The structures and properties of medium-sized protonated clusters in this temperature range are far from their global minimum cousins. The relevance of these findings to recent experiments and theoretical simulations is also discussed.     

ChemInform Abstract: Molecular Geometry and Fluxionality in the AX6E System ClF- 6.

The analysis of the calculated charge distribution and corresponding Laplacian distributions of the ClF₆^- ion in both its equilibrium geometry (C3v) and its fluxional transition state geometry (C2v) confirms the suggestion that AX6E systems should distort from Oh symmetry because of the presence of a stereochemically active lone pair.     

Contributions from atomic p(Se), d(Se), and f(Se) orbitals to absolute paramagnetic shielding tensors in neutral and charged SeHn and some oxides including the effect of methyl and halogen substitutions on sigmap(Se)

Contributions from atomic p(Se), d(Se), and f(Se) orbitals to sigmap(Se) are evaluated for neutral and charged Se*Hn (*=null, +, or -) and some oxides to build the image of the contributions. The effect of methyl and halogen substitutions is also examined employing RrSe*XxOo (*=null, +, or -) where R=H or Me; X=F, Cl, or Br. The p(Se) contributions are larger than 96 % for SeH- (Cinfinityv), SeH2 (C2v), SeH3 + (C3v), SeH3 + (D3h), and SeH4 (Td). Therefore, sigmap(Se) of these compounds can be analyzed based on p(Se). The p(Se) contributions are 79-75 % for SeH4 (TBP), SeH5 + (TBP), SeH5 + (SP), and SeH5 - (SP). Methyl and halogen substitutions increase the contributions by 1-2 % (per Me) and 4-7 % (per X), respectively. The contributions are 92-79 % for H2SeO (Cs), H2SeO2 (C2v), and H4SeO (C2v). The values are similarly increased by the substitutions. Consequently, sigmap(Se) of these compounds can be analyzed based on p(Se) with some corrections by d(Se). The p(Se) contribution of SeH6 (Oh) is 52 %: sigmap(Se: SeH6 (Oh)) must be analyzed based on both p(Se) and d(Se). The contributions for the Me and X derivatives of SeH(6) amount to 86-77 %. Therefore, sigmap(Se) of the derivatives can also be analyzed mainly based on p(Se) with some corrections by d(Se). Contributions from f(Se) are negligible. Contributions from 4p(Se) in vacant orbitals are also considered. A utility program derived from the Gaussian 03 (NMRANAL-NH03G) is applied to evaluate the contributions.     

Extensive all-atom Monte Carlo sampling and QM/MM corrections in the SAMPL4 hydration free energy challenge

We present our predictions for the SAMPL4 hydration free energy challenge. Extensive all-atom Monte Carlo simulations were employed to sample the compounds in explicit solvent. While the focus of our study was to demonstrate well-converged and reproducible free energies, we attempted to address the deficiencies in the general Amber force field force field with a simple QM/MM correction. We show that by using multiple independent simulations, including different starting configurations, and enhanced sampling with parallel tempering, we can obtain well converged hydration free energies. Additional analysis using dihedral angle distributions, torsion-root mean square deviation plots and thermodynamic cycles support this assertion. We obtain a mean absolute deviation of 1.7 kcal mol^?1 and a Kendall’s τ of 0.65 compared with experiment.     

Predictive concept for lone-pair distortions - DFT and vibronic model studies of AXn-(n-3) molecules and complexes (A = NIII to BiIII; X = F-I to I-I; n = 3-6)

The stereochemical and energetic consequences of the lone-pair effect in the title molecules and complexes have been studied by DFT calculations based on a vibronic coupling concept. The anionic complexes were examined as bare entities and, more realistically, in a polarizable charge-compensating solvent continuum. The tendency for distortions of AX3 compounds away from the high-symmetry parent geometry becomes more pronounced the larger the chemical hardness of a molecule and its constituents is; on the other hand, anionic complexes AXn-(n-3) (n = 4-6) become softer and less susceptible to distortion as compared to the corresponding AX3 molecule, the larger the coordination number and the anionic charge are. Thus, while all AX(3) compounds adopt the distorted C3v structure, only very few AX6(3-) species are calculated to deviate from the parent Oh geometry. If a complex possesses a low stabilization energy due to an unfavorable central ion/ligand size ratio, vibronic coupling may even lead to complete dissociation of one (SbF6(3-) --> SbF5(2-) + F-) or more (PF6(3-) --> PF4- + 2F-) ligands. The derived hardness rule perfectly covers the reported structural findings. The calculations indicate that the lone-pair effect is an orbital overlap phenomenon. The interpair repulsion within the valence shell, keeping the average bond distances constant, does not stabilize the distorted with respect to the parent geometry, in disagreement with the VSEPR model.     

Multiple‐Replica Exchange with Information Retrieval

The parallel tempering simulation method was recently extended to allow for possible exchanges between non‐adjacent replicas. We introduce a multiple‐exchange variant which naturally incorporates the information from all replicas when calculating statistical averages, building on the related virtual‐move method of Coluzza and Frenkel (ChemPhysChem 2005 , 6, 1779). The method is extensively tested on three model systems, namely, a Lennard‐Jones cluster exhibiting a finite size phase transition, the Lennard‐Jones fluid, and the 2D ferromagnetic Ising model. In all cases, the present method performs significantly better and converges faster than conventional parallel tempering Monte Carlo simulations. The standard deviations are also systematically decreased with respect to virtual moves.     

Protein structure prediction by tempering spatial constraints

Summary The probability to predict correctly a protein structure can be enhanced through introduction of spatial constraints – either from NMR experiments or from homologous structures. However, the additional constraints lead often to new local energy minima and worse sampling efficiency in simulations. In this work, we present a new parallel tempering variant that alleviates the energy barriers resulting from spatial constraints and therefore yields to an enhanced sampling in structure prediction simulations.     

Homoleptic mononuclear and binuclear osmium carbonyls Os(CO)n(n = 3-5) and Os2(CO)n (n = 8, 9): comparison with the iron analogues

The structures and energetics of the experimentally known Os(CO)n ( n = 3-5), Os2(CO)9, and Os2(CO)8 have been investigated using density functional theory. For Os(CO)5, the lowest-energy structure is the singlet D(3h) trigonal bipyramid. However, the C(4v) square pyramid for Os(CO)5 lies only approximately 1.5 kcal/mol higher in energy, suggesting extraordinary fluxionality. For the coordinatively unsaturated Os(CO)4 and Os(CO)3, a D(2d) strongly distorted tetrahedral structure and a Cs bent T-shaped structure are the lowest-energy structures, respectively. For Os2(CO)9, the experimentally observed singly bridged Os2(CO)8(mu-CO) structure is the lowest-energy structure. A triply bridged Os2(CO)6(mu-CO)3 structure analogous to the known Fe2(CO)9 structure is a transition state rather than a true minimum and collapses to the singly bridged global minimum structure upon following the corresponding normal mode. An unbridged (OC)5Os --> Os(CO)4 structure with a formal Os --> Os dative bond analogous to known stable complexes of the type (R3P)2(OC)3Os --> W(CO)5 is also found for Os2(CO)9 within 8 kcal/mol of the global minimum. The global minimum for the coordinatively unsaturated Os2(CO)8 is a singly bridged (OC)4Os(mu-CO)Os(CO)3 structure derived from the Os2(CO)9 global minimum by loss of a terminal carbonyl group. However, the unbridged structure for Os2(CO)8 observed in low-temperature matrix experiments lies only approximately 1 kcal/mol above this global minimum. In all cases, the triplet structures for these osmium carbonyls have significantly higher energies than the corresponding singlet structures.     

Hierarchical parallelization of divide-and-conquer density functional tight-binding molecular dynamics and metadynamics simulations

Massively parallel divide-and-conquer density functional tight-binding (DC-DFTB) molecular dynamics and metadynamics simulations are efficient approaches for describing various chemical reactions and dynamic processes of large complex systems via quantum mechanics. In this study, DC-DFTB simulations were combined with multi-replica techniques. Specifically, multiple walkers metadynamics, replica exchange molecular dynamics, and parallel tempering metadynamics methods were implemented hierarchically into the in-house Dcdftbmd program. Test simulations in an aqueous phase of the internal rotation of formamide and conformational changes of dialanine showed that the newly developed extensions increase the sampling efficiency and the exploration capabilities in DC-DFTB configuration space.     

A Hundred‐Year‐Old Experiment Re‐evaluated: Accurate Ab Initio Monte Carlo Simulations of the Melting of Radon

State‐of‐the‐art relativistic coupled‐cluster theory is used to construct many‐body potentials for the noble‐gas element radon to determine its bulk properties including the solid‐to‐liquid phase transition from parallel tempering Monte Carlo simulations through either direct sampling of the bulk or from a finite cluster approach. The calculated melting temperature are 200(3) K and 200(6) K from bulk simulations and from extrapolation of finite cluster values, respectively. This is in excellent agreement with the often debated (but widely cited) and only available value of 202 K, dating back to measurements by Gray and Ramsay in 1909.     

Density functional theory study of 11-atom germanium clusters: effect of electron count on cluster geometry

Density functional theory (DFT) at the hybrid B3LYP level has been applied to the germanium clusters Ge(11)(z) (z = -6, -4, -2, 0, +2, +4, +6) starting from eight different initial configurations. The global minimum within the Ge(11)(2-) set is an elongated pentacapped trigonal prism distorted from D(3)(h) to C(2v) symmetry. However, the much more spherical edge-coalesced icosahedron, also of C(2v) symmetry, expected by the Wade-Mingos rules for a 2n + 2 skeletal electron system and found experimentally in B(11)H(11)(2-) and isoelectronic carboranes, is of only slightly higher energy (+5.2 kcal/mol). Even more elongated D(3)(h) pentacapped trigonal prisms are the global minima for the electron-rich structures Ge(11)(4-) and Ge(11)(6-). For Ge(11)(4-) the C(5v) 5-capped pentagonal antiprism analogous to the dicarbollide ligand C(2)B(9)H(11)(2-) is of significantly higher energy (approximately 28 kcal/mol) than the D(3h) global minimum. The C(2v) edge-coalesced icosahedron is also the global minimum for the electron-poor Ge(11) similar to its occurrence in experimentally known 11-vertex "isocloso" metallaboranes of the type (eta(6)-arene)RuB(10)H(10). The lowest energy polyhedral structures computed for the more hypoelectronic Ge(11)(4+) and Ge(11)(6+) clusters are very similar to those found experimentally for the isoelectronic ions E(11)(7-) (E = Ga, In, Tl) and Tl(9)Au(2)(9-) in intermetallics in the case of Ge(11)(4+) and Ge(11)(6+), respectively. These DFT studies predict an interesting D(5h) centered pentagonal prismatic structure for Ge(11)(2+) and isoelectronic metal clusters.     

Multiparametric curve fitting XIV. Modus operandi of the least-squares algorithm MINOPT

Hybrid least-squares algorithm MINOPT for a nonlinear regression is introduced. MINOPT from CHEMSTAT package combines fast convergence of the Gauss-Newton method in a vicinity of minimum with good convergence of gradient methods for location far from a minimum. Quality of minimization and an accuracy of parameter estimates for six selected models are examined and compared with different derivative least-squares methods of five commercial regression packages.     

Structural transitions in the 309-atom magic number Lennard-Jones cluster

The thermal behavior of the 309-atom Lennard-Jones cluster, whose structure is a complete Mackay icosahedron, has been studied by parallel tempering Monte Carlo simulations. Surprisingly for a magic number cluster, the heat capacity shows a very pronounced peak before melting, which is attributed to several coincident structural transformation processes. The main transformation is somewhat akin to surface roughening and involves a cooperative condensation of vacancies and adatoms that leads to the formation of pits and islands one or two layers thick on the Mackay icosahedron. The second transition in order of importance involves a whole scale transformation of the cluster structure and leads to a diverse set of twinned structures that are assemblies of face-centered-cubic tetrahedra with six atoms along their edges, i.e., one atom more than the edges of the 20 tetrahedra that make up the 309-atom Mackay icosahedron. A surface reconstruction of the icosahedron from a Mackay to an anti-Mackay overlayer is also observed, but with a lower probability.