Chemical potential perturbation: extension of the method to lattice sum treatment of intermolecular potentials

The recently developed chemical potential perturbation (CPP) method [S. G. Moore and D. R. Wheeler, J. Chem. Phys. 134, 114514 (2011)] is extended to the lattice (Ewald) sum treatment of intermolecular potentials. The CPP method predicts chemical potentials for a range of composition points using the local (position-dependent) pressure tensor of an inhomogeneous system. When computing the local pressure tensor, one can use the Irving-Kirkwood (IK) or Harasima (H) contours of distributing the pressure. We compare these two contours and show that for a planar interface, the homogeneous pressure and resulting chemical potential can be approximated with the CPP method using either the IK or the H contour, though with the lattice sum method the H contour has much greater computational efficiency. The proposed methods are validated by calculating the chemical potentials of the Lennard-Jones fluid and extended simple point-charge (SPC/E) water, and results show a high level of agreement with respective equations of state.     

Field induced gradient simulations: a high throughput method for computing chemical potentials in multicomponent systems

We present a simulation method for direct computation of chemical potentials in multicomponent systems. The method involves application of a field to generate spatial gradients in the species number densities at equilibrium, from which the chemical potential of each species is theoretically estimated. A single simulation yields results over a range of thermodynamic states, as in high throughput experiments, and the method remains computationally efficient even at high number densities since it does not involve particle insertion at high densities. We illustrate the method by Monte Carlo simulations of binary hard sphere mixtures of particles with different sizes in a gravitational field. The results of the gradient Monte Carlo method are found to be in good agreement with chemical potentials computed using the classical Widom particle insertion method for spatially uniform systems.     

Analytic gradient and molecular dynamics simulations using the fragment molecular orbital method combined with effective potentials

The completely analytic energy gradients are derived and implemented for the two-body fragment molecular orbital (FMO2) method combined with the model core potentials (MCP) and effective fragment potentials (EFP). The many-body terms in EFP require solving coupled-perturbed Hartree-Fock equations, which are derived and implemented. The molecular dynamics (MD) simulations are performed using the FMO2/MCP method for the capped alanine decamer and with the FMO2/EFP method for the zwitterionic conformer of glycine tetramer immersed in the water layer of 6.0 Å (135 water molecules). The results of the MD simulations using the FMO2/EFP and FMO2/MCP gradients show that the total energy is conserved at the time steps less than 1 fs.     

A density-functional approach to polarizable models: a Kim-Gordon response density interaction potential for molecular simulations

A combined linear-response-frozen electron-density model has been implemented in a molecular-dynamics scheme derived from an extended Lagrangian formalism. This approach is based on a partition of the electronic charge distribution into a frozen region described by Kim-Gordon theory [J. Chem. Phys. 56, 3122 (1972); J. Chem. Phys. 60, 1842 (1974)] and a response contribution determined by the instantaneous ionic configuration of the system. The method is free from empirical pair potentials and the parametrization protocol involves only calculations on properly chosen subsystems. We apply this method to a series of alkali halides in different physical phases and are able to reproduce experimental structural and thermodynamic properties with an accuracy comparable to Kohn-Sham density-functional calculations.     

An efficient perturbation method to predict the functionally key sites of glutamine binding protein

Glutamine-Binding Protein (GlnBP) of Escherichia coli, an important member of the periplasmic binding protein family, is responsible for the first step in the active transport of glutamine across the cytoplasmic membrane. In this work, the functionally key regulation sites of GlnBP were identified by utilizing a perturbation method proposed by our group, in which the residues whose perturbations markedly change the binding free energy between GlnBP and glutamine are considered to be functionally key residues. The results show that besides the substrate binding sites, some other residues distant from the binding pocket, including the ones in the hinge regions between the two domains, the front- and back- door channels and the exposed region, are important for the function of glutamine binding and transport. The predicted results are well consistent with the theoretical and experimental data, which indicates that our method is an effective approach to identify the key residues important for both ligand binding and long-range allosteric signal transmission. This work can provide some insights into the function performance of GlnBP and the physical mechanism of its allosteric regulation.     

Application of MDGRAPE-3, a special purpose board for molecular dynamics simulations, to periodic biomolecular systems

We describe the application of a special purpose board for molecular dynamics simulations, named MDGRAPE-3, to the problem of simulating periodic bio-molecular systems. MDGRAPE-3 is the latest board in a series of hardware accelerators designed to calculate the nonbonding long-range interactions much more rapidly than normal processors. So far, MDGRAPEs were mainly applied to isolated systems, where very many nonbonded interactions were calculated without any distance cutoff. However, in order to regulate the density and pressure during simulations of membrane embedded protein systems, one has to evaluate interactions under periodic boundary conditions. For this purpose, we implemented the Particle-Mesh Ewald (PME) method, and its approximation with distance cutoffs and charge neutrality as proposed by Wolf et al., using MDGRAPE-3. When the two methods were applied to simulations of two periodic biomolecular systems, a single MDGRAPE-3 achieved 30-40 times faster computation times than a single conventional processor did in the both cases. Both methods are shown to have the same molecular structures and dynamics of the systems.     

An automated method to find transition states using chemical dynamics simulations

A procedure to automatically find the transition states (TSs) of a molecular system (MS) is proposed. It has two components: high‐energy chemical dynamics simulations (CDS), and an algorithm that analyzes the geometries along the trajectories to find reactive pathways. Two levels of electronic structure calculations are involved: a low level (LL) is used to integrate the trajectories and also to optimize the TSs, and a higher level (HL) is used to reoptimize the structures. The method has been tested in three MSs: formaldehyde, formic acid (FA), and vinyl cyanide (VC), using MOPAC2012 and Gaussian09 to run the LL and HL calculations, respectively. Both the efficacy and efficiency of the method are very good, with around 15 TS structures optimized every 10 trajectories, which gives a total of 7, 12, and 83 TSs for formaldehyde, FA, and VC, respectively. The use of CDS makes it a powerful tool to unveil possible nonstatistical behavior of the system under study. © 2014 Wiley Periodicals, Inc.     

A simple method to predict protein flexibility using secondary chemical shifts

Protein motions play a critical role in many biological processes, such as enzyme catalysis, allosteric regulation, antigen-antibody interactions, and protein-DNA binding. NMR spectroscopy occupies a unique place among methods for investigating protein dynamics due to its ability to provide site-specific information about protein motions over a large range of time scales. However, most NMR methods require a detailed knowledge of the 3D structure and/or the collection of additional experimental data (NOEs, T1, T2, etc.) to accurately measure protein dynamics. Here we present a simple method based on chemical shift data that allows accurate, quantitative, site-specific mapping of protein backbone mobility without the need of a three-dimensional structure or the collection and analysis of NMR relaxation data. Further, we show that this chemical shift method is able to quantitatively predict per-residue RMSD values (from both MD simulations and NMR structural ensembles) as well as model-free backbone order parameters.     

Alamethicin channels in a membrane: molecular dynamics simulations

Alamethicin (Alm) is a 20 residue peptide which forms a kinked alpha-helix in membrane and membrane-mimetic environments. Ion channels formed by intramembraneous aggregates of Alm are thought to be formed by bundles of approximately parallel Alm helices surrounding a central bilayer pore. Different channel conductance levels correspond to different numbers of helices per bundle, ranging from N = 5 to N > 8. Calculation of the predicted pKA values of the ring of Glu18 sidechains at the C-terminal mouth of the pore suggests that at neutral pH most or all of these sidechains will remain protonated. Nanosecond molecular dynamics (MD) simulations of N = 5, 6, 7 and 8 bundles of Alm helices in a POPC bilayer have been run, corresponding to a total simulation time of 4 ns. These simulations explore the stability and conformational dynamics of these helix bundle channels when embedded in a full phospholipid bilayer in an aqueous environment. The structural and dynamic properties of water in these model channels are examined. As in earlier in vacuo simulations (J. Breed, R. Sankararamakrishnan, I. D. Kerr and M. S. P. Sansom, Biophys. J., 1996, 70, 1643) the dipole moments of water molecules within the pores are aligned antiparallel to the helix dipoles. This helps to contribute to the stability of the helix bundles.     

A new perturbation theory for potentials with a soft core: Application to liquid sodium

An extension of the Weeks—Chandler—Andersen perturbation theory, proposed by Wright and Perram is applied to liquid Na By neglecting the potential beyond its first minimum, the structure factor is found to be given accurately except in the vicinity of the first peak.     

A density functional theory-based chemical potential equalisation approach to molecular polarizability

The electron density changes in molecular systems in the presence of external electric fields are modeled for simplicity in terms of the induced charges and dipole moments at the individual atomic sites. A chemical potential equalisation scheme is proposed for the calculation of these quantities and hence the dipole polarizability within the framework of density functional theory based linear response theory. The resulting polarizability is expressed in terms of the contributions from individual atoms in the molecule. A few illustrative numerical calculations are shown to predict the molecular polarizabilities in good agreement with available results. The usefulness of the approach to the calculation of intermolecular interaction needed for computer simulation is highlighted.     

Classifying large chemical data sets: using a regularized potential function method

In recent years classifiers generated with kernel-based methods, such as support vector machines (SVM), Gaussian processes (GP), regularization networks (RN), and binary kernel discrimination (BKD) have been very popular in chemoinformatics data analysis. Aizerman et al. were the first to introduce the notion of employing kernel-based classifiers in the area of pattern recognition. Their original scheme, which they termed the potential function method (PFM), can basically be viewed as a kernel-based perceptron procedure and arguably subsumes the modern kernel-based algorithms. PFM can be computationally much cheaper than modern kernel-based classifiers; furthermore, PFM is far simpler conceptually and easier to implement than the SVM, GP, and RN algorithms. Unfortunately, unlike, e.g., SVM, GP, and RN, PFM is not endowed with both theoretical guarantees and practical strategies to safeguard it against generating overfitting classifiers. This is, in our opinion, the reason why this simple and elegant method has not been taken up in chemoinformatics. In this paper we empirically address this drawback: while maintaining its simplicity, we demonstrate that PFM combined with a simple regularization scheme may yield binary classifiers that can be, in practice, as efficient as classifiers obtained by employing state-of-the-art kernel-based methods. Using a realistic classification example, the augmented PFM was used to generate binary classifiers. Using a large chemical data set, the generalization ability of PFM classifiers were then compared with the prediction power of Laplacian-modified naive Bayesian (LmNB), Winnow (WN), and SVM classifiers.     

Fluctuating charge normal modes: An algorithm for implementing molecular dynamics simulations with polarizable potentials

A new method for performing molecular dynamics simulations with fluctuating charge polarizable potentials is introduced. In fluctuating charge models, polarizability is treated by allowing the partial charges to be variables, with values that are coupled to charges on the same molecule as well as those on other molecules. The charges can be efficiently propagated in a molecular dynamics simulation using extended Lagrangian dynamics. By making a coordinate change from the charge variables to a set of normal mode charge coordinates for each molecule, a new method is constructed in which the normal mode charge variables uncouple from those on the same molecule. The method is applied to the TIP4P-FQ model of water and compared to other methods for implementing the dynamics. The methods are compared using different molecular dynamics time steps.     

NMR chemical shifts as a tool to analyze first principles molecular dynamics simulations in condensed phases: the case of liquid water

We present (1)H NMR chemical shift calculations of liquid water based on first principles molecular dynamics simulations under periodic boundary conditions. We focus on the impact of computational parameters on the structural and spectroscopic data, which is an important question for understanding how sensitive the computed (1)H NMR resonances are upon variation of the simulation setup. In particular, we discuss the influence of the exchange-correlation functional and the size of the basis set, the choice for the fictitious electronic mass and the use of pseudopotentials for the nuclear magnetic resonance (NMR) calculation on one hand and the underlying Car-Parrinello-type molecular dynamics simulations on the other hand. Our findings show that the direct effect of these parameters on (1)H shifts is not big, whereas the indirect dependence via the structural data is more important. The (1)H NMR chemical shifts clearly reflect the induced structural changes, illustrating once again the sensitivity of (1)H NMR observables on small changes in the local chemical structure of complex hydrogen-bonded liquids.     

RRS-PBC: a molecular approach for periodic systems

Technically, when dealing with a perfect crystal, methods in k-(reciprocal) space that impose periodic boundary conditions(PBC) in conjunction with plane-wave basis sets are widely used. Chemists, however, tend to think of a solid as a giant molecule, which offers a molecular way to describe a solid by using a finite cluster model(FCM). However, FCM may fail to simulate a perfect crystal due to its inevitable boundary effects. We propose an RRS-PBC method that extracts the k-space information of a perfect crystalline solid out of a reduced real space(RRS) of an FCM. We show that the inevitable boundary effects in an FCM are eliminated naturally to achieve converged high-quality band structures.     

A recent development in computational chemistry: chemical reactions from first principles molecular dynamics simulations

First-principles molecular dynamics simulations have recently been found an effective tool to study a large variety of chemical problems. Finite temperature simulations reveal unique information, including explicit dynamical effects and the evaluation of proper free energy differences. Moreover, dynamics simulations reveal information on the flexibility of molecular systems, and elucidate, often otherwise inaccessible, mechanistic details of chemical reactions. In addition this methodology allows the study of larger, periodic, systems, revealing computationally unique information which may be directly compared to experiments on realistic chemical systems. A variety of examples will be given, although most focus on the important field of catalysis.     

Spin-labeled alkylphospholipids in a dipalmitoylphosphatidylcholine bilayer: molecular dynamics simulations

Molecular dynamics simulation has been performed to investigate the structural properties of perifosine and its synthetic spin-labeled alkylphospholipid analogues. The conformations adopted by these compounds in water and in a dipalmitoylphosphatidylcholine bilayer as a function of the presence and position of the N-oxyl-4',4'-dimethyloxazolidine ring (doxyl group) have been investigated by all-atom molecular dynamics. No predominant conformation was observed in water, but the molecules adopt specific orientations and conformations in the lipid bilayer. As is expected, alkyl chains tend to insert into the hydrophobic core, while charged groups stay at the lipid-water interface. A doxyl group in the middle of the alkyl chain moves up to the interface region, thus preventing adoption of the extended conformation. Compounds with a doxyl group close to the polar head group adopt conformations similar to that of unlabeled perifosine within the first nanoseconds of simulation. When the doxyl group is at the end of alkyl chain, the spin-labeled molecule needs more time to reach equilibrium. These results indicate a considerable effect of the doxyl position within the alkyl chain on its localization in the lipid bilayer and can be extended further to other similar spin probes used in the electron paramagnetic resonance spectroscopy of biological membranes.