Molecular docking is a computational approach for predicting the most probable position of ligands in the binding sites of macromolecules and constitutes the cornerstone of structure-based computer-aided drug design. Here, we present a new algorithm called Attracting Cavities that allows molecular docking to be performed by simple energy minimizations only. The approach consists in transiently replacing the rough potential energy hypersurface of the protein by a smooth attracting potential driving the ligands into protein cavities. The actual protein energy landscape is reintroduced in a second step to refine the ligand position. The scoring function of Attracting Cavities is based on the CHARMM force field and the FACTS solvation model. The approach was tested on the 85 experimental ligand–protein structures included in the Astex diverse set and achieved a success rate of 80% in reproducing the experimental binding mode starting from a completely randomized ligand conformer. The algorithm thus compares favorably with current state-of-the-art docking programs. © 2015 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

A new algorithm called Attracting Cavities, which allows molecular docking to be performed by only simple energy minimizations, is presented. The approach consists in transiently replacing the rough protein potential energy landscape by a smooth attracting landscape driving the ligands into protein cavities. The actual protein energy landscape is ultimately reintroduced to refine the ligand pose. The approach achieved a success rate of 80% on the Astex diverse set starting from completely randomized ligand conformers.

The extensibility of force field is a key to solve the missing parameter problem commonly found in force field applications. The extensibility of conventional force fields is traditionally managed in the parameterization procedure, which becomes impractical as the coverage of the force field increases above a threshold. A hierarchical atom-type definition (HAD) scheme is proposed to make extensible atom type definitions, which ensures that the force field developed based on the definitions are extensible. To demonstrate how HAD works and to prepare a foundation for future developments, two general force fields based on AMBER and DFF functional forms are parameterized for common organic molecules. The force field parameters are derived from the same set of quantum mechanical data and experimental liquid data using an automated parameterization tool, and validated by calculating molecular and liquid properties. The hydration free energies are calculated successfully by introducing a polarization scaling factor to the dispersion term between the solvent and solute molecules. © 2015 Wiley Periodicals, Inc.

A hierarchical definition of atom type is proposed to solve the problem of “missing parameters”. The extensible atom type definitions make the force field extensible. Two general force fields are parameterized for some common organic molecules. Parameters are derived from ab initial data and experimental liquid data. Calculation results show good accuracy of the parameters.

This article describes a revised version 56A6CARBO_R of the GROMOS 56A6CARBO force field for hexopyranose-based carbohydrates. The simulated properties of unfunctionalized hexopyranoses are unaltered with respect to 56A6CARBO. In the context of both O1-alkylated hexopyranoses and oligosaccharides, the revision stabilizes the regular 4C1 chair for α-anomers, with the opposite effect for β-anomers. As a result, spurious ring inversions observed in α(14)-linked chains when using the original 56A6CARBO force field are alleviated. The 4C1 chair is now the most stable conformation for all d-hexopyranose residues, irrespective of the linkage type and anomery, and of the position of the residue along the chain. The methylation of a d-hexopyranose leads to a systematic shift in the ring-inversion free energy (4C1 to 1C4) by 7–8 kJ mol−1, positive for the α-anomers and negative for the β-anomers, which is qualitatively compatible with the expected enhancement of the anomeric effect upon methylation at O1. The ring-inversion free energies for residues within chains are typically smaller in magnitude compared to those of the monomers, and correlate rather poorly with the latter. This suggests that the crowding of ring substituents upon chain formation alters the ring flexibility in a nonsystematic fashion. In general, the description of carbohydrate chains afforded by 56A6CARBO_R suggests a significant extent of ring flexibility, i.e., small but often non-negligible equilibrium populations of inverted chairs, and challenges the “textbook” picture of conformationally locked carbohydrate rings. © 2015 Wiley Periodicals, Inc.

A revised version 56A6CARBO_R of the GROMOS 56A6CARBO force field for hexopyranose-based carbohydrates is proposed. The revision significantly improves the description of ring-conformational properties for residues within chains. In the revised force field, the 4C1 chair is the most stable conformation for all d-hexopyranose residues, irrespective of the linkage type and anomery, and of the position of the residue along the chain. Additionally, the influence of the functionalization type (alkylation or glycosylation) of the anomeric oxygen atom on the ring-inversion properties is discussed.

Combined quantum mechanical calculations and classical molecular dynamics simulations were conducted to investigate the hydration properties of carboxybetaine zwitterion brushes with varying separation distances between the quaternary ammonium cation and carboxylic anion. The brushes consist of zwitterion trimers and are investigated to mimic interacting zwitterion chains grafted on a substrate as well as polymers with interacting zwitterion side chains. Our results show that the values of both positive and negative charges, their separation distances as well as chain interactions appear to play a critical role in the hydration properties of the zwitterions. The overall hydration property of these zwitterions is dictated by the competition between the strong hydration of the charged groups and the dehydration of the hydrocarbon chains. The strongest hydration occurs when the CH2− unit in the hydrocarbon chain reaches 6–8 for these trimers. Further increase in the hydrocarbon chain length to 10–14 leads to significant and sudden dehydration of the trimers. The water structure and the water residence time surrounding the zwitterions also demonstrate substantial alteration at this length scale. This hydrophilic-to-hydrophobic transition is induced by the hydrophobic interactions of the trimer chains. Our hydration results could explain the observed trend of the superiority of the methylated carbohydrates and poly(ethylene glycol) as antifouling materials compared to corresponding hydroxyl-terminated compounds. © 2015 Wiley Periodicals, Inc.

The average pair distribution functions (g(r)) correlates the C atoms on the hydrocarbon chains between the cationic and anionic groups and the O atoms in water as a function of the carbon spacer chain length (CSL) on a carboxybetaine trimer. The sudden decrease in the g(r) function correlates to the hydrophilic-to-hydrophobic transition when the CSL increases from 12 to 14.

When an electric field is applied to an insulating membrane, movement of charged particles through a nanopore is induced. The measured ionic current reports on biomolecules passing through the nanopore. In this work, we explored the kinetics of DNA unzipping in a nanopore using our coarse-grained model (Stachiewicz and Molski, J. Comput. Chem. 2015, 36, 947). Coarse graining allowed a more detailed analysis for a wider range of parameters than all-atom simulations. Dependence of the translocation mode (unzipping or distortion) on the pore diameter was examined, and the threshold voltages were estimated. We determined the potential of mean force, position-dependent diffusion coefficient, and position-dependent effective charge for the DNA unzipping. The three molecular profiles were correlated with the ionic current and molecular events. On the unzipping/translocation force profile, two energy maxima were found, one of them corresponding to the unzipping, and the other to the translocation barriers. The unzipping kinetics were further explored using Brownian dynamics. © 2015 Wiley Periodicals, Inc.

When an electric field is applied to an insulating membrane, movement of charged particles through a nanopore is induced. In this work, the kinetics of DNA unzipping in a nanopore is explored using a coarse-grained model. Dependence of the translocation mode on the pore diameter is examined, and the threshold voltages are estimated. The potential of mean force, position-dependent diffusion coefficient and position-dependent effective charge are determined for the DNA unzipping.

Noncovalent functionalization of buckybowls sumanene (S), corannulene (R), and coronene (C) with greenhouse gases (GGs) such as CO2, CH4 (M), and C2H2 (A) has been studied using hybrid density functional theory. The propensity and preferences of these small molecules to interact with the concave and convex surfaces of the buckybowls has been quantitatively estimated. The results indicate that curvature plays a significant role in the adsorption of these small molecules on the π surface and it is observed that buckybowls have higher binding energies (BEs) compared with their planar counterpart coronene. The concave surface of the buckybowl is found to be more feasible for adsorption of small molecules. BEs of small molecules towards π systems is CO2 > A > M and the BEs of π systems toward small molecules is S > R > C. Obviously, the binding preference is dictated by the way in which various noncovalent interactions, such as π···π, lone pair···π, and CH···π manifest themselves on carbaneous surfaces. To delineate the intricate details of the interactions, we have employed Bader's quantum theory of atoms in molecule and localized molecular orbital energy decomposition analysis (LMO-EDA). LMO-EDA, which measures the contribution of various components and traces the physical origin of the interactions, indicates that the complexes are stabilized largely by dispersion interactions. © 2015 Wiley Periodicals, Inc.

Computational studies reveal that monolayer carbonaceous surfaces, such as sumanene, corannulene and coronene, are optimal for adsorption of small gas molecules such as CO2, CH4 and C2H2. The concave surface of the buckybowl is preferred over the convex alternative for the adsorption as well as the selective binding of small molecules. Dispersion is largely responsible for binding of small molecules, buckybowls and coronene.

A zone-folding (ZF) approach is applied for the estimation of the phonon contributions to thermodynamic properties of carbon-and ZrS2-based nanotubes (NTs) of hexagonal morphology with different chiralities. The results obtained are compared with those from the direct calculation of the thermodynamic properties of NTs using PBE0 hybrid exchange-correlation functional. The phonon contribution to the stability of NTs proved to be negligible for the internal energy and small for the Helmholtz free energy. It is found that the ZF approach allows us an accurate estimation of phonon contributions to internal energy, but slightly overestimates the phonon contributions to entropy. © 2015 Wiley Periodicals, Inc.

In the case of nanotubes rolled up from the layers of layered compounds (such as graphite, ZrS2, or V2O5) the phonon contributions to the heat capacity and internal energy calculated directly and estimated with using the appropriate 2D layer supercell (within the zone-folding approach) remain very close to each other for temperatures up to 600 K.

A set of 42 molecules with N-F, O-F, N-Cl, P-F, and As-F bonds has been investigated in the search for potential bond anomalies, which lead to reverse bond length–bond strength (BLBS) relationships. The intrinsic strength of each bond investigated has been determined by the local stretching force constant obtained at the CCSD(T)/aug-cc-pVTZ level of theory. N-F or O-F bond anomalies were found for fluoro amine radicals, fluoro amines, and fluoro oxides, respectively. A rationale for the deviation from the normal Badger-type inverse BLBS relation is given and it is shown that electron withdrawal accompanied by strong orbital contraction and bond shortening is one of the prerequisites for a bond anomaly. In the case of short electron-rich bonds such as N-F or O-F, anomeric delocalization of lone pair electrons in connection with lone pair repulsion are decisive whether a bond anomaly can be observed. This is quantitatively assessed with the help of the CCSD(T) local stretching force constants, CCSD(T) charge distributions, and G4 bond dissociation energies. Bond anomalies are not found for fluoro phosphines and fluoro arsines because the bond weakening effects are no longer decisive. © 2015 Wiley Periodicals, Inc.

Reverse bond length/bond strength relations have been reported in the literature, in particular for chemical bonds between electron rich atoms, e.g. NF or OF bonds. In this work, a comprehensive rational is derived covering all electronic and electrostatic factors that may lead to shorter but weaker bonds. A key feature is the use of a qualified bond strength measure based on vibrational spectroscopy.

We investigated the performance of the density functional theory (DFT) functionals PBE, PBE0, M06, and M06-L for describing the molecular and dissociative adsorption of O2 onto pure and doped Al(111) surfaces. Adsorption of O2 was studied at the perfect Al(111) surface and compared with the case where an additional Al atom was present as an adatom. Additionally, we studied how these functionals perform when different dopants are present at the Al(111) surface in two distinct geometries: as an adatom or as a substitutional atom replacing an Al atom. The performance of the different functionals is greatly affected by the surface geometry. The inclusion of Hartree-Fock exchange in the functional leads to slight differences in adsorption energies for molecular adsorption of O2. These differences become very pronounced for dissociative adsorption, with the hybrids PBE0 and M06 predicting more exergonic adsorption than PBE and M06-L. Furthermore, PBE0 and M06 predicted trends in adsorption energies for defective and perfect surfaces which are in line with the experimental knowledge of the effects of surface defects in adsorption energies. The predictions of the non-hybrids PBE and M06-L point in the opposite direction. The analysis of the contributions of the van der Waals (vdW) forces to the adsorption energies reveals that the PBE and PBE0 functionals have similar difficulties in describing vdW interactions for molecular adsorption of O2 while the M06 functional can give a description of these forces with an accuracy which is at least similar to that of the correction of the D3 type. © 2015 Wiley Periodicals, Inc.

Within density functional theory (DFT), different functionals produce very large differences in the description of the molecular and dissociative adsorption of O2 at pure and doped aluminum surfaces. These differences are here reported and discussed in view of the type of physical descriptors of the electron density incorporated in each functional.

A new type of reaction pathway which involves a nontotally symmetric trifurcation was found and investigated for a typical SN2-type reaction, NC- + CH3X NCCH3 + X- (X = F, Cl). A nontotally symmetric valley-ridge inflection (VRI) point was located along the C3v reaction path. For X = F, the minimum energy path (MEP) starting from the transition state (TS) leads to a second-order saddle point with C3v symmetry, which connects three product minima of Cs symmetry. For X = Cl, four product minima have been observed, of which three belong to Cs symmetry and one to C3v symmetry. The branching path from the VRI point to the lower symmetry minima was determined by a linear interpolation technique. The branching mechanism is discussed based on the reaction path curvature and net atomic charges, and the possibility of a nonotally symmetric n-furcation is discussed. © 2015 Wiley Periodicals, Inc.

A new type of reaction pathway which involves a nontotally symmetric trifurcation was found and investigated for the SN2-type reaction, NC- + CH3X NCCH3 + X- (X = F, Cl). A nontotally symmetric valley-ridge inflection (VRI) point was located along the C3v reaction path. The branching path from the VRI point to the lower symmetry minima was determined. The possibility of a nontotally symmetric n-furcation is discussed.

Molecular mechanics (MM4) studies have been carried out on the catenane (C13H26)2, specifically 13-13D2. The structure obtained is in general agreement with second-order perturbation theory. More importantly, the MM4 structure allows a breakdown of the energy of the molecule into its component classical parts. This allows an understanding of why the structure is so distorted, in terms of CC bonding and nonbonding interactions, van der Waals repulsion, CCC and CCH angle bending, torsional energies, stretch-bend, torsion-stretch, and bend–torsion–bend interactions. Clearly, the hole in 113-membered ring is too small for the other ring to fit through comfortably. There are too many atoms trying to fit into the limited space at the same time, leading to large van der Waals repulsions. The rings distort in such a way as to enlarge this available space, and lower the total energy of the molecule. While the distortions are spread around the rings, one of the nominally tetrahedral CCC bond angles in each ring is opened to 147.9° by MM4 (146.8° by MP2). The stability of the compound is discussed in terms of the strain energy. © 2015 Wiley Periodicals, Inc.

Catenanes have become important molecular units for the design of new materials for 21st century technology. The origins of strain in the simplest viable saturated hydrocarbon knots is investigated here. The combination of quantum chemistry with molecular mechanics provides many new insights.

DAMQT-2.1.0 is a new version of DAMQT package which includes topographical analysis of molecular electron density (MED) and molecular electrostatic potential (MESP), such as mapping of critical points (CPs), creating molecular graphs, and atomic basins. Mapping of CPs is assisted with algorithmic determination of Euler characteristic in order to provide a necessary condition for locating all possible CPs. Apart from the mapping of CPs and determination of molecular graphs, the construction of MESP-based atomic basin is a new and exclusive feature introduced in DAMQT-2.1.0. The GUI in DAMQT provides a user-friendly interface to run the code and visualize the final outputs. MPI libraries have been implemented for all the tasks to develop the parallel version of the software. Almost linear scaling of computational time is achieved with the increasing number of processors while performing various aspects of topography. A brief discussion of molecular graph and atomic basin is provided in the current article highlighting their chemical importance. Appropriate example sets have been presented for demonstrating the functions and efficiency of the code. © 2015 Wiley Periodicals, Inc.

DAMQT-2.1.0, a new version of DAMQT package, includes topographical analysis of molecular electron density and molecular electrostatic potential viz. mapping of critical points, creating molecular graphs and atomic basins. Determination of MESP-based topo-graph and MESP-based atomic basins are new and exclusive features introduced in DAMQT-2.1.0. The DAMQT package is provided with a user-friendly graphical user interface to run the code and visualize the final outputs.

A large number of fully halogenated benzene derivatives containing the fluorine, chlorine, bromine, and iodine atoms have been experimentally synthesized both as single- and co-crystals (e.g., Desiraju et al., Chem. Eur. J. 2006, 12, 2222), yet the natures of the halogen ··· halogen interactions between the vicinal halogens in these compounds within the intramolecular domain are undisclosed. Given a fundamental understanding of these interactions is incredibly important in many areas of chemical, biological, supramolecular, and material sciences, we present here our newly discovered theoretical results that delineate whilst the nature of an F···F interaction in a pair of two adjacent fluorine atoms in either of the hexafluorobenzene and 1,4-dibromotetrafluorobenzene compounds examined is almost unclear, each of the latter three hexahalogenated benzene derivatives (viz., C6Cl6, C6Br6, and C6I6), and each of the seven of their fully mixed hexahalogenated benzene analogues, are found to be stabilized by means of a number of halogen···halogen interactions, each a form of long-range attraction within the intramolecular domain. The Molecular Electrostatic Surface Potential model was found to be unsurprisingly unsuitable in unraveling any of the aforesaid attractions between the halogen atoms. However, such interactions successfully enunciated by a set of noncovalent interaction descriptors of geometrical, topological, and electrostatic origins. These latter properties were extracted combining the results of the Density Functional Theory electronic structure calculations with those revealed from Atoms in Molecules, and Reduced Density Gradient charge density-based topological calculations, and are expounded in detail to formalize the conclusions. © 2015 Wiley Periodicals, Inc.

Combined DFT, and QTAIM- and RDG-based topological charge density studies were performed to uncover the nature of the previously overlooked halogen–halogen intramolecular interactions between vicinal halogens in four fully halogenated and eight of their mixed benzene derivatives. The efficacy and usefulness of the latter two theoretical models, in contrast to the popular molecular electrostatic surface potential and isodesmic reaction models, to quantify intramolecular noncovalent interactions have demonstrated.

On page 2271 (DOI: 10.1002/jcc.24195), Ivan Carnimeo, Chiara Cappelli, and Vincenzo Barone present and validate a polarizable QM/MM scheme for the calculation of excited state properties of large molecular systems, providing both the theoretical framework and pilot test cases. Analytical post-SCF gradients for excited states at the TD-DFT level, as well as for ground state MP2 and B2PLYP Hamiltonians, have been implemented, thus allowing geometry optimization and the evaluation of vibronic spectra of complex molecular systems. In the image a graphical representation of a typical molecular system, which can be treated with this new approach, is shown.

The benzene-benzene interaction (Bz-Bz), tested on ab initio calculations of the intermolecular potential energy and given in a suitable analytical form, is used to formulate the multidimensional potential energy surface of Bzn clusters of growing size. On page 2291 (DOI: 10.1002/jcc.24201), Massimiliano Bartolomei, Fernando Pirani, and Jorge M.C. Marques adopt an evolutionary algorithm to identify the energy and structure of low lying isomers of the clusters, including the main features of the most stable configurations.

Using current state-of-the-art computational approaches (DFT, QTAIM, RDG-NCI, DI, etc.), Pradeep R. Varadwaj, Arpita Varadwaj, and Bih-Yaw Jin demonstrate on page 2328 (DOI: 10.1002/jcc.24211) that vicinal halogen atoms in each of the three fully halogenated benzene molecules involving the latter three halogens X (X = Cl, Br, I), as well as eight of their fully mixed halogenated benzene derivatives including the fluorine atom, are linked with each other via the X···X intramolecular interaction. This leads to an interpretation that the presence of such prototypical interactions are essentially important as they might contribute to the overall stabilities of each of such compounds.

Intermolecular ternary complexes composed of: (1) the centrally placed trifluoroacetonitrile or its higher analogs with central carbon exchanged by silicon or germanium (M = C, Si, Ge), (2) the benzonitrile molecule or its para derivatives on one side, and (3) the boron trifluoride of trichloride molecule (X = F, Cl) on the opposite side as well as the corresponding intermolecular tetrel- and triel-bonded binary complexes, were investigated by symmetry-adapted perturbation theory (SAPT) and the supermolecular Møller-Plesset method (MP2) at the complete basis set limit for optimized geometries. A character of interactions was studied by quantum theory of atoms-in-molecules (QTAIM). A comparison of interaction energies and QTAIM bond descriptors for dimers and trimers reveals that tetrel and triel bonds increase in their strength if present together in the trimer. For the triel-bonded complex, this growth leads to a change of the bond character from closed-shell to partly covalent for Si or Ge tetrel atoms, so the resulting bonding scheme corresponds to a preliminary stage of the 2 reaction. Limitations of the Lewis theory of acids and bases were shown by its failure in predicting the stability order of the triel complexes. The necessity of including interaction energy terms beyond the electrostatic component for an elucidation of the nature of σ- and π-holes was presented by a SAPT energy decomposition and by a study of differences in monomer electrostatic potentials obtained either from isolated monomer densities, or from densities resulting from a perturbation with the effective field of another monomer. © 2015 Wiley Periodicals, Inc.

Interplay of noncovalent interactions for ternary complexes is important from both practical (designing new materials) and theoretical points of view. In this contribution we employ high-quality computational methods to study cooperativity of two bonding types: tetrel (being a preliminary step in the SN2 reaction) and triel (interesting due to the hypovalency of boron). We found that triel bondings change their character from noncovalent to partly covalent under influence of the third molecule in the complex.

Hartree–Fock and density functional theory with the hybrid B3LYP and general gradient KT2 exchange-correlation functionals were used for nonrelativistic and relativistic nuclear magnetic shielding calculations of helium, neon, argon, krypton, and xenon dimers and free atoms. Relativistic corrections were calculated with the scalar and spin-orbit zeroth-order regular approximation Hamiltonian in combination with the large Slater-type basis set QZ4P as well as with the four-component Dirac–Coulomb Hamiltonian using Dyall's acv4z basis sets. The relativistic corrections to the nuclear magnetic shieldings and chemical shifts are combined with nonrelativistic coupled cluster singles and doubles with noniterative triple excitations [CCSD(T)] calculations using the very large polarization-consistent basis sets aug-pcSseg-4 for He, Ne and Ar, aug-pcSseg-3 for Kr, and the AQZP basis set for Xe. For the dimers also, zero-point vibrational (ZPV) corrections are obtained at the CCSD(T) level with the same basis sets were added. Best estimates of the dimer chemical shifts are generated from these nuclear magnetic shieldings and the relative importance of electron correlation, ZPV, and relativistic corrections for the shieldings and chemical shifts is analyzed. © 2015 Wiley Periodicals, Inc.

The amount of relativistic contribution to nuclear shielding is the same for noble gas dimers and free atoms. For Ne to Xe, the relativistic contribution is more important than both electron correlation and vibrational corrections. However, for the chemical shifts the relativistic corrections almost cancel, so that ZPVC become more important than relativistic corrections for Ne2, Kr2, and Ar2.

A computational protein design method is extended to allow Monte Carlo simulations where two ligands are titrated into a protein binding pocket, yielding binding free energy differences. These provide a stringent test of the physical model, including the energy surface and sidechain rotamer definition. As a test, we consider tyrosyl-tRNA synthetase (TyrRS), which has been extensively redesigned experimentally. We consider its specificity for its substrate l-tyrosine (l-Tyr), compared to the analogs d-Tyr, p-acetyl-, and p-azido-phenylalanine (ac-Phe, az-Phe). We simulate l- and d-Tyr binding to TyrRS and six mutants, and compare the structures and binding free energies to a more rigorous “MD/GBSA” procedure: molecular dynamics with explicit solvent for structures and a Generalized Born + Surface Area model for binding free energies. Next, we consider l-Tyr, ac- and az-Phe binding to six other TyrRS variants. The titration results are sensitive to the precise rotamer definition, which involves a short energy minimization for each sidechain pair to help relax bad contacts induced by the discrete rotamer set. However, when designed mutant structures are rescored with a standard GBSA energy model, results agree well with the more rigorous MD/GBSA. As a third test, we redesign three amino acid positions in the substrate coordination sphere, with either l-Tyr or d-Tyr as the ligand. For two, we obtain good agreement with experiment, recovering the wildtype residue when l-Tyr is the ligand and a d-Tyr specific mutant when d-Tyr is the ligand. For the third, we recover His with either ligand, instead of wildtype Gln. © 2015 Wiley Periodicals, Inc.

A computational protein design method is extended to allow Monte Carlo simulations where two ligands are titrated into a protein binding pocket, yielding binding free energy differences. These provide a stringent test of the physical model, including the energy surface and sidechain rotamer definition. As a test, tyrosyl-tRNA synthetase and its specificity for its substrate l-Tyr, compared to d-Tyr, p-acetyl-, and p-azido-phenylalanine is considered.