Towards a force field based on density fitting

Total intermolecular interaction energies are determined with a first version of the Gaussian electrostatic model (GEM-0), a force field based on a density fitting approach using s-type Gaussian functions. The total interaction energy is computed in the spirit of the sum of interacting fragment ab initio (SIBFA) force field by separately evaluating each one of its components: electrostatic (Coulomb), exchange repulsion, polarization, and charge transfer intermolecular interaction energies, in order to reproduce reference constrained space orbital variation (CSOV) energy decomposition calculations at the B3LYP/aug-cc-pVTZ level. The use of an auxiliary basis set restricted to spherical Gaussian functions facilitates the rotation of the fitted densities of rigid fragments and enables a fast and accurate density fitting evaluation of Coulomb and exchange-repulsion energy, the latter using the overlap model introduced by Wheatley and Price [Mol. Phys. 69, 50718 (1990)]. The SIBFA energy scheme for polarization and charge transfer has been implemented using the electric fields and electrostatic potentials generated by the fitted densities. GEM-0 has been tested on ten stationary points of the water dimer potential energy surface and on three water clusters (n = 16,20,64). The results show very good agreement with density functional theory calculations, reproducing the individual CSOV energy contributions for a given interaction as well as the B3LYP total interaction energies with errors below kBT at room temperature. Preliminary results for Coulomb and exchange-repulsion energies of metal cation complexes and coupled cluster singles doubles electron densities are discussed.     

Numerical fitting of molecular properties to Hermite Gaussians

A procedure is presented to fit gridded molecular properties to auxiliary basis sets (ABSs) of Hermite Gaussians, analogous to the density fitting (DF) method (Dunlap; et al. J. Chem. Phys. 1979, 71, 4993). In this procedure, the ab initio calculated properties (density, electrostatic potential, and/or electric field) are fitted via a linear- or nonlinear-least-squares procedure to auxiliary basis sets (ABS). The calculated fitting coefficients from the numerical grids are shown to be more robust than analytic density fitting due to the neglect of the core contributions. The fitting coefficients are tested by calculating intermolecular Coulomb and exchange interactions for a set of dimers. It is shown that the numerical instabilities observed in DF are caused by the attempt of the ABS to fit the core contributions. In addition, this new approach allows us to reduce the number of functions required to obtain an accurate fit. This results in decreased computational cost, which is shown by calculating the Coulomb energy of a 4096 water box in periodic boundary conditions. Using atom centered Hermite Gaussians, this calculation is only 1 order of magnitude slower than conventional atom-centered point charges.     

Poisson-transformed density fitting in relativistic four-component Dirac-Kohn-Sham theory

We present recent developments in the implementation of the density fitting approach for the Coulomb interaction within the four-component formulation of relativistic density functional theory [Belpassi et al., J. Chem. Phys. 124, 124104 (2006)]. In particular, we make use of the Poisson equation to generate suitable auxiliary basis sets and simplify the electron repulsion integrals [Manby and Knowles, Phys. Rev. Lett. 87, 163001 (2001)]. We propose a particularly simple and efficient method for the generation of accurate Poisson auxiliary basis sets, based on already available standard Coulomb fitting sets. Just as is found in the nonrelativistic case, we show that the number of standard auxiliary fitting functions that need to be added to the Poisson-generated functions in order to achieve a fitting accuracy equal or, in some cases, better than that of the standard procedure is remarkably small. The efficiency of the present implementation is demonstrated in a detailed study of the spectroscopic properties and energetics of several gold containing systems, including the Au dimer and the CsAu molecule. The extraction reaction of a H(2)O molecule from a Au(H(2)O)(9) (+) cluster is also calculated as an example of mixed heavy-light-atom molecular systems. The scaling behavior of the algorithm implemented is illustrated for some closed shell gold clusters up to Au(5) (+). The increased sparsity of the Coulomb matrices involved in the Poisson fitting is identified, as are potential computational applications and the use of the Poisson fitting for the relativistic exchange-correlation problem.     

Quantum-chemical study of CHCl3-SO2 association

CHCl(3)-SO(2) association is studied by high-level quantum-chemical calculations of stationary points of the dimer electronic potential-energy hypersurface, including correlated second-order Moller-Plesset and CCSD(T) calculations with basis sets up to 6-311++G(d,p). During geometry optimization, frequency, and energy calculations, a self-written computer code embedding the GAMESS ab initio program suite applies counterpoise correction of the basis set superposition error. A CH...O hydrogen-bonded complex (DeltaE(0)=-8.73 kJmol) with a 2.4 A intermolecular H...O distance and two very weak van der Waals complexes (DeltaE(0)=-3.78 and -2.94 kJmol) are located on the counterpoise-corrected potential-energy surface. The intermolecular interactions are characterized by Kitaura-Morokuma interaction energy decompositions and Mulliken electron population analyses. The unusual hydrogen bond is distinguished by a CH-bond contraction, a pronounced enhancement of the IR intensity and a shift to higher frequency ("blueshift") of the CH-stretching vibration compared to the CHCl(3) monomer. Spectroscopy and association in liquid solution is also discussed; our results provide an alternative explanation for features in the CH-stretching vibration spectrum of chloroform dissolved in liquid sulfur dioxide which have been attributed previously to an intermolecular Fermi resonance.     

Validation of intermolecular transfer integral and bandwidth calculations for organic molecular materials

We present an interpretation of the intermolecular transfer integral that is independent from the origin of the energy scale allowing convergence studies of this important parameter of organic molecular materials. We present extensive numerical studies by using an ethylene pi dimer to investigate the dependence of transfer integrals on the level of theory and intermolecular packing. Transfer integrals obtained from semiempirical calculations differ substantially from one another and from ab initio results. The ab initio results are consistent across all the levels used including Hartree-Fock, outer valence Green's function, and various forms of density functional theory (DFT). Validation of transfer integrals and bandwidths is performed by comparing the calculated values with the experimental values of tetrathiafulvalene-tetracyanoquinodimethane (TTF-TCNQ), bis[1,2,5]thiadiazolo-p-quinobis(1,3-dithiole), (BTQBT) K-TCNQ, and hexagonal graphite. DFT in one of its presently popular forms, such as Perdew-Wang functionals (PW91), in combination with sufficient basis sets provides reliable transfer integrals, and therefore can serve as a basis for energy band calculations for soft organic materials with van der Waals gaps.     

A self-consistent coulomb bath model using density fitting

A self-consistent Coulomb bath model is presented to provide an accurate and efficient way of performing calculations for interfragment electrostatic and polarization interactions. In this method, a condensed-phase system is partitioned into molecular fragment blocks. Each fragment is embedded in the Coulomb bath due to other fragments. Importantly, the present Coulomb bath is represented using a density fitting method in which the electron densities of molecular fragments are fitted using an atom-centered auxiliary basis set of Gaussian type. The Coulomb bath is incorporated into an effective Hamiltonian for each fragment, with which the electron density is optimized through an iterative double self-consistent field (DSCF) procedure to realize the mutual many-body polarization effects. In this work, the accuracy of interfragment interaction energies enumerated using the Coulomb bath is tested, showing a good agreement with the exact results from an energy decomposition analysis. The qualitative features of many-body polarization effects are visualized by electron density difference plots. It is also shown that the present DSCF method can yield fast and robust convergence with near-linear scaling in performance with increase in system size.     

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.     

Critical examination of the supermolecule density functional theory calculations of intermolecular interactions

The results of calculations employing twelve different combinations of exchange and correlation functionals are compared with results of ab initio calculations for two different configurations of the water dimer and three different configurations of the thymine-adenine complex. None of the density functional theory (DFT) treatments could properly reproduce the results of coupled-cluster calculations for all configurations examined. The DFT approaches perform well when the interaction energy is dominated by the electrostatic component and the dispersion energy is less important. Two mechanisms that compensate for the missing dispersion component were identified. The first one is the decrease of the magnitude of the intermolecular exchange-repulsion and the second one is the increase of the magnitude of the attractive deformation energy. For some functionals both effects are observed together, but for some other ones only the second effect occurs. The three correlation functionals that were examined were found to make only very small contributions to the deformation energy. The examination of angular and distance dependence of the interactions shows that the currently available DFT approaches are not suitable for developing intermolecular potential energy surfaces. They could however be used to find global minima on potential energy surfaces governed by intermolecular electrostatic interactions. Additional single point ab initio calculations are recommended as the means of validating optimized structures.     

Evaluation of the charge penetration energy between non-orthogonal molecular orbitals using the Spherical Gaussian Overlap approximation

An overlap dependent formula for evaluating the charge penetration energy between non-orthogonal molecular orbitals is derived using the Spherical Gaussian Overlap approximation. When combined with an accurate multipole representation of the electrostatic energy, such as in the effective fragment potential method, ab initio electrostatic energies are generally reproduced to within 0.2 kcal/mol for a variety of molecular dimers and basis sets. The only larger error is for the DMSO dimer, where the electrostatic energy is overestimated by 0.7 kcal/mol.     

Calculation of the O-H stretching vibrational overtone spectrum of the water dimer

The O-H stretching vibrational overtone spectrum of the water dimer has been calculated with the dimer modeled as two individually vibrating monomer units. Vibrational term values and absorption intensities have been obtained variationally with a computed dipole moment surface and an internal coordinate Hamiltonian, which consists of exact kinetic energy operators within the Born-Oppenheimer approximation of the monomer units. Three-dimensional ab initio potential energy and dipole moment surfaces have been calculated using the internal coordinates of the monomer units using the coupled cluster method including single, double, and perturbative triple excitations [CCSD(T)] with the augmented correlation consistent valence triple zeta basis set (aug-cc-pVTZ). The augmented correlation consistent valence quadruple zeta basis set (aug-cc-pVQZ), counterpoise correction, basis set extrapolation to the complete basis set limit, relativistic corrections, and core and valence electron correlations effects have been included in one-dimensional potential energy surface cuts. The aim is both to investigate the level of ab initio and vibrational calculations necessary to produce accurate results when compared with experiment and to aid the detection of the water dimer under atmospheric conditions.     

Interaction energies in hydrogen-bonded systems: a testing ground for subsystem formulation of density-functional theory

The formalism based on the total energy bifunctional (E[rhoI,rhoII]) is used to derive interaction energies for several hydrogen-bonded complexes (water dimer, HCN-HF, H2CO-H2O, and MeOH-H2O). Benchmark ab initio data taken from the literature were used as a reference in the assessment of the performance of gradient-free [local density approximation (LDA)] and gradient-dependent [generalized gradient approximation (GGA)] approximations to the exchange-correlation and nonadditive kinetic-energy components of E[rhoI,rhoII]. On average, LDA performs better than GGA. The average absolute error of calculated LDA interaction energies amounts to 1.0 kJ/mol. For H2CO-H2O and H2O-H2O complexes, the potential-energy curves corresponding to the stretching of the intermolecular distance are also calculated. The positions of the minima are in a good agreement (less than 0.2 A) with the reference ab initio data. Both variational and nonvariational calculations are performed to assess the energetic effects associated with complexation-induced deformations of molecular electron densities.     

Electron density fitting for the Coulomb problem in relativistic density-functional theory

A density fitting approach for the Coulomb matrix representation within the four-component formulation of relativistic density-functional theory is presented. Our implementation, which uses G-spinor basis sets, shares all the advantages of those found in nonrelativistic quantum chemistry. We show that very accurate Coulomb energies may be obtained using a modest number of scalar auxiliary basis functions for molecules containing heavy atoms. The efficiency of this new implementation is demonstrated in a detailed study of the spectroscopic properties of the gold dimer, and its scaling behavior has been tested by calculations of some closed-shell gold clusters (Au2, Au3+, Au4, Au5+). The algorithm is found to scale as O(N3), just as it does in the nonrelativistic case, and represents a dramatic improvement in efficiency over the conventional approach in the calculation of the Coulomb matrix, with computation times that are reduced to less than 3% for Au2 and up to 1% in the case of Au5+.     

Role of electrostatic interactions for the domain shapes of Langmuir monolayers of monoglycerol amphiphiles

The role of electrostatic interaction in the domain morphology of amide, ether, ester, and amine monoglycerol monolayers (abbreviated as ADD, ETD, ESD, and AMD, respectively) with systematic variation in the molecular structure of the headgroup region is investigated. Experimental studies using Brewster angle microscopy (BAM) and grazing incidence X-ray diffraction (GIXD) show that the characteristic features of the condensed monolayer phase, such as domain morphology, crystallinity, and lattice parameters, are very different for these monoglycerols. Therefore, the intermolecular interactions of the four amphiphilic monoglycerols are investigated in detail. First, the dipole moments of four monoglycerols of similar structure but with different functional groups are calculated by a semiempirical quantum mechanical technique. The dipole moments for monoglycerols follow the sequence AMD < ETD < ESD < ADD for the population of conformers of compounds investigated. The dipolar repulsion energies for the amphiphilic monoglycerols are also calculated for different possible mutual orientations between the dipoles. The calculated dipolar energies also follow the same trend for different possible headgroup orientations. These results can explain the domain shape of the monoglycerols observed experimentally. Second, ab initio calculations on the basis of the HF/6-31G** method are performed for representative monoglycerol headgroup segments. The results show that the intermolecular interaction energy related to dimer formation follows the order ETD < ESD < AMD < ADD segments, similar to that observed in experiment except in the case of the AMD segment. The relative importance of intra- and intermolecular hydrogen bonding in dimers is analyzed. The enhanced role of the intermolecular interaction relative to intramolecular interaction in the case of AMD contributes to the relatively high intermolecular interaction energy for the particular conformation of the dimer of AMD segment as observed from ab initio calculation. The present work shows that the variations in headgroup molecular structure alter drastically the domain shape, and the theoretical calculations conclusively reveal the important role of the electrostatic interactions for the mesoscopic domain architecture.     

Experimental/theoretical electrostatic properties of a styrylquinoline-type HIV-1 integrase inhibitor and its progenitors

We have established that polyhydroxylated styrylquinolines are potent inhibitors of HIV-1 integrase (IN). Among them, we have identified (E)-8-hydroxy-2-[2-(4,5-dihydroxy-3-methoxyphenyl)-ethenyl]-7-quinolinecarboxylic acid (1) as a promising lead. Previous molecular dynamics simulations and docking procedures have shown that the inhibitory activity involves one or two metal cations (Mg2+), which are present in the vicinity of the active center of the enzyme. However, such methods are generally based on a force-field approach and still remain not as reliable as ab initio calculations with extended basis sets on the whole system. To go further in this area, the aim of the present study was to evaluate the predictive ability of the electron density and electrostatic properties in the structure-activity relationships of this class of HIV-1 antiviral drugs. The electron properties of the two chemical progenitors of 1 were derived from both high-resolution X-ray diffraction experiments and ab initio calculations. The twinning phenomenon and solvent disorder were observed during the crystal structure determination of 1. Molecule 1 exhibits a planar s-trans conformation, and a zwitterionic form in the crystalline state is obtained. This geometry was used for ab initio calculations, which were performed to characterize the electronic properties of 1. The electron densities, electrostatic potentials, and atomic charges of 1 and its progenitors are here compared and analyzed. The experimental and theoretical deformation density bond peaks are very comparable for the two progenitors. However, the experimental electrostatic potential is strongly affected by the crystal field and cannot straightforwardly be used as a predictive index. The weak difference in the theoretical electron densities between 1 and its progenitors reveals that each component of 1 conserves its intrinsic properties, an assumption reinforced by a 13C NMR study. This is also shown through an excellent correlation of the atomic charges for the common fragments. The electrostatic potential minima in zwitterionic and nonzwitterionic forms of 1 are discussed in relation with the localization of possible metal chelation sites.     

Progress in the understanding of drug-receptor interactions, part 2: experimental and theoretical electrostatic moments and interaction energies of an angiotensin II receptor antagonist (C30H30N6(O)3S)

A combined experimental and theoretical charge density study of an angiotensin II receptor antagonist (1) is presented focusing on electrostatic properties such as atomic charges, molecular electric moments up to the fourth rank and energies of the intermolecular interactions, to gain an insight into the physical nature of the drug-receptor interaction. Electrostatic properties were derived from both the experimental electron density (multipole refinement of X-ray data collected at T=17 K) and the ab initio wavefunction (single molecule and fully periodic calculations at the DFT level). The relevance of SO and SN intramolecular interactions on the activity of 1 is highlighted by using both the crystal and gas-phase geometries and their electrostatic nature is documented by means of QTAIM atomic charges. The derived electrostatic properties are consistent with a nearly spherical electron density distribution, characterised by an intermingling of electropositive and -negative zones rather than by a unique electrophilic region opposed to a nucleophilic area. This makes the first molecular moment scarcely significant and ill-determined, whereas the second moment is large, significant and highly reliable. A comparison between experimental and theoretical components of the third electric moment shows a few discrepancies, whereas the agreement for the fourth electric moment is excellent. The most favourable intermolecular bond is show to be an NHN hydrogen bond with an energy of about 50 kJ mol(-1). Key pharmacophoric features responsible for attractive electrostatic interactions include CHX hydrogen bonds. It is shown that methyl and methylene groups, known to be essential for the biological activity of the drug, provide a significant energetic contribution to the total binding energy. Dispersive interactions are important at the thiophene and at both the phenyl fragments. The experimental estimates of the electrostatic contribution to the intermolecular interaction energies of six molecular pairs, obtained by a new model proposed by Spackman, predict the correct relative electrostatic energies with no exceptions.     

Intermolecular coulomb couplings from ab initio electrostatic potentials: application to optical transitions of strongly coupled pigments in photosynthetic antennae and reaction centers

An accurate and numerically efficient method for the calculation of intermolecular Coulomb couplings between charge densities of electronic states and between transition densities of electronic excitations is presented. The coupling of transition densities yields the F?rster type excitation energy transfer coupling, and from the charge density coupling, a shift in molecular excitation energies results. Starting from an ab initio calculation of the charge and transition densities, atomic partial charges are determined such as to fit the resulting electrostatic potentials of the different states and the transition. The different intermolecular couplings are then obtained from the Coulomb couplings between the respective atomic partial charges. The excitation energy transfer couplings obtained in the present TrEsp (transition charge from electrostatic potential) method are compared with couplings obtained from the simple point-dipole and extended dipole approximations and with those from the ab initio transition density cube method of Krüger, Scholes, and Fleming. The present method is of the same accuracy as the latter but computationally more efficient. The method is applied to study strongly coupled pigments in the light-harvesting complexes of green sulfur bacteria (FMO), purple bacteria (LH2), and higher plants (LHC-II) and the "special pairs" of bacterial reaction centers and reaction centers of photosystems I and II. For the pigment dimers in the antennae, it is found that the mutual orientation of the pigments is optimized for maximum excitonic coupling. A driving force for this orientation is the Coulomb coupling between ground-state charge densities. In the case of excitonic couplings in the "special pairs", a breakdown of the point-dipole approximation is found for all three reaction centers, but the extended dipole approximation works surprisingly well, if the extent of the transition dipole is chosen larger than assumed previously. For the "special pairs", a large shift in local transition energies is found due to charge density coupling.     

On the hydrogen bonding between water and ammonia molecules

The hydrogen bond N·HO between the water and ammonia molecules has been investigated ab initio using the SCF LCAO MO method. The minimal and extended basis sets of Slater type orbitals were used. It was found that the energy of the hydrogen bond is equal to 6.44 kcal/mole and the equilibrium separation of the oxygen and nitrogen atoms in the dimer is 5.72 au. At this intermolecular distance there is only one minimum in the potential energy curve for the motion of proton.     

The electronic structure of 2-chloro-2′-deoxyadenosine (drug) studied by Cl-NQR spectroscopy and ab initio calculations

Electronic structure and conformation of 2-chloro-2′-deoxyadenosine (2CldAdo) and 2-chloroadenosine (2ClAde) have been studied by ³⁵Cl-NQR spectroscopy and ab initio calculations by the B3LYP/6-31G^* method.

The usefulness of NQR spectroscopy and ab initio calculations for the investigation of redistribution of electron density and structural changes in 2CldAdo and 2ClAde was assessed and a comparison of the electron distribution in the unsubstituted 2ClAde (monomer, dimer) and the substituted 2ClAdo was made.     

An empirical potential for the O-H·O Hydrogen bond: Part I. Comparison with ab initio molecular orbital results

An empirical potential energy function has been devised for the O-H·O hydrogen bond, for use with the MMI force field. The energy of the hydrogen bond is described as the sum of van der Waals, electrostatic and Morse components. The function has been used to calculate the potential energy hypersurface of the water dimer, and the results are compared with published ab initio molecular orbital studies. Satisfactory agreement is obtained except for orientations involving very short H·H contacts. The geometry and hydrogen bond energy of the equilibrium linear form of (H₂O)₂ are calculated to be r(O·O) = 2.84 Å, θ = 36°, ΔE = ?5.35 kcal mol^?1, which are close to the values obtained by experiment, and from molecular orbital calculations. The relative importance of the electrostatic component of the empirical hydrogen bond energy is consistent with molecular orbital energy decomposition studies. The empirical function has also been used to calculate the energy of the water trimer in orientations which serve as models for the crystallographic bifurcated hydrogen bond. The results indicate that, in these orientations, the trimer is typically 0–3 kcal mol^?1 more stable than the dimer, a result which is consistent with ab initio calculations.