Energy decomposition in molecular complexes: implications for the treatment of polarization in molecular simulations

This study examines the contribution of electrostatic and polarization to the interaction energy in a variety of molecular complexes. The results obtained from the Kitaura-Morokuma (KM) energy decomposition analysis at the HF/6-31G(d) level indicate that, for intermolecular distances around the equilibrium geometries, the polarization energy can be determined as the addition of the polarization energies of interacting blocks, as the mixed polarization term is typically negligible. Comparison of KM and QM/MM results shows that the electrostatic energy determined in the KM method is underestimated (in absolute value) by QM/MM methods. The reason of such underestimation can be attributed to the simplified representation of treating the interaction between overlapping charge distribution by the interaction of a QM molecule with a set of point charges. Nevertheless, the polarization energies calculated by KM and QM/MM methods are in close agreement. Finally, a consistent, automated strategy to derive charge distributions that include implicitly polarization effects in pairwise, additive force fields is presented. The strategy relies in the simultaneous fitting of electrostatic and polarization energies computed by placing a suitable perturbing particle at selected points around the molecule. The suitability of these charges to describe molecular interactions is discussed.     

Quantum mechanics/molecular mechanics electrostatic embedding with continuous and discrete functions

A quantum mechanics/molecular mechanics (QM/MM) implementation that uses the Gaussian electrostatic model (GEM) as the MM force field is presented. GEM relies on the reproduction of electronic density by using auxiliary basis sets to calculate each component of the intermolecular interaction. This hybrid method has been used, along with a conventional QM/MM (point charges) method, to determine the polarization on the QM subsystem by the MM environment in QM/MM calculations on 10 individual H(2)O dimers and a Mg(2+)-H(2)O dimer. We observe that GEM gives the correct polarization response in cases when the MM fragment has a small charge, while the point charges produce significant over-polarization of the QM subsystem and in several cases present an opposite sign for the polarization contribution. In the case when a large charge is located in the MM subsystem, for example, the Mg(2+) ion, the opposite is observed at small distances. However, this is overcome by the use of a damped Hermite charge, which provides the correct polarization response.     

A minimal implementation of the AMBER-GAUSSIAN interface for ab initio QM/MM-MD simulation

For applying to a number of theoretical methodologies based on an ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics method connecting AMBER9 with GAUSSIAN03, we have developed an AMBER-GAUSSIAN interface (AG-IF), which can be one of the simplest architectures. In the AG-IF, only a few subroutines addition is necessary to retrieve the QM/MM energy and forces, obtained by GAUSSIAN, for solving a set of Newtonian equations of motion in AMBER. It is, therefore, easy to be modified for individual applications since AG-IF utilizes most of those functions originally equipped not only in AMBER but also in GAUSSIAN. In the present minimal implementation, only AMBER is modified, whereas GAUSSIAN is left unchanged. Moreover, a different method of calculating electrostatic forces of MM atoms interacting with QM region is proposed. Using the AG-IF, we also demonstrate three examples of application: (i) the QM versus MM comparison in the radial distribution function, (ii) the free energy gradient method, and (iii) the charge from interaction energy and forces.     

Calculation of the polychlorinated terphenyl congener distribution and gas chromatographic retention behavior of selected single standards on four stationary liquid phases

The distribution of the 8557 theoretically possible structures of polychlorinated terphenyls (PCT) over the different degrees of chlorination and basic terphenyl backbones has been calculated and is presented in detail. As congeners missing chlorine atoms in an ortho position are suspected to be as noxious to organisms as dioxins, the number of these congeners is specified as well. The gas chromatographic retention behavior of 14 PCT single standards is studied on four stationary liquid phases for capillary columns (HP-Ultra 1, HP-5MS, RSL-300 and Polycarbonate) relative to 2,2′,3,4,4′,5,5′-heptachlorobiphenyl and shows no consequent separation of PCT in terms of the degree of chlorination. Individual congeners have been assessed for their suitability as internal standards in PCT analysis.     

Interactions between chlorinated dioxins and a positively charged molecular probe: New molecular interaction potential

Interaction with the ligand binding domain of receptors for natural chemicals present one potential mechanism for the biological effects of environmental chemicals. Evidence suggests that the electrostatic interaction between the ligand and the receptor is an important component for binding to some of the relevant receptors. The presence of charged residues near the binding site suggests that the charge distribution of the free ligand may be different from the charge distribution of the ligand as it approaches the binding domain of the protein. In this study a new type of potential is computed for a series of dibenzo-p-dioxin (dioxin) ligands. This quantum mechanically computed potential results from interaction between the ligand and a trimethyl ammonium probe at a set of grid points. This interaction potential is compared with the molecular electrostatic potential computed from the wave function of the isolated ligands. Three types of local minima are found: (1) above the oxygen; (2) above the conjugated ring; and (3) above the chlorine(s). The molecular electrostatic potential emphasizes the minima associated with the chlorine atoms and, in that potential, the minima associated with the oxygen atoms disappear with chlorination. In the new potential, the minima over the oxygen atoms are maintained even in tetrachlorodioxin. As chlorination is increased the differences between the two potentials increases. The new potential shows the influence of the π-cation interaction, which is largest when there is little substitution on the ring. The presence of the probe induces a dipole component of 1 debye perpendicular to the plane of the ligand. Local minima in the interaction potential are then used as starting structures for the determination of the most stable ligand–probe complexes. The most stable structures are obtained from the minima associated with the oxygen atoms. These structures are stabilized by a hydrogen bond formation between the probe and the oxygen and the molecule is bent by 30° about the O(SINGLE BOND)O axis. For this series of molecules, the new potential retains some of the features that determine the hydrogen bond whereas the molecular electrostatic potential does not. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 673–684, 1998     

Isomer specific determination of polychlorinated technical mixtures such as polychloronaphthalenes by GC/MS

Summary A procedure has been developed which allows the determination of the concentration and composition of polychlorinated technical mixtures in environmental samples avoiding the use of standard compounds as much as possible. A common MS response of the congeners of each degree of chlorination is assumed and determined by comparison of the FID and MS response of the respective congeners in a technical mixture. A measure for the variability of the assumed common MS response is derived so that its validity is proved in an objective manner. The method presented here allows the isomer specific determination of those polychlorinated technical mixtures in environmental samples which could not so far be determined in this way for lack of standards. The polychlorinated naphthalenes are quoted as illustration and a few data of residues in environmental samples are given.     

Development of a quantification methodology for polychlorinated biphenyls by using Kanechlor products as the secondary reference standard

Kanechlor (KC)-300, 400, 500 and 600, Japanese polychlorinated biphenyl (PCB) products, and their equivalent mixture were analyzed by using a gas chromatograph (GC) equipped with an SE-54 capillary column/electron capture detector (ECD) and a GC/mass spectrometer in the selected ion monitoring mode (MS-SIM). All peaks were assigned to the composing congeners based on the data on peak assignment of Clophen A-30, 40, 50, 60 and Aroclor 1016, 1242, 1254, 1260 [1] and on the relative retention time values of 209 PCB congeners [2]. The weight percentage of the congener(s) which corresponds to each peak in the mass chromatograms was calculated by comparison of its height with that of certified reference standard with the same molecular weight. Each weight percentage of PCB congener(s) corresponding to each ECD peak was obtained by summing up the percent contribution values of the PCB congeners co-eluting. The results showed that it was possible to use KC products and their equivalent mixture as secondary reference standards for congener-specific PCB quantification.     

Calculation of the polychlorinated terphenyl congener distribution and gas chromatographic retention behavior of selected single standards on four stationary liquid phases

The distribution of the 8557 theoretically possible structures of polychlorinated terphenyls (PCT) over the different degrees of chlorination and basic terphenyl backbones has been calculated and is presented in detail. As congeners missing chlorine atoms in an ortho position are suspected to be as noxious to organisms as dioxins, the number of these congeners is specified as well. The gas chromatographic retention behavior of 14 PCT single standards is studied on four stationary liquid phases for capillary columns (HP-Ultra 1, HP-5MS, RSL-300 and Polycarbonate) relative to 2,2′,3,4,4′,5,5′-heptachlorobiphenyl and shows no consequent separation of PCT in terms of the degree of chlorination. Individual congeners have been assessed for their suitability as internal standards in PCT analysis. Received: 16 January 1998 / Revised: 3 April 1998 / Accepted: 7 April 1998     

Effect of polarization on the opsin shift in rhodopsins. 2. Empirical polarization models for proteins

The explicit treatment of polarization as a many-body interaction in condensed-phase systems represents a current problem in empirical force-field development. Although a variety of efficient models for molecular polarization have been suggested, polarizable force fields are still far from common use nowadays. In this work, we consider interactive polarization models employing Thole's short-range damping scheme and assess them for application on polypeptides. Despite the simplicity of the model, we find mean polarizabilities and anisotropies of amino acid side chains in excellent agreement with MP2/cc-pVQZ benchmark calculations. Combined with restrained electrostatic potential (RESP) derived atomic charges, the models are applied in a quantum-mechanical/molecular-mechanical (QM/MM) approach. An iterative scheme is used to establish a self-consistent mutual polarization between the QM and MM moieties. This ansatz is employed to study the influence of the protein polarizability on calculated optical properties of the protonated Schiff base of retinal in rhodopsin (Rh), bacterio-rhodopsin (bR), and pharaonis sensory rhodopsin II (psRII). The shifts of the excitation energy due to the instantaneous polarization response of the protein to the charge transfer on the retinal chromophore are quantified using the high level ab initio multireference spectroscopy-oriented configuration interaction (SORCI) method. The results are compared with those of previously published QM1/QM2/MM models for bR and psRII.     

Conformational dependence of charges in protein simulations

We have studied the conformational dependence of molecular mechanics atomic charges for proteins by calculating the charges fitted to the quantum mechanical (QM) electrostatic potential (ESP) for all atoms in complexes between avidin and seven biotin analogues for 20 snapshots from molecular dynamics simulations. We have studied how various other charge sets reproduce those charges. The QM charges, even if averaged over all snapshots or all residues, in general have a larger magnitude than standard Amber charges, indicating that the restraint toward zero in the restrained ESP method is too strong. This has a significant influence on the electrostatic conformational energies and the interaction energy between the biotin ligand and the protein, giving a difference between the QM and Amber charges of 43 and 8 kJ/mol for the negatively charged and neutral biotin analogues, respectively (3-4%). However, this energy difference is strongly reduced if the solvation energy (calculated by the Poisson-Boltzmann or Generalized Born methods) is added, viz., to 7 kJ/mol for charged and 3 kJ/mol for uncharged ligand. In fact, charges need to be recalculated with a QM method only for residues within 7 or 4 A of the ligand, if the error should be less than 4 kJ/mol. Unfortunately, the QM charges do not give significantly better MM/PBSA estimates of ligand-binding affinities than standard Amber charges.     

Micellar electrokinetic chromatography of polychlorinated biphenyl congeners using a polymeric surfactant as the pseudostationary phase

Micellar electrokinetic chromatography (MEKC) of highly hydrophobic compounds is generally difficult using sodium dodecyl sulfate micellar solutions. The polymeric surfactant, polysodium undecyl sulfate (poly-SUS) has been used to separate moderately to highly hydrophobic polychlorinated biphenyl (PCB) congeners by MEKC in the absence of cyclodextrins. Parameters such as concentration of acetonitrile (ACN), polymeric surfactant concentration, and the effect of pH were examined. Optimum MEKC conditions to get baseline resolution of nine PCBs was 7.5 mM borate in 40% (v/v) ACN fraction buffered at pH 9.2 using 0.5% (w/v) poly-SUS. The applied voltage was 30 kV and the temperature was maintained at 25 degrees C. Elution order for each PCB congener was found to be dependent on the degree of chlorination and hydrophobic character.     

Improved Procedure for Single DB-XLB Column GC-MS-SIM Quantitation of PCB Congener Distributions and Characterization of Two Different Preparations Sold as “Aroclor 1254”

Complete PCB congener distributions in a panel of Aroclor mixtures were previously obtained by combining data from several HRGC systems. In that study quantitation of minor components may have been unreliable due to single level calibration against high levels of individual congener standards. Two lots of Aroclor 1254 had markedly different congener distributions. In this study, the design and performance of a congener-specific PCB analysis method employing GC-MS-SIM detection of congeners separated on a DB-XLB capillary column are discussed. Quantitation is carried out against a 6-level inclusive standard curve of a mixture of 144 congeners found in Aroclors. A separate procedure to measure trace levels of PCB 126 in Aroclors using the same system, combined with levels initially acquired for other congeners, facilitates estimation of TEQ values (Toxic Equivalencies of the PCB mixtures to 2,3,7,8-TCDD). PCB congener profiles of 15 Aroclor 1254 mixtures are presented. These profiles show that the less common, high TEQ variety of Aroclor 1254 was manufactured by an atypical, two-stage chlorination process that was apparently used during the final 1% of Aroclor 1254 production (ca. 1974–1976).     

The large quadrupole of water molecules

Many quantum mechanical calculations indicate water molecules in the gas and liquid phase have much larger quadrupole moments than any of the common site models of water for computer simulations. Here, comparisons of multipoles from quantum mechanical∕molecular mechanical (QM∕MM) calculations at the MP2∕aug-cc-pVQZ level on a B3LYP∕aug-cc-pVQZ level geometry of a waterlike cluster and from various site models show that the increased square planar quadrupole can be attributed to the p-orbital character perpendicular to the molecular plane of the highest occupied molecular orbital as well as a slight shift of negative charge toward the hydrogens. The common site models do not account for the p-orbital type electron density and fitting partial charges of TIP4P- or TIP5P-type models to the QM∕MM dipole and quadrupole give unreasonable higher moments. Furthermore, six partial charge sites are necessary to account reasonably for the large quadrupole, and polarizable site models will not remedy the problem unless they account for the p-orbital in the gas phase since the QM calculations show it is present there too. On the other hand, multipole models by definition can use the correct multipoles and the electrostatic potential from the QM∕MM multipoles is much closer than that from the site models to the potential from the QM∕MM electron density. Finally, Monte Carlo simulations show that increasing the quadrupole in the soft-sticky dipole-quadrupole-octupole multipole model gives radial distribution functions that are in good agreement with experiment.     

Interactions between phosphate and water in solution: a natural bond orbital based analysis in a QM/MM framework

The natural bond orbital (NBO) and natural energy decomposition analysis (NEDA) calculations are used to analyze the interaction between mono-methyl phosphate-ester (MMP) and its solvation environment in a combined quantum mechanical/molecular mechanical (QM/MM) framework. The solute-solvent configurations are generated using a specific parametrization of the self-consistent-charge density functional tight-binding (SCC-DFTB) model for the MMP and TIP3P for water. The NBO and NEDA calculations are done with several QM/MM partitioning schemes with HF/6-31+G** as the QM level. Regardless of the size of the QM region, a notable amount of charge transfer is observed between MMP and the neighboring water molecules and the charge-transfer interactions are, in the NEDA framework, as important as the electric (electrostatic and polarization) components. This work illustrates that NBO based analyses are effective tools for probing intermolecular interactions in condensed phase systems.     

Geometry optimization based on linear response free energy with quantum mechanical/molecular mechanical method: applications to Menshutkin-type and Claisen rearrangement reactions in aqueous solution

The authors present a method based on a linear response theory that allows one to optimize the geometries of quantum mechanical/molecular mechanical (QM/MM) systems on the free energy surfaces. Two different forms of linear response free energy functionals are introduced, and electronic wave functions of the QM region, as well as the responses of electrostatic and Lennard-Jones potentials between QM and MM regions, are self-consistently determined. The covariant matrix relating the QM charge distribution to the MM response is evaluated by molecular dynamics (MD) simulation of the MM system. The free energy gradients with respect to the QM atomic coordinates are also calculated using the MD trajectory results. They apply the present method to calculate the free energy profiles of Menshutkin-type reaction of NH3 with CH3Cl and Claisen rearrangement of allyl vinyl ether in aqueous solution. For the Menshutkin reaction, the free energy profile calculated with the modified linear response free energy functional is in good agreement with that by the free energy perturbation calculations. They examine the nonequilibrium solvation effect on the transmission coefficient and the kinetic isotope effect for the Claisen rearrangement.     

Electrostatic features of molecular recognition by cyclic urea mimics of chymotrypsin

The molecular electrostatic-potential pattern was used to investigate the electrostatic features of molecular recognition by two cyclic urea mimics of the active site of α-chymotrypsin. The structures of the mimics were obtained by molecular-mechanics evaluation of the conformational potential-energy surface of the molecules. Calculations were done by using two different atomic point-charge sets in order to assess the effect of charge on the electrostatic potential pattern. The molecules studies were: (1) a “full” mimic of chymotrypsin containing the hydroxyl, imidazole, and carboxylate anion functionalities typical of the active site of the enzyme, and (2) a “partial” mimic with only the hydroxyl and imidazole functional groups. Comparison of the molecular electrostatic-potential patterns of the two mimics in both charge sets showed that the largest differences were due to the structural addition of the carboxylate anion, rather than any particular differences in the choice of atomic point charge. For the full mimic, the pattern was essentially dominated by the negative charge on the carboxylate. Small structural changes which optimized the orientation of the catalytic components had little effect on the electrostatic potential pattern of the molecule. In the absence of the anionic functionality, greater differences were noted in the electrostatic potential pattern of the partial mimic in the two charge sets. The choice of atomic point charge was seen to influence the hydrogen-bonding pattern of the hydroxyl and imidazole moieties, resulting in differences in the spatial orientation of the electrostatic potential minima. In general, both charge sets produced molecular electrostatic-potential patterns which indicated that long-range electrostatic interactions would direct the cationic end of the substrate into the electron-rich binding site. However, specific local features of the electrostatic potential pattern were found to depend on point-charge set through the influence of charge on the hydrogen-bonding pattern.     

Intraprotein electrostatics derived from first principles: divide-and-conquer approaches for QM/MM calculations

Two divide-and-conquer (DAQ) approaches for building multipole-based molecular electrostatic potentials of proteins are presented and evaluated for use in QM/MM calculations. One approach is a further development of the neutralization method of Bellido and Rullmann (J Comput Chem 1989, 10, 479-487) while the other is based on removing part of the electron density before performing the multipole expansion. Both methods create systems with integer charges without using charge renormalization. To determine their performance in terms of location of cuts and distance to QM region, the new DAQ approaches are tested in calculations of the proton affinity of N(zeta) of Lys55 in the inhibitor turkey ovomucoid third domain. Finally, the two methods are used to build a variety of MM regions, applied to calculations of the pK(a) of Lys55, and compared to other computational methodologies in which force field charges are employed.