ICSM: An order [formula omitted] method for calculating electrostatic interactions added to TINKER

Generalized image charge solvation model for electrostatic interactions in molecular dynamics simulations of aqueous solutions

Animportant and computationally demanding part of molecular dynamicssimulations is the calculation of long-range electrostatic interactions.Today, the prevalent method to compute these interactions is particlemesh Ewald (PME). The PME implementation in the GROMACS moleculardynamics package is extremely fast on individual GPU nodes. However,for large scale multinode parallel simulations, PME becomes the mainscaling bottleneck as it requires all-to-all communication betweenthe nodes; as a consequence, the number of exchanged messages scalesquadratically with the number of involved nodes in that communicationstep. To enable efficient and scalable biomolecular simulations onfuture exascale supercomputers, clearly a method with a better scalingproperty is required. The fast multipole method (FMM) is such a method.As a first step on the path to exascale, we have implemented a performance-optimized,highly efficient GPU FMM and integrated it into GROMACS as an alternativeto PME. For a fair performance comparison between FMM and PME, wefirst assessed the accuracies of the methods for various sets of inputparameters. With parameters yielding similar accuracies for both methods,we determined the performance of GROMACS with FMM and compared itto PME for exemplary benchmark systems. We found that FMM with a multipoleorder of 8 yields electrostatic forces that are as accurate as PMEwith standard parameters. Further, for typical mixed-precision simulationsettings, FMM does not lead to an increased energy drift with multipoleorders of 8 or larger. Whereas an ≈50 000 atom simulationsystem with our FMM reaches only about a third of the performancewith PME, for systems with large dimensions and inhomogeneous particledistribution, e.g., aerosol systems with water droplets floating ina vacuum, FMM substantially outperforms PME already on a single node.

OxyR is a transcriptional factor, which activates transcription of antioxidant genes. In this study, we constructed a structural model for the full-length OxyR from Shewanella oneidensis MR-1 using threading and comparative modeling techniques. To further study the conformational changes, we also performed molecular dynamics simulations on the activation domain of OxyR. Molecular dynamics simulations were performed using GROMACS force field under periodic boundary conditions. The particle mesh Ewald (PME) method was used to treat long-range electrostatic interactions. The simulation results show that the oxidized form is very stable while the reduced form is quite flexible. Our results suggest that the reduced form provides structural flexibility for disulfide bond formation and which in turn regulates its function.

In this paper, we report a detailed comparison between the popularly used cutoff and Particle Mesh Ewald (PME) methods in terms of the time complexity and the energy convergence of the long-range electrostatic and van der Waals interaction calculations. For the comparison, we performed various calculations on various representative biological molecules, including seven peptides and proteins, eleven oligonucleotides, and three conformations of a nucleotide-sugar. The results provide useful insights into the appropriate choice of the methods (i.e. the cutoff or PME) and that of the cutoff values for the calculations on different kinds of molecules. It has also been demonstrated that for some cases using different cutoff values for calculating the electrostatic and van der Waals interaction energies will be computationally more efficient.

In this article we present two 1000 ps molecular dynamics simulations on the ratμ-glutathione S-transferase dimeric enzyme in complex with the product 1-(S-glutathionyl)-2,4-dinitrobenzene, in a periodic box with explicit solvent molecules, and investigate the effect of long-range electrostatics models on the structure and dynamics of the dimer and its components. One simulation used the standard cutoff method (10Å), whilst the other used the particle-mesh Ewald (PME) method. We monitored the root mean-square atomic deviation (RMSD) from the initial crystal structure to examine the convergence of both simulations, as well as several other structural parameters such as the distance between active sites, rigid body rotation between domains in subunits, radius of gyration, B-factors, number of hydrogen bonds and salt bridges and solvent-accessible surface area. For example, with the PME method, the dimer structure remains much closer to the initial crystallographic structure with an average RMSD of 1.3Å ± 0.1Å and 1.0Å ± 0.1Å for all heavy and backbone atoms, respectively, in the last 200 ps; the respective values for the cutoff simulation are 4.7Å ± 0.3Å and 4.2Å ± 0.3Å. The large deviations observed in the cutoff simulation severely affected the stability of the enzyme dimer and its complex with the bound product. This finding is contrary to that found in a similar study of the monomelic protein ubiquitin [Fox, T. & Kollman, P. A. Proteins Struct. Fune. Genet. 25, 315–334 (1996)]. Unlike the earlier published work, the present study provides evidence that the standard cutoff method is not generally valid for the study of protein complexes, or their subunits.

Ab initio force field for simulations of proteins and nucleic acids: Application to molecular dynamics simulations in aqueous solution

Ubiquitin, a ubiquitous small protein found in all living organisms, is crucial for tagging proteins earmarked for degradation and holds pivotal importance in biomedicine. Protein functionality is intricately linked to its structure. To comprehend the impact of diverse temperatures on ubiquitin protein structure, our study delved into the energy landscape, hydrogen bonding, and overall structural stability of ubiquitin protein at varying temperatures. Through meticulous analysis of root mean square deviation and root mean square fluctuation, we validated the robustness of the simulation conditions employed. Within our simulated system, the bonding energy and electrostatic potential energy exhibited linear augmentation, while the van der Waals energy demonstrated a linear decline. Additionally, our findings highlighted that the α-Helix secondary structure of the ubiquitin protein gradually transitions toward helix destabilization under high-temperature conditions. The secondary structure of ubiquitin protein experiences distinct changes under varying temperatures. The outcomes of our molecular simulations offer a theoretical framework that enhances our comprehension of how temperature impacts the structural stability of ubiquitin protein. These insights contribute not only to a deeper understanding of iniquity's behavior but also hold broader implications in the realm of biomedicine and beyond. All the MD simulations were performed using the GROMACS software with GROMOS96 force field and SPC for water. The ubiquitin protein was put in the center of a cubic box with a length of 8nm, a setting that allowed > 0.8nm in the minimal distance between the protein surface and the box wall. To remove the possible coordinate collision of the configurations, in the beginning, the steepest descent method was used until the maximum force between atoms was under 100kJ/mol·nm with a 0.01nm step size. Minimization was followed by 30ps of position-restrained MD simulation. The protein was restrained to its initial position, and the solvent was freely equilibrated. The product phase was obtained with the whole system simulated for 10ns without any restraint using an integral time step of 1fs with different temperatures. The cutoff for short-range electronic interaction was set to 1.5nm. The long-range interactions were treated with a particle-mesh Ewald (PME) method with a grid width of 1.2nm.

Replica exchange molecular dynamics simulations were performed to investigate the effects of different electrostatic treatments on the structure and thermodynamics of a small beta-hairpin forming peptide. Three different electrostatic schemes were considered: regular cutoffs, generalized reaction field (GRF), and particle mesh Ewald (PME), with the peptide modeled using OPLS/AA all-atom force field with explicit TIP3P water. Both the GRF and PME methods yielded results consistent with experiment, with free energy surfaces displaying a single minimum corresponding to the native beta-hairpin structure. In contrast, use of straight cutoffs led to the population of an additional local minimum corresponding to nonhairpin conformations that compete with the formation of the native beta-hairpin at low temperatures. This extra minimum would not be apparent in conventional constant-temperature molecular dynamics simulations run for a few nanoseconds. This result points to the critical need of careful sampling of conformational space to assess the quality of different numerical treatments of long-range forces. While differences emerged in the nature of the unfolded states populated using PME and GRF approaches, simulations on the beta-hairpin forming peptide and on two additional control peptides indicate that the GRF treatment of electrostatics offers a satisfactory compromise between accuracy and computational speed for the identification of low-energy conformations. A GRF-based approach emerges as a viable means for treating larger biological systems that would be prohibitively costly to simulate using PME methods.

A Cuda Fast Multipole Method with Highly Efficient M2L Far Field Evaluation

Though computational molecular dynamics is an effective tool for nano-scale phenomenon analysis, computational costs associated with its computer simulation are extremely high. There are two major computational steps associated with computer simulation of dynamics of molecular structures. They are calculation of interatomic forces and formation and solution of the equations of motion. Currently, these two computational steps are treated separately in most commonly-used methods. For example, Fast Multipole Method (FMM) and Cell Multipole Method (CMM) have been used for calculation of interatomic forces, and Cartesian Coordinate Method (CCM) and Internal Coordinate Molecular Dynamics Method (ICMD) have been created for the formation and solution of equations of motion of an atomistic molecular system. In this paper, a new procedure is presented through a proper integration between multibody molecular algorithms (MMA) and fast multipole methods to improve computational efficiency for computer simulation of the dynamical behaviors of multibody molecular structures in polymers and biopolymers. For the computational costs associated with interatomic forces, a fast multipole method is used to calculate the interatomic forces due to the potentials. For the computational costs associated with formation and solution of equations of motion, a multibody molecular algorithm developed by the author in his previous work will be utilized to integrate with fast multipole methods. The algorithm significantly improves computational efficiency when comparing with its counterpart procedures. The fast multipole method begins by scaling all atoms into a box with coordinate ranges to ensure numerical stability of subsequent operations. The parent box is then divided into half in the direction of each Cartesian axis and each child box is then subdivided to form a computational family tree. The flow of calculations is carried out along the tree structure with five passes. The fast multipole method has been improved and modified to achieve better effectiveness and higher efficiency since it was created. The multibody molecular algorithm starts with numbering subsets, forming bond graph, and developing three computing passes along the tree structure of an atomistic molecular system. Computing data flows in the fast multipole method and the multibody molecular algorithm will properly line up with the parent-child recursive relationship along the configuration of the tree structure due to linear recursive natures of both fast multipole method and multibody molecular algorithm. Then the time spent on the recursive simulation passes in the fast multipole method for computing forces may overlap with the time spent on the three recursive computational passes in the multibody molecular algorithm for forming and solving equations of motion.

An Atomic Model for the Pleated β-Sheet Structure of A β Amyloid Protofilaments

Ionoids: Macrocyclic, organic and small inorganic building blocks electrostatically self-assemble in solution. The formed nano-objects feature a mobile interior. The molecules inside the aggregates seem to rotate freely on certain spatial positions. The hydrodynamic radius of the nano-assemblies, however, is precisely defined. Many societal challenges at the beginning of the 21st century lead to an apparent and growing need for functional materials and novel ways of materials synthesis and assembly. Rising to the challenge, the utilization of small, self-assembling building blocks for the bottom-up construction of new types of polymers and nanostructures has enjoyed increasing popularity among materials researchers in the recent past. Supramolecular materials like foldamers, surface films, nanoparticles, etc. are created by exploiting noncovalent forces1 leading to an ordered arrangement of nanoscale building blocks.2 In the search for new polymers based on noncovalent molecular forces, we are motivated by the idea of supramolecular or even polymer-like structures by self-assembly of small ionic monomers, merely formed from electrostatic and solvent interactions. We, in particular, focus on applying intrinsically nondirected electrostatic interactions to create nonetheless well- organized supramolecular structures.3 To achieve this aim, we here make use of a recently highlighted, multicationic molecular box (14+; Figure 1 a), developed by Sessler and co-workers, which has already been applied in self-assembly applications of metal–organic rotaxane frameworks.4 The molecular box 14+ seems to facilitate structure generation through self- assembly, indicated by, for example, the surprising ability to build supramolecular necklaces through incorporation of electron-rich aromatic guest molecules.5 However, not focused on short-ranged interactions but on new applications of long-ranged electrostatic forces at the nanoscale,6 we combine 14+ with small di-anionic salts, (KSO3)2CH2 (22−) and (KSO3)2NO (32−), to generate a defined assembly of cations and anions in solution (Figure 1 a). However, self-assembly due to electrostatic attraction is complex to describe qualitatively as well as quantitatively, because not only Coulomb forces but also entropy changes due to counterion release, solvation of ions, depletion forces, etc., contribute to the total free energy of a system. A central question in any application of solution-based ionic self-assembly therefore is: does any kind of ordered arrangement of the ionic building blocks take place?7 Intuitively one would assume the distribution of ions in solution to be of random, rather homogeneous nature. a) Molecular structure of 14+ (left) and 22−/32− (top right). b) A typical monoexponential fit of an autocorrelation function gained from DLS on 22−/14+ ionoids. c) Hydrodynamic radii of the 22−/14+ ionoids plotted against time, as derived from DLS measurements. t=0 days refers to the date of initial preparation. The upper and lower lines indicate error estimates derived from the fits of the autocorrelation functions. In the following, we show dynamic light scattering (DLS) data, which clearly illustrate that low-concentrated mixtures of anionic 22− and cationic 14+ self-assemble in solution to highly defined, monodisperse ionic clusters. In this special case, the distribution of ions in solution is therefore not at all random but well-ordered due to long-range electrostatic correlations. Performing molecular dynamics (MD) and Monte Carlo (MC) simulations, we further elucidate the internal constitution of the ionic clusters and show that a long-ranged order may be induced by electrostatic correlations of the cationic macrocycles. The simulations indicate that the monodisperse objects are constituted of relatively loosely bound ions that, however, occupy certain fixed coordinates in the clusters for at least orders of nanoseconds. The reliability of the computational data is additionally supported by double electron–electron resonance (DEER)8 spectroscopy data. Taking long-ranged (DLS, MC) and molecular (MD, DEER) insights together, we demonstrate, to the best of our knowledge for the first time, that size-controlled self-assembly through electrostatic interaction between small ions in solution is possible. The so formed objects can be regarded as loosely bound ion-based colloids, thus we propose the name ionoids. These ionoids are peculiar as they result from ion–ion correlations of moderately charged constituents at low concentrations, while long-range electrostatic correlations have so far only been discussed in the context of macroions.9 DLS experiments were performed on a system containing a mixture of 22−/14+ at an approximate ratio of 3 mM:1 mM in DMSO/88 % aqueous glycerol (1:1). Note that the exact concentration of 22− within the light-scattering sample is difficult to determine, as the solubility of the pure component 22−, in contrast to 32−, is only 3 mM in our solvent mixture. The measurements revealed monodisperse small objects in the solution, after an initial period of several days where no correlation could be detected. The fit of one correlation function measured at a scattering angle of 90 ° is depicted in Figure 1 b. Note the monoexponential nature of this fit, which indicates monodispersity of the aggregates in the solution. The intercept of the correlation function is unusually low, due to the small signal-to-noise ratio of <3, that is, the scattered intensity from the pure solvent mixture is nearly 40 % of the overall scattered intensity. As, in the absence of 22−, we found no correlation at all but only the DLS signature of the pure solvent mixture, it is at hand to deduce that these small objects stem from self- assembly of 22−–14+ clusters. The hydrodynamic radius of the self-assembled 22−–14+ cluster is quite small, that is, only 6.2 nm±0.2 nm can be observed 10 days after sample preparation, thus, growing very slow (+2 nm after 140 days; see Figure 1 c) and fostering an analogy with charge-stabilized colloids. It is important to note that all correlation functions measured over time can be fitted by a single exponential decay, and the resulting diffusion coefficient is independent of the scattering angle, which underlines the monodispersity in aggregate size. Thus, one can conclude that molecules 22− are necessary to overcome like-charge repulsion and to bring together 3–4 molecules 14+ (estimated from the size of approx. 6–7 nm) in an extraordinarily uniform manner. To support the reliability of the DLS data and to elucidate the internal constitution of the ionic clusters, we conducted all-atom MD simulations corresponding to a 32−/14+ system at a molar ratio of 6:2. In order to check the reliability of the MD simulation, we supplement it with DEER data and therefore replaced 22− by paramagnetic 32−. (This change may have an influence on special interactions in the ionoids, yet long-range electrostatic interactions remain unchanged. Note that the concentration in the simulations was three times higher than in the DEER measurements, to reduce the computation time). Starting from a random distribution of all molecules, we let the system equilibrate for 0.2 ns, and subsequently we extracted the trajectories for 8 ns with the aim of calculating the radial distribution function, g(r), of 32−. Figure 2 a shows a snapshot of the simulation after 2 ns, illustrating the molecular arrangement of the ions in such an ionoid-type cluster. Figure 2 b depicts the radial distance distribution of 32− summed up over 8 ns of simulation and over different time intervals of 1 ns. Prominent features in the overall g(r) at 1.9 nm, 2.4 nm, 3.5 nm are observable. Thus, the distribution of 32− in the simulated system is far from being random. Also, the distributions summed up in each case over 1 ns of simulation indicate favored population of distinct distances in the 32−/14+ system. At first sight, it is astonishing that the maxima in g(r) occur with a spacing of approximately 0.5 nm, indicating a periodic arrangement of the spin probes (Figure 2 b) . a) Snapshot after 2 ns of MD of the system 32−/14+ at 6:2 molar ratio. b) Pair correlation function of 32− molecules gained from MD simulation of the system 32−/14+ (6:2) summed up over different time intervals. The overall distance distribution is shown in orange (—). Note that over time, certain maxima in the distributions appear frequently, as indicated by the dashed lines. The black distribution depicts the distances gained from DEER on the same system (—). c) The Monte Carlo-derived radial distribution function for 32− calculated for two 14+ and six 32− molecules (—). The box size was 263 nm3. As can be seen, there are no distinct maxima in the distance distribution. Hence, the aggregation of 32− and 14+ is not based on Coulomb interaction and excluded volume effects alone. Solvent entropy may possibly contribute significantly to the generation of the ionoids, as the MD-derived distributions lead to reasonable agreements with the experimental data. The 14+–14+ distribution (—) reveals a long-range correlation at approximately 15 nm. For more details, see the Supporting Information. Although the inner constitution of the 32−–14+ aggregates could be assessed with MD, the global size of the so formed aggregates cannot be deduced from the simulation and compared with RH because of a box size of only 83 nm3. Consequently, we performed Monte Carlo (MC) simulations with a box size of 263 nm3 to elucidate the individual cluster size. These simulations included electrostatic potentials as well as excluded volume effects. However, the MC simulation showed that, due to mere electrostatic interactions and excluded volume, no distinct short-range association takes place (Figure 2 c). Yet, the MC simulations do reveal a surprising long-range correlation of 14+ at about 15 nm, suggesting a defined size of the observed ionic clusters, electrostatic and excluded volume effects to be essential for the formation of the observed ionic clusters (Figure 2 c and the Supporting Information), because the experimentally determined radius (6.2 nm–8.2 nm) and the size suggested from MC simulations (R≈7.5 nm) are in good agreement. Yet the distinct 32−–32− distances, as gained from MD, seem to be a consequence of additional entropic or solvent effects next to mere electrostatic effects. To further support the reliability of the computational data, we performed electron paramagnetic resonance (EPR)-based nanoscale distance measurement (DEER) experiments on a freeze-quenched 32−/14+ system to extract the experimental distance distribution of 32−. The resulting distance distribution is depicted in Figure 2 b (black curve). It shows the population of all 32−–32− distances of 6 mM molecules 32− around 2 mM molecules 14+ and can be compared directly with the g(r) gained from MD. Taking into account that the simulation features significantly higher concentrations than the experimentally measured system (the 14+ concentration was three times higher than in the experiment to increase the speed of calculation) and the restricted number of simulated 14+ and 32− molecules, the simulation-derived and the experimental distribution overall are in reasonable agreement. However, there are no distinct dipolar modulations observable in the DEER time traces (see the Supporting Information). Hence, the extracted distances (by Tikhonov regularization) have to be treated with caution. Yet, as they are clearly indicative of nonhomogeneous/nonrandom distributions of 32− and in addition are even in reasonably good agreement with the simulated distributions, one can assume that the distances extracted from the MD simulation resemble the situation inside the ionoids to a certain degree. Thus, judging from the simulation, the situation as shown in Figure 2 a can be taken as a representative snapshot of the distribution of ions in the solution. As it does not change significantly over time, it seems like the molecules 14+ are forced into a short-distance arrangement in order to minimize the potential electrostatic energy of the whole system. Divalent molecules 32−, thereby, could be acting as 'ionic glue' between macrocycles 14+. In doing so, the positions of local maxima in the distance distribution as derived from the MD simulation do not change significantly over time. This means that the spatial conformational ensemble covered seems to be relatively small. Once 32− is 'trapped' in an ensemble, it remains at a certain position for at least several nanoseconds. Continuous wave (CW)-EPR data at room temperature heavily support this hypothesis: the 32− probes still rotate relatively freely, although their rotational mobility is slightly restricted when small amounts of 14+ are added to a 32− solution (see the Supporting Information for the CW-EPR spectra and corresponding spectral simulations). At the same time, they still favor certain coordinates in the ionoids because of electrostatic attraction to 14+ and repulsion among themselves. It is likely that the combination of larger organic and small, yet doubly charged, inorganic building blocks is necessary for the emergence of ordered ionic assemblies, which persist in solution over time. Also, solvent effects could play a role, as DMSO is known to feature anomalous properties because of self-association.10 Taking all the experimental and computational data into account, one can conclude that the self-assembled objects feature a defined radius, although, on a molecular scale, their interior seems not to be strongly fixed: 'gluing' 32− molecules can still rotate relatively undisturbed. Yet, on the supramolecular scale, certain defined coordinates inside the aggregates are preferentially populated for at least some nanoseconds. Hence, the position of the 'glue′ is not arbitrary. From a physico-chemical point of view it is highly intriguing that an intrinsically nondirected force like electrostatic attraction/repulsion leads to the formation of highly defined nanoscopic objects. One can even control the initial size of the self-assembled ionoids precisely and keep them stable in solution for months. Thus, by choosing the right ions and counterions, their distribution in solution is not more or less homogeneous as one would expect intuitively for electrolyte solutions.10–11 Contrarily, they assemble into discrete arrangements as has been observed for macroions such as in DNA-condensation years ago .12 As the observed ionoids are based on noncovalent interaction among constituents, they probably are sensitive towards even minor changes in the solvent properties, ionic strength of the solution, and concentration of the constituents. These factors are likely to influence their composition and size. Furthermore, we expect that the observed growth of the ionoids will be fostered by a decrease in solvent viscosity, and that their occurrence is limited to certain solvents. Altogether, we have presented experimental and computational evidence for the formation of electrostatically self-assembled supramolecular structures in solution. The nano-assemblies have a well-defined size, and molecules in their interior populate, at least transiently, certain fixed positions in the cluster, while keeping a certain rotational freedom. The so formed ionoids can possibly open new routes towards the design and self-assembly of novel materials by exploiting electrostatic interactions between small molecules which have not attracted much attention in this field of research so far. In addition, these nano-objects are intriguing from a theoretical point of view, as long-range electrostatic correlations in solution have been postulated but, to the best of our knowledge, have not yet been observed directly for low-molecular electrolytes (unlike macroions).13 From the slow but steady size increase over the course of two months (Figure 1 c) one could even speculate that ionoids could be some sort of solvated precursor for precipitation or even crystallization nuclei. More experimental evidence would be needed to further substantiate such a picture, which currently is a mere speculation. Considering the small size of the ionoids, it is likely that X-ray and neutron scattering experiments will yield substantial insight into the specific properties of these solution structures. Such measurements are currently on our agenda. The responsiveness of ionoids towards stimuli such as concentration, solvent-mixture or charge ratio could also become interesting in the course of the search for new responsive materials. We thank Christian Bauer for technical support and Prof. Hans W. Spiess for continuing support. D.K. acknowledges support by the Gutenberg Academy of the University of Mainz. This work was financially supported by the Deutsche Forschungsgemeinschaft (DFG) under grant number HI 1094/2-1 and the Max Planck Graduate Center with the University of Mainz (MPGC). Detailed facts of importance to specialist readers are published as "Supporting Information". Such documents are peer-reviewed, but not copy-edited or typeset. They are made available as submitted by the authors. Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.

Computational methods for the study of collision processes in glow discharge plasmas

An efficient implementation of the canonical molecular dynamics simulation using the reversible reference system propagator algorithm (r-RESPA) combined with the particle mesh Ewald method (PMEM) and with the macroscopic expansion of the fast multipole method (MEFMM) was examined. The performance of the calculations was evaluated for systems with 3000, 9999, 30,000, 60,000, and 99,840 particles. For a given accuracy, the optimal conditions for minimizing the CPU time for the implementation of the Ewald method, the PMEM, and the MEFMM were first analyzed. Using the optimal conditions, we evaluated the performance and the reliability of the integrated methods. For all the systems examined, the r-RESPA with the PMEM was about twice as fast as the r-RESPA with the MEFMM. The difference arose from the difference in the numerical complexities of the fast Fourier transform in the PMEM and from the transformation of the multipole moments into the coefficients of the local field expansion in the MEFMM. Compared with conventional methods, such as the velocity-verlet algorithm with the Ewald method, significant speedups were obtained by the integrated methods; the speedup of the calculation was a function of system size, and was a factor of 100 for a system with 3000 particles and increased to a factor of 700 for a system with 99,840 particles. These integrated calculations are, therefore, promising for realizing large-scale molecular dynamics simulations for complex systems. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 201–217, 2000

Version 9 of the Amber simulation programs includes a new semi-empirical hybrid QM/MM functionality. This includes support for implicit solvent (generalized Born) and for periodic explicit solvent simulations using a newly developed QM/MM implementation of the particle mesh Ewald (PME) method. The code provides sufficiently accurate gradients to run constant energy QM/MM MD simulations for many nanoseconds. The link atom approach used for treating the QM/MM boundary shows improved performance, and the user interface has been rewritten to bring the format into line with classical MD simulations. Support is provided for the PM3, PDDG/PM3, PM3CARB1, AM1, MNDO, and PDDG/MNDO semi-empirical Hamiltonians as well as the self-consistent charge density functional tight binding (SCC-DFTB) method. Performance has been improved to the point where using QM/MM, for a QM system of 71 atoms within an explicitly solvated protein using periodic boundaries and PME requires less than twice the cpu time of the corresponding classical simulation.