A united-residue force field for off-lattice protein-structure simulations. II. Parameterization of short-range interactions and determination of weights of energy terms by Z-score optimization

Continuing our work on the determination of an off-lattice united-residue force field for protein-structure simulations, we determined and parameterized appropriate functional forms for the local-interaction terms, corresponding to the rotation about the virtual bonds (U_tor), the bending of virtual-bond angles (U_b), and the energy of different rotameric states of side chains (U_rot). These terms were determined by applying the Boltzmann principle to the distributions of virtual-bond torsional and virtual-bond angles and side-chain rotameric states, respectively, calculated from a data base of 195 high-resolution nonhomologous proteins. The complete energy function was constructed by combining the individual energy terms with appropriate weights. The weights were determined by optimizing the so-called Z-score value (which is the normalized difference between the energy of the native structure and the mean energy of non-native structures) of the histidine-containing phosphocarrier protein from Streptococcus faecalis (1PTF; an 88-residue α + β protein). To accomplish this, a database of C^α patterns was created using high-resolution nonhomologous protein structures from the Protein Data Bank, and the distributions of energy components of 1PTF were obtained by threading its sequence through ∼500 randomly chosen C^α-patterns from the X-ray structures in the PDB, followed by energy minimization, with the energy function incorporating initially guessed weights. The resulting minimized energies were used to optimize the Z-score value of 1PTF as a function of the weights of the various energy terms, and the new weights were used to generate new energy-component distributions. The process was iterated, until the weights used to generate the distributions and the optimized weights were self-consistent. The potential function with the weights of the various energy terms obtained by optimizing the Z-score value for 1PTF was found to locate the native structures of other test proteins (within an average RMS deviation of 3 Å): calcium-binding protein (4ICB), ubiquitin (1UBQ), α-spectrin (1SHG), major cold-shock protein (1MJC), and cytochrome b₅ (3B5C) (which included α and β structures) as distinctively lowest in energy in similar threading experiments. © 1997 by John Wiley & Sons, Inc. J Comput Chem 18: 874–887, 1997     

Molecular dynamics with the united-residue model of polypeptide chains. I. Lagrange equations of motion and tests of numerical stability in the microcanonical mode

The Lagrange formalism was implemented to derive the equations of motion for the physics-based united-residue (UNRES) force field developed in our laboratory. The C(alpha)...C(alpha) and C(alpha)...SC (SC denoting a side-chain center) virtual-bond vectors were chosen as variables. The velocity Verlet algorithm was adopted to integrate the equations of motion. Tests on the unblocked Ala(10) polypeptide showed that the algorithm is stable in short periods of time up to the time step of 1.467 fs; however, even with the shorter time step of 0.489 fs, some drift of the total energy occurs because of momentary jumps of the acceleration. These jumps are caused by numerical instability of the forces arising from the U(rot) component of UNRES that describes the energetics of side-chain-rotameric states. Test runs on the Gly(10) sequence (in which U(rot) is not present) and on the Ala(10) sequence with U(rot) replaced by a simple numerically stable harmonic potential confirmed this observation; oscillations of the total energy were observed only up to the time step of 7.335 fs, and some drift in the total energy or instability of the trajectories started to appear in long-time (2 ns and longer) trajectories only for the time step of 9.78 fs. These results demonstrate that the present U(rot) components (which are statistical potentials derived from the Protein Data Bank) must be replaced with more numerically stable functions; this work is under way in our laboratory. For the purpose of our present work, a nonsymplectic variable-time-step algorithm was introduced to reduce the energy drift for regular polypeptide sequences. The algorithm scales down the time step at a given point of a trajectory if the maximum change of acceleration exceeds a selected cutoff value. With this algorithm, the total energy is reasonably conserved up to a time step of 2.445 fs, as tested on the unblocked Ala(10) polypeptide. We also tried a symplectic multiple-time-step reversible RESPA algorithm and achieved satisfactory energy conservation for time steps up to 7.335 fs. However, at present, it appears that the reversible RESPA algorithm is several times more expensive than the variable-time-step algorithm because of the necessity to perform additional matrix multiplications. We also observed that, because Ala(10) folds and unfolds within picoseconds in the microcanonical mode, this suggests that the effective (event-based) time unit in UNRES dynamics is much larger than that of all-atom dynamics because of averaging over the fast-moving degrees of freedom in deriving the UNRES potential.     

Implementation of the analytic energy gradient for the combined time-dependent density functional theory/effective fragment potential method: application to excited-state molecular dynamics simulations

Excited-state 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 linear-response time-dependent 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 excited-state MD of aqueous formaldehyde in the n-π* state. The calculated π*→n transition energy and solvatochromic shift are in good agreement with other theoretical results.     

Molecular dynamics with the united-residue model of polypeptide chains. II. Langevin and Berendsen-bath dynamics and tests on model alpha-helical systems

The implementation of molecular dynamics (MD) with our physics-based protein united-residue (UNRES) force field, described in the accompanying paper, was extended to Langevin dynamics. The equations of motion are integrated by using a simplified stochastic velocity Verlet algorithm. To compare the results to those with all-atom simulations with implicit solvent in which no explicit stochastic and friction forces are present, we alternatively introduced the Berendsen thermostat. Test simulations on the Ala(10) polypeptide demonstrated that the average kinetic energy is stable with about a 5 fs time step. To determine the correspondence between the UNRES time step and the time step of all-atom molecular dynamics, all-atom simulations with the AMBER 99 force field and explicit solvent and also with implicit solvent taken into account within the framework of the generalized Born/surface area (GBSA) model were carried out on the unblocked Ala(10) polypeptide. We found that the UNRES time scale is 4 times longer than that of all-atom MD simulations because the degrees of freedom corresponding to the fastest motions in UNRES are averaged out. When the reduction of the computational cost for evaluation of the UNRES energy function is also taken into account, UNRES (with hydration included implicitly in the side chain-side chain interaction potential) offers about at least a 4000-fold speed up of computations relative to all-atom simulations with explicit solvent and at least a 65-fold speed up relative to all-atom simulations with implicit solvent. To carry out an initial full-blown test of the UNRES/MD approach, we ran Berendsen-bath and Langevin dynamics simulations of the 46-residue B-domain of staphylococcal protein A. We were able to determine the folding temperature at which all trajectories converged to nativelike structures with both approaches. For comparison, we carried out ab initio folding simulations of this protein at the AMBER 99/GBSA level. The average CPU time for folding protein A by UNRES molecular dynamics was 30 min with a single Alpha processor, compared to about 152 h for all-atom simulations with implicit solvent. It can be concluded that the UNRES/MD approach will enable us to carry out microsecond and, possibly, millisecond simulations of protein folding and, consequently, of the folding process of proteins in real time.     

Theoretical study of Rb2 in He(N): potential energy surface and Monte Carlo simulations

An analytical potential energy surface for a rigid Rb? in the 3Σ(u)? state interacting with one helium atom based on accurate ab initio computations is proposed. This 2-dimensional potential is used, together with the pair approximation approach, to investigate Rb? attached to small helium clusters He(N) with N = 1-6, 12, and 20 by means of quantum Monte Carlo studies. The limit of large clusters is approximated by a flat helium surface. The relative orientation of the dialkali axis and the helium surface is found to be parallel. Dynamical investigations of the pendular and of the in-plane rotation of the rigid Rb? molecule on the surface are presented.     

Temperature dependence of the v(NH) region in the Raman spectrum of imidazole single crystals

The Raman spectrum of imidazole single crystals has been studied in the v(NH) region. The temperature dependence of the band intensity has been discussed. A possible interpretation of the observed bands is proposed.     

Frozen density functional free energy simulations of redox proteins: computational studies of the reduction potential of plastocyanin and rusticyanin

The evaluation of reduction potentials of proteins by ab initio approaches presents a major challenge for computational chemistry. This is addressed in the present investigation by reporting detailed calculations of the reduction potentials of the blue copper proteins plastocyanin and rusticyanin using the QM/MM all-atom frozen density functional theory, FDFT, method. The relevant ab initio free energies are evaluated by using a classical reference potential. This approach appears to provide a general consistent and effective way for reproducing the configurational ensemble needed for consistent ab initio free energy calculations. The FDFT formulation allows us to treat a large part of the protein quantum mechanically by a consistently coupled QM/QM/MM embedding method while still retaining a proper configurational sampling. To establish the importance of proper configurational sampling and the need for a complete representation of the protein+solvent environment, we also consider several classical approaches. These include the semi-macroscopic PDLD/S-LRA method and classical all-atom simulations with and without a polarizable force field. The difference between the reduction potentials of the two blue copper proteins is reproduced in a reasonable way, and its origin is deduced from the different calculations. It is found that the protein permanent dipole tunes down the reduction potential for plastocyanin compared to the active site in regular water solvent, whereas in rusticyanin it is instead tuned up. This electrostatic environment, which is the major effect determining the reduction potential, is a property of the entire protein and solvent system and cannot be ascribed to any particular single interaction.     

The inversion of spectroscopic data to the diatomic effective Hamiltonian containing the generalized potential energy function. Ground state of PbO

A method of inversion of spectroscopic data of a diatomic molecule to the effective Hamiltonian containing adiabatic and non-adiabatic corrections is reported. The method is based on the use of a previously suggested generalized potential energy function with correct long-range part. The potential energy function for the X¹Σ⁺ state of PbO is calculated by inversion of the infrared and microwave spectra.     

Measurements and simulations of high energy O(3P) + Ar(1S) angular scattering: single and multi-collision regimes

We present differential angular cross sections for O(3P) + Ar(1S) scattering at collision energies near 90 kcal mol(-1) (approximately 8 km s(-1) relative velocity) from molecular beam measurements and high-level theoretical calculations. Beams of hyperthermal O(3P) are now being used to investigate novel gas-phase and gas-surface chemistries, and the comparison of theory and measurements on this simple system will be a stringent test of the experimental methodology. Potential energy curves were generated for O(3P) + Ar(1S) using a large cc-pVQZ basis within a valence multi-configuration plus perturbation theory treatment. These curves were then used in quantum scattering calculations to generate differential cross sections. Agreement between experiment and theory is excellent. In addition to these comparisons, the cross sections were used in direct simulation Monte Carlo calculations to investigate effects of increasing the Ar flux above the "single-collision" regime. As the Ar flux increases, the observed differential angular cross sections change in two ways. In addition to the main "single-scatter" peak along the incident O-atom beam direction, a secondary O-atom peak appears in the direction of the incident Ar beam, and the multiple-scattered O-atom translational energy starts to reflect the energy of the relatively slow moving Ar beam.     

CCl和CCl2(1A1)的结构与势能函数

利用Gaussian03软件包,采用多种方法和多种基组对CCl和CCl2分子的基态结构进行优化计算,优选出B3P86/6-311+ G(3 df)方法对CCl分子进行计算得到基态为X2Ⅱ、键长RCCl=0.164 42 nm,谐振频率we=886.3062 cm-1;优选出B3P86/6-311G( 2df)方法对CC...     

Small imines and oximes as model compounds in the optimization of a consistent force field potential energy function

Parameters have been optimized for a potential energy function to be used in molecular mechanics calculations of small imines as a preliminary step to calculations on larger systems. A consistent force field (CFF) program was used, and a new damping algorithm due to Sundius was introduced in the optimization procedure. Optimization of parameters has been done on structural and vibrational data of five small imines and one oxime. The quality of the derived potential energy function is examined by calculations on larger oximes.     

Efficiency of alchemical free energy simulations. I. A practical comparison of the exponential formula, thermodynamic integration, and Bennett's acceptance ratio method

We investigate the relative efficiency of thermodynamic integration, three variants of the exponential formula, also referred to as thermodynamic perturbation, and Bennett's acceptance ratio method to compute relative and absolute solvation free energy differences. Our primary goal is the development of efficient protocols that are robust in practice. We focus on minimizing the number of unphysical intermediate states (λ-states) required for the computation of accurate and precise free energy differences. Several indicators are presented which help decide when additional λ-states are necessary. In all tests Bennett's acceptance ratio method required the least number of λ-states, closely followed by the "double-wide" variant of the exponential formula. Use of the exponential formula in only strict "forward" or "backward" mode was not found to be competitive. Similarly, the performance of thermodynamic integration in terms of efficiency was rather poor. We show that this is caused by the use of the trapezoidal rule as method of numerical quadrature. A systematic study focusing on the optimization of thermodynamic integration is presented in a companion paper.     

Temperature and strain dependence of the raman depolarization in poly(dimethylsiloxane)

Large changes in depolarization ratio of Raman bands in poly(dimethylsiloxane) have been observed as a function of temperature and strain. Analysis in terms of rotational isomeric state populations yields a configurational energy difference of 1000 ± 140 cal/mole.     

Comment on "Free energy simulations of single and double ion occupancy in gramicidin A" [J. Chem. Phys. 126, 105103 (2007)

In a recent article published by Bastug and Kuyucak [J. Chem. Phys.126, 105103 (2007)] investigated the microscopic factors affecting double ion occupancy in the gramicidin channel. The analysis relied largely on the one-dimensional potential of mean force of ions along the axis of the channel (the so-called free energy profile of the ion along the channel axis), as well as on the calculation of the equilibrium association constant of the ions in the channel binding sites. It is the purpose of this communication to clarify this issue.     

Temperature dependence of the free energy of sorption ofn-alkanes and energy contribution of their methylene groups on columns with fullerene C60

The partial molar free energy, enthalpy, and entropy of sorption of C₁₁−C₂₃ n-alkanes were calculated on the basis of the GC data obtained on the glass capillary column coated with fullerene C₆₀ (Ful-60) as stationary phase. The thermodynamic parameters ofn-alkane sorption on a column with Ful-60 and a fused silica capillary column with polydimethylsiloxane OV-1 were determined and compared. The enthalpy-entropy compensation effect for the sorption ofn-alkanes on Ful-60 and OV-1 was found. A linear dependence of the partial molar free energy ofn-alkane sorption on the temperature of analysis and carbon chain length was found. The free energy contributions of the methylene groups were calculated, and their temperature dependences were studied. The differences in the temperature dependences of the energy contributions of methylene groups ofn-alkanes on Ful-60 and OV-1 were revealed. The entropy contribution is 68–82% of the enthalpy contribution which indicates a substantial role of the number of contacts with Ful-60 in retention ofn-alkanes. The ability of Ful-60 for dispersive interactions is similar to those of nonpolar liquid phases and substantially differs from that for carbon adsorbents. Fullerene columns were shown to be convenient for analysis of highly boiling organic substances in aqueous and organic solutions. Translated fromIzvestiya Akademii Nauk. Seriya Khimicheskaya, No. 8, pp. 1490–1495, August, 1999.     

Accurate ab initio potential energy curve of F2. I. Nonrelativistic full valence configuration interaction energies using the correlation energy extrapolation by intrinsic scaling method

The recently introduced method of correlation energy extrapolation by intrinsic scaling (CEEIS) is used to calculate the nonrelativistic electron correlations in the valence shell of the F(2) molecule at 13 internuclear distances along the ground state potential energy curve from 1.14 A to 8 A, the equilibrium distance being 1.412 A. Using Dunning's correlation-consistent double-, triple-, and quadruple-zeta basis sets, the full configuration interaction energies are determined, with an accuracy of about 0.3 mhartree, by successively generating up to octuple excitations with respect to multiconfigurational reference functions that strongly change along the reaction path. The energies of the reference functions and those of the correlation energies with respect to these reference functions are then extrapolated to their complete basis set limits. The applicability of the CEEIS method to strongly multiconfigurational reference functions is documented in detail.     

Temperature and pressure dependence of the reaction of OH and CO: Master equation modeling on a high‐level potential energy surface

The temperature and pressure dependence of the reaction of OH + CO has been modeled using the (energy‐resolved) master equation and RRKM theory. These calculations are based on the coupled‐cluster potential energy surface of Yu and co‐workers (Chem Phys Lett 349, 547–554, 2001). As is well known, this reaction shows a strong non‐Arrhenius behavior at moderate and low temperatures because of the stabilization of the HOCO intermediate. Kinetic simulations are in excellent agreement with experiments at temperatures above 300 K, but the agreement is only modest at temperatures below 250 K. Our calculations indicate that the contribution of tunneling to the rate constant is marginal, given the small energy difference between the transition states corresponding to formation and decomposition of the HOCO intermediate. Parametric fits to the calculated rate constants are provided for modeling purposes. © 2003 Wiley Periodicals, Inc. Int J Chem Kinet 35: 464–474, 2003