Comparing search strategies for finding global optima on energy landscapes

We provide some tests of the convex global underestimator (CGU) algorithm, which aims to find global minima on funnel-shaped energy landscapes. We use two different potential functions—the reduced Lennard–Jones cluster potential, and the modified Sun protein folding potential, to compare the CGU algorithm with the simplest versions of the traditional trajectory-based search methods, simulated annealing (SA), and Monte Carlo (MC). For both potentials, the CGU reaches energies lower on the landscapes than both SA and MC, even when SA and MC are given the same number of starting points as in a full CGU run or when all methods are given the same amount of computer time. The CGU consistently finds the global minima of the Lennard–Jones potential for all cases with up to at least n=30 degrees of freedom. Finding the global or near-global minimum in the CGU method requires polynomial time [scaling between O(n³) and O(n⁴)], on average. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1527–1532, 1999     

Free energy of solvation of a small Lennard–Jones particle

In molecular dynamics (MD) and Monte Carlo (MC) free energy calculations, the choices of the thermodynamic paths from state a to state b affect the accuracy of the result and the efficiency of the programs. Most of the problems occur at the initial stages of growing in a new particle into a solvent. Based on statistical mechanical perturbation theory, an accurate and efficient direct calculation of inserting a small Lennard–Jones particle into solvent is derived. This eliminates the need for calculation of the initial stages of growing in a new particle by MD or MC simulation. Examples are given to show the utility of direct calculation. The recommended procedure is to use direct calculation for a small Lennard–Jones particle and then use MD or MC simulations to calculate the ΔG of changing the small Lennard–Jones particle into the target molecule. © 1994 by John Wiley & Sons, Inc.     

A Molecular Model of the Elastic Behavior of Polymethylene Chains

In this paper, the elastic behavior of a polymethylene (PM) chain is investigated by using a realistic rotational‐isomeric‐state (RIS) model. In our calculation, the non‐local interactions between pairs of segments in a polymethylene chain are also considered, and the Lennard‐Jones (L‐J) potential function is adopted. Chain dimensions and thermodynamics statistical properties of PM chains with various elongation ratios λ are calculated. We find, that the elastic force increases slowly with elongation ratio for small λ, and abruptly for large λ. In the meantime, the energy contribution to elastic force is negative and significant, especially for large λ. Our calculations may provide some insight into the macroscopic phenomena of rubber elasticity.     

Viscosity Prediction for Oxygen,Nitrogen and Their Mixtures at Zero and Moderately Dense Regimes via Semi‐Empirically Based Assessment

The viscosity coefficients for the gaseous states of N₂ and O₂ and their mixtures are determined at zero and moderately density regimes. The Lennard‐Jones 12–6 (LJ 12–6) potential energy function is used as the initial model potential required y the technique. The interaction potential energies from the inversion procedure reproduce the viscosity commensurate to the best measurements. The initial density dependence of gaseous viscosity coefficient according to the Rainwater‐Friend theory, which was given by Najafi et al., has been considered for pure N₂ and pure O₂.     

Exponential repulsion improves structural predictability of molecular docking

Molecular docking is a powerful tool for theoretical prediction of the preferred conformation and orientation of small molecules within protein active sites. The obtained poses can be used for estimation of binding energies, which indicate the inhibition effect of designed inhibitors, and therefore might be used for in silico drug design. However, the evaluation of ligand binding affinity critically depends on successful prediction of the native binding mode. Contemporary docking methods are often based on scoring functions derived from molecular mechanical potentials. In such potentials, nonbonded interactions are typically represented by electrostatic interactions between atom‐centered partial charges and standard 6–12 Lennard–Jones potential. Here, we present implementation and testing of a scoring function based on more physically justified exponential repulsion instead of the standard Lennard–Jones potential. We found that this scoring function significantly improved prediction of the native binding modes in proteins bearing narrow active sites such as serine proteases and kinases. © 2016 Wiley Periodicals, Inc.     

Effects of the Bead‐Solvent Interaction on the Dynamics of Macromolecules, 1

Summary: Hamiltonian dynamics and a chain model are used to study the dynamics of macromolecules immersed in a solution. From the Hamiltonian of the overall system, “macromolecule + solvent,” a master and a Fokker‐Planck equation are then derived for the phase‐space distribution of the macromolecule. In the Fokker‐Planck equation, all the information about the interaction among the beads of the macromolecule as well as the effects of the surrounding solvent is described by friction tensors, which are expressed in terms of the bead‐solvent interaction and the dynamic structure factor of the solvent. To explore the influence of the bead‐solvent potential on the dynamics of macromolecules, the friction tensors are calculated for a dumbbell molecule and for three choices of the interaction (Yukawa, Born‐Mayer, and Lennard‐Jones). Expressions are derived, in particular, for the friction tensor coefficients of the center‐of‐mass and the relative coordinates of the dumbbell. For the long‐term behaviour of the internal momentum autocorrelation function, moreover, an “algebraic decay” is found, in contrast to the (unphysical) exponential decay as known from phenomenological theory.

Yukawa, Born‐Mayer and Lennard‐Jones bead‐solvent interaction potentials.     

The SAAP force field: Development of the single amino acid potentials for 20 proteinogenic amino acids and Monte Carlo molecular simulation for short peptides

Molecular simulation by using force field parameters has been widely applied in the fields of peptide and protein research for various purposes. We recently proposed a new all‐atom protein force field, called the SAAP force field, which utilizes single amino acid potentials (SAAPs) as the fundamental elements. In this article, whole sets of the SAAP force field parameters in vacuo, in ether, and in water have been developed by ab initio calculation for all 20 proteinogenic amino acids and applied to Monte Carlo molecular simulation for two short peptides. The side‐chain separation approximation method was employed to obtain the SAAP parameters for the amino acids with a long side chain. Monte Carlo simulation for Met‐enkephalin (CHO‐Tyr‐Gly‐Gly‐Phe‐Met‐NH₂) by using the SAAP force field revealed that the conformation in vacuo is mainly controlled by strong electrostatic interactions between the amino acid residues, while the SAAPs and the interamino acid Lennard‐Jones potentials are predominant in water. In ether, the conformation would be determined by the combination of the three components. On the other hand, the SAAP simulation for chignolin (H‐Gly‐Tyr‐Asp‐Pro‐Glu‐Thr‐Gly‐Thr‐Trp‐Gly‐OH) reasonably reproduced a native‐like β‐hairpin structure in water although the C‐terminal and side‐chain conformations were different from the native ones. It was suggested that the SAAP force field is a useful tool for analyzing conformations of polypeptides in terms of intrinsic conformational propensities of the single amino acid units. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Orientation of Polymer Chains in Dilute Solution under Shear: Effect of Chain Model and Excluded Volume

Summary: Polymer orientation in dilute solutions undergoing shear flow is investigated computationally by means of the Brownian dynamics simulation technique applied to the bead‐spring chain model. The dependence of the degree of orientation on the shear intensity is evaluated through a quantity called orientation resistance. All simulations were performed using non‐preaveraged hydrodynamic interaction (HI). The spring type (Gaussian or FENE) is shown to strongly determine the shear flow behavior of the chain orientation. Solvent quality (Θ, good or bad), represented by a suitable Lennard‐Jones intramolecular potential, does not affect the flow behavior but influences the values of the orientation resistance. Hence, the orientability of the polymer molecule is, in a way, related to the flow intensity.

Evolution of m_G (orientational resistance parameter, open circles are simulation, dashed line is Gaussian approximation) and m_τ (filled circles are simulation, dotted line is Gaussian approximation) with β for ideal Gaussian chains with N = 15.     

Evaluating thermodynamic integration performance of the new amber molecular dynamics package and assess potential halogen bonds of enoyl‐ACP reductase (FabI) benzimidazole inhibitors

Thermodynamic integration (TI) can provide accurate binding free energy insights in a lead optimization program, but its high computational expense has limited its usage. In the effort of developing an efficient and accurate TI protocol for FabI inhibitors lead optimization program, we carefully compared TI with different Amber molecular dynamics (MD) engines (sander and pmemd), MD simulation lengths, the number of intermediate states and transformation steps, and the Lennard‐Jones and Coulomb Softcore potentials parameters in the one‐step TI, using eleven benzimidazole inhibitors in complex with Francisella tularensis enoyl acyl reductase (FtFabI). To our knowledge, this is the first study to extensively test the new AMBER MD engine, pmemd, on TI and compare the parameters of the Softcore potentials in the one‐step TI in a protein‐ligand binding system. The best performing model, the one‐step pmemd TI, using 6 intermediate states and 1 ns MD simulations, provides better agreement with experimental results (RMSD = 0.52 kcal/mol) than the best performing implicit solvent method, QM/MM‐GBSA from our previous study (RMSD = 3.00 kcal/mol), while maintaining similar efficiency. Briefly, we show the optimized TI protocol to be highly accurate and affordable for the FtFabI system. This approach can be implemented in a larger scale benzimidazole scaffold lead optimization against FtFabI. Lastly, the TI results here also provide structure‐activity relationship insights, and suggest the parahalogen in benzimidazole compounds might form a weak halogen bond with FabI, which is a well‐known halogen bond favoring enzyme. © 2015 Wiley Periodicals, Inc.     

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.     

Free volume properties of model fluids and polymers: Shape and connectivity

The Cavity Energetic Sizing Algorithm (CESA) method of in 't Veld (J Phys Chem B 2000, 104, 12028) is extended to characterize the nonspherical nature of free volume. The new technique is introduced with reference to simple model fluids (water, hard spheres, and a Lennard‐Jones fluid) and then applied to polymers of interest to membrane scientists. A set of shape parameters is introduced, characterizing nanopores in terms of surface area, volume, radius of gyration, and span. Results are presented for Lennard‐Jones fluid and hard sphere fluid, and for the high free volume polymers (poly‐trimethyl‐silyl‐propane) poly(1‐trimethylsilyl‐1‐propyne) (PTMSP) and a random copolymer of 2,2‐bis(trifluoromethyl)‐4,5‐difluoro‐1,3‐dioxole (TFE/BDD). PTMSP is observed to have an average free volume cluster span of 1.43 nm, compared to TFE/BDD with an average cluster span of 0.98 nm, consistent with the markedly higher permeability of CO₂ observed in PTMSP. An additional method for measuring free volume is introduced, similar to a method introduced by Greenfield and Theodorou (Macromolecules 1993, 26, 5461; Mol Simul 1997, 19, 329; Macromolecules 1998, 31, 7068; 2001, 34, 8541), which measures free volume relative to a specific probe. The method captures 1–3 times the fractional cavity volume captured by CESA. Free volume measurements are presented for a set of polysulfones with respect to noble gas probes (J Chem Phys 2005, 122, 84906; J Mol Struct 2005, 739, 173). © 2006 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys 44:1385–1393, 2006     

Second virial coefficients of Lennard–Jones chains

The second virial coefficients B₂ of Lennard–Jones chain fluids were calculated through Monte Carlo integration as a function of chain length m (up to 48 segments) and temperature. We found that at a fixed temperature the second virial coefficient decreases with chain length. At low temperatures, the virial coefficient changes sign from positive to negative as m increases. The simulation data also provide an estimate for the theta temperature T_Θ at which the attractive and repulsive interactions cancel each other for dilute solutions. It is found that the theta temperature T_Θ for Lennard–Jones chains with m>32 is 4.65 independent of chain length m. A comparison of simulated values of B₂ with those evaluated from two different perturbation theories for chain fluid shows that these approximate theories underestimate the second virial coefficients of Lennard–Jones chains.     

Reparameterizations of the χ Torsion and Lennard‐Jones σ Parameters Improve the Conformational Characteristics of Modified Uridines

The currently available force field parameters for modified RNA residues in AMBER show significant deviations in conformational properties from experimental observations. The examination of the transferability of the recently revised torsion parameters revealed that there was an overall improvement in the conformational properties for some of the modifications but the improvements were still insufficient in describing the sugar pucker preferences (J. Chem. Inf. Model. 2014, 54, 1129–1142). Here, we report an approach for the development and fine tuning of the AMBER force field parameters for 2‐thiouridine, 4‐thiouridine, and pseudouridine with diverse conformational preferences. The χ torsion parameters were reparameterized at the individual nucleoside level. The effect of combining the revised γ torsion parameter and modifying the Lennard‐Jones σ parameters were also tested by directly comparing the conformational preferences obtained from our extensive molecular dynamics simulations with those from experimental observations. © 2016 Wiley Periodicals, Inc.     

Systematic calculation of molecular vibrational spectra through a complete Morse expansion

We propose an accurate and efficient method to compute vibrational spectra of molecules, based on exact diagonalization of an algebraically calculated matrix based on powers of Morse coordinate. The present work focuses on the 1D potential of diatomic molecules: as typical examples, we apply this method to the standard Lennard‐Jones oscillator, and to the ab initio potential of the H₂ molecule. Global cm^?1 accuracy is exhibited through the H₂ spectrum, obtained through the diagonalization of a 30 × 30 matrix. This theory is at the root of a new method to obtain globally accurate vibrational spectral data in the context of the multi‐dimensional potential of polyatomic molecules, at an affordable computational cost. © 2006 Wiley Periodicals, Inc. Int J Quantum Chem, 2007     

Polarizable empirical force field for nitrogen‐containing heteroaromatic compounds based on the classical Drude oscillator

The polarizable empirical CHARMM force field based on the classical Drude oscillator has been extended to the nitrogen‐containing heteroaromatic compounds pyridine, pyrimidine, pyrrole, imidazole, indole, and purine. Initial parameters for the six‐membered rings were based on benzene with nonbond parameter optimization focused on the nitrogen atoms and adjacent carbons and attached hydrogens. In the case of five‐member rings, parameters were first developed for imidazole and transferred to pyrrole. Optimization of all parameters was performed against an extensive set of quantum mechanical and experimental data. Ab initio data were used for the determination of initial electrostatic parameters, the vibrational analysis, and in the optimization of the relative magnitudes of the Lennard‐Jones (LJ) parameters, through computations of the interactions of dimers of model compounds, model compound‐water interactions, and interactions of rare gases with model compounds. The absolute values of the LJ parameters were determined targeting experimental heats of vaporization, molecular volumes, heats of sublimation, crystal lattice parameters, and free energies of hydration. Final scaling of the polarizabilities from the gas‐phase values by 0.85 was determined by reproduction of the dielectric constants of pyridine and pyrrole. The developed parameter set was extensively validated against additional experimental data such as diffusion constants, heat capacities, and isothermal compressibilities, including data as a function of temperature. © 2008 Wiley Periodicals, Inc. J Comput Chem 2009     

PVT properties of linear and dendritic polymers

The aim of this article is to examine the limits of applicability of the Simha‐Somcynsky (S‐S) equation of state (EOS) by comparing the pressure‐volume‐temperature (PVT) data and the derivatives (compressibility, κ, and thermal expansion coefficient, α) of anionic linear polystyrene (PS) with poly(benzyl ether) dendrimers (PBED). Fitting the PVT of PBED data to the S‐S EOS was similarly satisfactory as that of PS and the computed Lennard‐Jones (L‐J) interaction parameters showed similar errors of ca. 1%. Next, the experimental derivatives, α and κ of PS and PBED were compared with these functions computed from the S‐S EOS—good agreement was obtained for α at ambient pressure, P, indicating validity of the S‐S theory at least up to the first derivative. While the predicted κ = κ(P) dependence for PS and a linear PBED homologue was correct, for dendrimers the compressibility was higher at low pressure and it was lower at high P than theory predicts. Also the extracted values of the L‐J repulsion volume, v*, between a segment pair was smaller than expected. The specific architecture of dendrimer molecules is responsible for this behavior, since their 3D configuration is significantly different from the S‐S model with uniform segmental density and oxygen bonds in the main and side chains add flexibility. © 2009 NRC Canada. J Polym Sci Part B: Polym Phys 48: 322–332, 2010     

Equilibrium sampling of self-associating polymer solutions: a parallel selective tempering approach

We present a novel simulation algorithm based on tempering a fraction of relaxation-limiting interactions to accelerate the process of obtaining uncorrelated equilibrium configurations of self-associating polymer solutions. This approach consists of tempering (turning off) the attractive interactions for a fraction of self-associating groups determined by a biasing field h. A number of independent configurations (replicas) with overlapping Hamiltonian distributions in the expanded (NVTh) ensemble with constant NVT but different biasing fields, forming a chain of Hamiltonians, were simulated in parallel with occasional attempts to exchange the replicas associated with adjacent fields. Each field had an associated distribution of tempered interactions, average fraction of tempered interactions, and structural decorrelation time. Tempering parameters (number of replicas, fields, and exchange frequencies) were chosen to obtain the highest efficiency in sampling equilibrium configurations of a self-association polymer solution based on short serial simulation runs and a statistical model. Depending on the strength of the relaxation-limiting interactions, system size, and thermodynamic conditions, the algorithm can be orders of magnitude more efficient than conventional canonical simulation and is superior to conventional temperature parallel tempering.     

Description of self-diffusion coefficients of gases,liquids and fluids at high pressure based on molecular simulation data

The purpose of this work is to evaluate the potential of modeling the self-diffusion coefficient (SDC) of real fluids in all fluid states based on Lennard–Jones analytical relationships involving the SDC, the temperature, the density and the pressure. For that, we generated an equation of state (EOS) that interrelates the self-diffusion coefficient, the temperature and the density of the Lennard–Jones (LJ) fluid. We fit the parameters of such LJ–SDC–EOS using recent wide ranging molecular simulation data for the LJ fluid. We also used in this work a LJ pressure–density–temperature EOS that we combined with the LJ–SDC–EOS to make possible the calculation of LJ–SDC values from given temperature and pressure. Both EOSs are written in terms of LJ dimensionless variables, which are defined in terms of the LJ parameters ɛ and σ. These parameters are meaningful at molecular level. By combining both EOSs, we generated LJ corresponding states charts which make possible to conclude that the LJ fluid captures the observed behavioral patterns of the self-diffusion coefficient of real fluids over a wide range of conditions. In this work, we also performed predictions of the SDC of real fluids in all fluid states. For that, we assumed that a given real fluid behaves as a Lennard–Jones fluid which exactly matches the experimental critical temperature T_c and the experimental critical pressure P_c of the real fluid. Such an assumption implies average true prediction errors of the order of 10% for vapors, light supercritical fluids, some dense supercritical fluids and some liquids. These results make possible to conclude that it is worthwhile to use the LJ fluid reference as a basis to model the self-diffusion coefficient of real fluids, over a wide range of conditions, without resorting to non-LJ correlations for the density–temperature–pressure relationship. The database considered here contains more than 1000 experimental data points. The database practical reduced temperature range is from 0.53 to 2.4, and the practical reduced pressure range is from 0 to 68.4.     

Characterization of polyolefin melts using the polymer reference interaction site model integral equation theory with a single‐site united atom model

A general approach, based on the polymer reference interaction site model (PRISM) integral equation theory, suitable for characterizing arbitrarily complex polyolefin melts is described. We tested the method by calculating the melt structures of linear polyethylene (PE) and isotactic polypropylene (iPP) and the spinodal decomposition temperatures for PE/iPP blends. The computational expense of the PRISM calculation was reduced with a single‐site united atom model in which the polyolefin CH, CH₂, and CH₃ groups were approximated as chemically equivalent sites with spherically symmetric energetic interactions. The site–site interactions were defined by a potential function comprising a hard core with an attractive Lennard–Jones term. These energetic parameters were optimized with a central composite design strategy that enabled a simultaneous fit of experimental melt density and structure factor data. Values were obtained for PE and iPP individually and for common universal parameters that could potentially be used for all polyolefins. The rotational isomeric state–metropolis Monte Carlo (RMMC) technique was used to generate sets of conformers at specified temperatures covering the melt‐temperature range of the polymers. The characteristic ratio was used to assess the quality of the conformers and the RMMC method. Values of 9.68 for PE and 9.27 for iPP were obtained. The single‐chain structure factors calculated by the RMMC method were used to calculate the total structure factor for each melt. These were validated against published X‐ray diffraction results. © 2001 John Wiley & Sons, Inc. J Polym Sci Part B: Polym Phys 39: 1803–1814, 2001