Ab Initio Studies of π···π Interactions: The Effects of Quadruple Excitations†

Accurate and efficient open-source implementation of domain-based local pair natural orbital (DLPNO) coupled-cluster theory using a t1-transformed Hamiltonian

We present an efficient, open-source formulation for coupled-cluster theory through perturbative triples with domain-based local pair natural orbitals [DLPNO-CCSD(T)]. Similar to the implementation of the DLPNO-CCSD(T) method found in the ORCA package, the most expensive integral generation and contraction steps associated with the CCSD(T) method are linear-scaling. In this work, we show that the t1-transformed Hamiltonian allows for a less complex algorithm when evaluating the local CCSD(T) energy without compromising efficiency or accuracy. Our algorithm yields sub-kJ mol-1 deviations for relative energies when compared with canonical CCSD(T), with typical errors being on the order of 0.1 kcal mol-1, using our TightPNO parameters. We extensively tested and optimized our algorithm and parameters for non-covalent interactions, which have been the most difficult interaction to model for orbital (PNO)-based methods historically. To highlight the capabilities of our code, we tested it on large water clusters, as well as insulin (787 atoms).

Tensor Hypercontraction Form of the Perturbative Triples Energy in Coupled-Cluster Theory

We present the working equations for a reduced-scaling method of evaluating the perturbative triples (T) energy in coupled-cluster theory, through the tensor hypercontraction (THC) of the triples amplitudes (t_ijk^abc). Through our method, we can reduce the scaling of the (T) energy from the traditional O(N7) to a more modest O(N5). We also discuss implementation details to aid future research, development, and software realization of this method. Additionally, we show that this method yields submillihartree (mEh) differences from CCSD(T) when evaluating absolute energies and sub-0.1 kcal/mol energy differences when evaluating relative energies. Finally, we demonstrate that this method converges to the true CCSD(T) energy through the systematic increasing of the rank or eigenvalue tolerance of the orthogonal projector, as well as exhibiting sublinear to linear error growth with respect to system size.

Machine Learning of Coupled Cluster (T)-Energy Corrections via Delta (Δ)-Learning

Accurate thermochemistry is essential in many chemical disciplines, such as astro-, atmospheric, or combustion chemistry. These areas often involve fleetingly existent intermediates whose thermochemistry is difficult to assess. Whenever direct calorimetric experiments are infeasible, accurate computational estimates of relative molecular energies are required. However, high-level computations, often using coupled cluster theory, are generally resource-intensive. To expedite the process using machine learning techniques, we generated a database of energies for small organic molecules at the CCSD(T)/cc-pVDZ, CCSD(T)/aug-cc-pVDZ, and CCSD(T)/cc-pVTZ levels of theory. Leveraging the power of deep learning by employing graph neural networks, we are able to predict the effect of perturbatively included triples (T), that is, the difference between CCSD and CCSD(T) energies, with a mean absolute error of 0.25, 0.25, and 0.28 kcal mol^-1 (R² of 0.998, 0.997, and 0.998) with the cc-pVDZ, aug-cc-pVDZ, and cc-pVTZ basis sets, respectively. Our models were further validated by application to three validation sets taken from the S22 Database as well as to a selection of known theoretically challenging cases.

Structures, binding energies and non-covalent interactions of furan clusters

Understanding of the furan solvent is subjected to the knowledge of the structures of the furan clusters and interactions taking place therein. Although, furan clusters can be very important to determine the dynamics and the properties of the furan solvent, there has been only a few investigations reported on furan dimer. In this work, we have explored the potential energy surfaces (PESs) of the furan clusters using two incremental levels of theory. Structures have been initially generated using classical molecular dynamics followed by full optimization at the MP2/aug-cc-pVDZ level of theory. The results show that the most stable structure of the furan dimer has a stacking configuration while that of the trimer has a cyclic configuration. We have noted that the structures of the furan tetramer have no definite configurations. In addition, we have performed a quantum theory of atoms in molecule (QTAIM) analysis to identify all possible non-covalent interactions of the furan clusters. The results show that six different types of non-covalent interactions can be identified in furan clusters. We have noted that the CH⋯C and CH⋯O hydrogen bondings are the strongest non-covalent interactions while the H⋯H bonding interaction is found to be the weakest. Furthermore, we have assessed the performance of ten DFT functionals in calculating the binding energies of the furan clusters. The ten DFT functionals (M05, M05-2X, M06, M06-2X, M08HX, PBE0, ωB97XD, PW6B95D3, APFD and MN15) have been benchmarked to DLPNO-CCSD(T)/CBS. The functionals M05-2X and M06 are recommended for further affordable investigations of the furan clusters.

Prototypical π-π dimers re-examined by means of high-level CCSDT(Q) composite ab initio methods

The benzene-ethene and parallel-displaced (PD) benzene-benzene dimers are the most fundamental systems involving π-π stacking interactions. Several high-level ab initio investigations calculated the binding energies of these dimers using the coupled-cluster with singles, doubles, and quasi-perturbative triple excitations [CCSD(T)] method at the complete basis set [CBS] limit using various approaches such as reduced virtual orbital spaces and/or MP2-based basis set corrections. Here, we obtain CCSDT(Q) binding energies using a Weizmann-3-type approach. In particular, we extrapolate the self-consistent field (SCF), CCSD, and (T) components using large heavy-atom augmented Gaussian basis sets [namely, SCF/jul-cc-pV{5,6}Z, CCSD/jul-cc-pV{Q,5}Z, and (T)/jul-cc-pV{T,Q}Z]. We consider post-CCSD(T) contributions up to CCSDT(Q), inner-shell, scalar-relativistic, and Born-Oppenheimer corrections. Overall, our best relativistic, all-electron CCSDT(Q) binding energies are ∆E_e,all,rel = 1.234 (benzene-ethene) and 2.550 (benzene-benzene PD), ∆H₀ = 0.949 (benzene-ethene) and 2.310 (benzene-benzene PD), and ∆H₂₉₈ = 0.130 (benzene-ethene) and 1.461 (benzene-benzene PD) kcal mol^-1. Important conclusions are reached regarding the basis set convergence of the SCF, CCSD, (T), and post-CCSD(T) components. Explicitly correlated calculations are used as a sanity check on the conventional binding energies. Overall, post-CCSD(T) contributions are destabilizing by 0.028 (benzene-ethene) and 0.058 (benzene-benzene) kcal mol^-1, and thus, they cannot be neglected if sub-chemical accuracy is sought (i.e., errors below 0.1 kcal mol^-1). CCSD(T)/aug-cc-pwCVTZ core-valence corrections increase the binding energies by 0.018 (benzene-ethene) and 0.027 (benzene-benzene PD) kcal mol^-1. Scalar-relativistic and diagonal Born-Oppenheimer corrections are negligibly small. We use our best CCSDT(Q) binding energies to evaluate the performance of MP2-based, CCSD-based, and lower-cost composite ab initio procedures for obtaining these challenging π-π stacking binding energies.

Compressed intramolecular dispersion interactions

The feasibility of the compression of localized virtual orbitals is explored in the context of intramolecular long-range dispersion interactions. Singular value decomposition (SVD) of coupled cluster doubles amplitudes associated with the dispersion interactions is analyzed for a number of long-chain systems, including saturated and unsaturated hydrocarbons and a silane chain. Further decomposition of the most important amplitudes obtained from these SVDs allows for the analysis of the dispersion-specific virtual orbitals that are naturally localized. Consistent with previous work on intermolecular dispersion interactions in dimers, it is found that three important geminals arise and account for the majority of dispersion interactions at the long range, even in the many body intramolecular case. Furthermore, it is shown that as few as three localized virtual orbitals per occupied orbital can be enough to capture all pairwise long-range dispersion interactions within a molecule.

A Linear-Scaling Method for Noncovalent Interactions: An Efficient Combination of Absolutely Localized Molecular Orbitals and a Local Random Phase Approximation Approach

A novel method for the accurate and efficient calculation of interaction energies in weakly bound complexes composed of a large number of molecules is presented. The new ALMO+RPAd method circumvents the prohibitive scaling of coupled cluster singles and doubles while still providing similar accuracy across a diverse range of weakly bound chemical systems. Linear-scaling procedures for the Fock build are given utilizing absolutely localized molecular orbitals (ALMOs), resulting in the a priori exclusion of basis set superposition errors. A bespoke data structure and algorithm using density fitting are described, leading to linear scaling for the storage and computation of the two-electron integrals. Electron correlation is included through a new, linear-scaling pairwise local random phase approximation approach, including exchange interactions, and decomposed into purely dispersive excitations (RPAxd). Collectively, these allow meaningful decomposition of the interaction energy into physically distinct contributions: electrostatic, polarization, charge transfer, and dispersion. Comparison with symmetry-adapted perturbation theory shows good qualitative agreement. Tests on various dimers and the S66 benchmark set demonstrate results within 0.5 kcal mol^-1 of coupled cluster singles and doubles results. On a large cluster of water molecules, we achieve calculations involving over 3500 orbital and 12,000 auxiliary basis functions in under 10 min on a single CPU core.

Benchmarking the Effective Fragment Potential Dispersion Correction on the S22 Test Set

The usual modeling of dispersion interactions in density functional theory (DFT) is often limited by the use of empirically fitted parameters. In this study, the accuracies of the popular empirical dispersion corrections and the first-principles derived effective fragment potential (EFP) dispersion correction are compared by computing the DFT-D and HF-D equilibria interaction energies and intermolecular distances of the S22 test set dimers. Functionals based on the local density approximation (LDA) and generalized gradient approximation (GGA), as well as hybrid functionals, are compared for the DFT-D calculations using coupled cluster CCSD(T) at the complete basis set (CBS) limit as the reference method. In general, the HF-D(EFP) method provides accurate equilibrium dimerization energies and intermolecular distances for hydrogen-bonded systems compared to the CCSD(T)/CBS reference data without using any empirical parameters. For dispersion-dominant and mixed systems, the structures and interaction energies obtained with the B3LYP-D(EFP) method are similar to or better than those obtained with the other DFT-D and HF-D methods. Overall, the first-principles derived -D(EFP) correction presents a robust alternative to the empirical -D corrections when used with the B3LYP functional for dispersion-dominant and mixed systems or with Hartree-Fock for hydrogen-bonded systems.

Intermolecular interactions of oligothienoacenes: Do S⋯S interactions positively contribute to crystal structures of sulfur-containing aromatic molecules?

Intermolecular interactions in the crystals of tetra- and penta-thienoacene were studied using ab initio molecular orbital calculations for evaluating the magnitude of characteristic S⋯S interactions with great attention paid to their origin. The interactions between the π-stacked neighboring molecules are significantly greater than those between the neighboring molecules exhibiting the S⋯S contact, although it has sometimes been claimed that the S⋯S interactions play important roles in adjusting the molecular arrangement of sulfur-containing polycyclic aromatic molecules in the crystals owing to short S⋯S contacts. The coupled cluster calculations with single and double substitutions with noniterative triple excitation interaction energies at the basis set limit estimated for the π-stacked and S⋯S contacted neighboring molecules in the tetrathienoacene crystal are -11.17 and -4.27 kcal/mol, respectively. Those for π-stacked molecules in the pentathienoacene crystal is -14.38 kcal/mol, while those for S⋯S contacted molecules are -7.02 and -6.74 kcal/mol. The dispersion interaction is the major source of the attraction between the π-stacked and S⋯S contacted molecules, while the orbital-orbital interactions are repulsive: The orbital-orbital interactions, which are significant for charge carrier transport properties, are not much more than the results of the short S⋯S contact caused by the strong dispersion interactions. Besides, the intermolecular interaction energy calculated for a trithienoacene dimer has strong orientation dependence.

Self-Assembled Dehydro[24]annulene Monolayers at the Liquid/Solid Interface: Toward On-Surface Synthesis of Tubular π-Conjugated Nanowires

We have studied the self-assembly behavior of dehydro[24]annulene (D24A) derivatives 1, 2a-2d, and 3a-3c at the liquid/solid interface using scanning tunneling microscopy (STM). Both the relative placement and the nature of the four D24A substituents strongly influence the self-assembly pattern. Overall, the eight D24A derivatives examined in this study display seven types of 2D packing patterns. The D24A derivatives 1, 2a, and 3a have either two or four stearate groups and adopt face-on configurations of their macrocyclic cores with respect to the highly oriented pyrolytic graphite (HOPG) surface. Their 2D packing pattern is determined by the interchain spacings and number of stearate substituents. The D24A derivatives 2b-2d and 3b-3c bear hydrogen-bonding carbamate groups to further strengthen intermolecular interactions. Face-on patterns were also observed for most of these compounds, while an unstable edge-on self-assembly was observed in the case of 2b at room temperature. Stable edge-on self-assemblies of D24A derivatives were sought for this work as an important stepping stone to achieving the on-surface topochemical polymerization of these carbon-rich macrocycles into tubular π-conjugated nanowires. The overall factors determining the 2D packing patterns of D24As at the liquid/solid interface are discussed on the basis of theoretical simulations, providing useful guidelines for controlling the self-assembly pattern of future D24A macrocycles.

Big Changes for Small Noncovalent Dimers: Revisiting the Potential Energy Surfaces of (P2)2 and (PCCP)2 with CCSD(T) Optimizations and Vibrational Frequencies

This article details the re-examination of low-lying stationary points on the potential energy surfaces (PESs) of two challenging noncovalent homogeneous dimers, (P2)2 and (PCCP)2. The work was motivated by the rather large differences between MP2 and CCSD(T) energetics that were recently reported for these systems (J. Comput. Chem. 2014, 35, 479-487). The current investigation reveals significant qualitative and quantitative changes when the CCSD(T) method is used to characterize the stationary points instead of MP2. For example, CCSD(T) optimizations and harmonic vibrational frequency computations with the aug-cc-pVTZ basis set indicate that the parallel-slipped (PS) structure is the only P2 dimer stationary point examined that is a minimum (zero imaginary frequencies, ni = 0), whereas prior MP2 computations indicated that it was a transition state (ni = 1). Furthermore, the L-shaped structure of (P2)2 was the only minimum according to MP2 computations, but it collapses to the PS structure on the CCSD(T)/aug-cc-pVTZ PES. For the larger PCCP dimer, the CCSD(T) computations reveal that four rather than just two of the six stationary points characterized are minima. A series of explicitly correlated single-point energies were computed for all of the optimized structures to estimate the MP2 and CCSD(T) electronic energies at the complete basis set limit. CCSDT(Q) computations were also performed to assess the effects of dynamical electron correlation beyond the CCSD(T) level. For both (P2)2 and (PCCP)2, dispersion remains the dominant attractive component to the interaction energy according to symmetry-adapted perturbation theory analyses, and it is also the most challenging component to accurately evaluate.

Benzene Dimer: High-Level Wave Function and Density Functional Theory Calculations

High-level OVOS (optimized virtual orbital space) CCSD(T) interaction energy calculations (up to the aug-cc-pVQZ basis set) and various extrapolations toward the complete basis set (CBS) limit are presented for the most important structures on the benzene dimer potential energy surface. The geometries of these structures were obtained via an all-coordinate gradient geometry optimization using the DFT-D/BLYP method, covering the empirical dispersion correction fitted exclusively for this system. The fit was carried out against two estimated CCSD(T)/CBS potential energy curves corresponding to the distance variation between two benzene rings for the parallel-displaced (PD) and T-shaped (T) structures. The effect of the connected quadruple excitations on the interaction energy was estimated using the CCSD(TQf) method in a 6-31G*(0.25) basis set, destabilizing the T and T-shaped tilted (TT) structures by ≈0.02 kcal/mol and the PD structure by ≈0.04 kcal/mol. Our best CCSD(T)/CBS results show, within the error bars of the applied methodology, that the energetically lowest-lying structure is the TT structure, which is nearly 0.1 kcal/mol more stable than the almost isoenergetic PD and T structures. The specifically parametrized DFT-D/BLYP method leads to a correct energy ordering of the structures, with the errors being smaller by 0.2 kcal/mol with respect to the most accurate CCSD(T) values.

Appointing silver and bronze standards for noncovalent interactions: a comparison of spin-component-scaled (SCS), explicitly correlated (F12), and specialized wavefunction approaches

A systematic examination of noncovalent interactions as modeled by wavefunction theory is presented in comparison to gold-standard quality benchmarks available for 345 interaction energies of 49 bimolecular complexes. Quantum chemical techniques examined include spin-component-scaling (SCS) variations on second-order perturbation theory (MP2) [SCS, SCS(N), SCS(MI)] and coupled cluster singles and doubles (CCSD) [SCS, SCS(MI)]; also, method combinations designed to improve dispersion contacts [DW-MP2, MP2C, MP2.5, DW-CCSD(T)-F12]; where available, explicitly correlated (F12) counterparts are also considered. Dunning basis sets augmented by diffuse functions are employed for all accessible ζ-levels; truncations of the diffuse space are also considered. After examination of both accuracy and performance for 394 model chemistries, SCS(MI)-MP2/cc-pVQZ can be recommended for general use, having good accuracy at low cost and no ill-effects such as imbalance between hydrogen-bonding and dispersion-dominated systems or non-parallelity across dissociation curves. Moreover, when benchmarking accuracy is desirable but gold-standard computations are unaffordable, this work recommends silver-standard [DW-CCSD(T**)-F12/aug-cc-pVDZ] and bronze-standard [MP2C-F12/aug-cc-pVDZ] model chemistries, which support accuracies of 0.05 and 0.16 kcal/mol and efficiencies of 97.3 and 5.5 h for adenine·thymine, respectively. Choice comparisons of wavefunction results with the best symmetry-adapted perturbation theory [T. M. Parker, L. A. Burns, R. M. Parrish, A. G. Ryno, and C. D. Sherrill, J. Chem. Phys. 140, 094106 (2014)] and density functional theory [L. A. Burns, Á. Vázquez-Mayagoitia, B. G. Sumpter, and C. D. Sherrill, J. Chem. Phys. 134, 084107 (2011)] methods previously studied for these databases are provided for readers' guidance.

A theoretical benchmark study of the spectroscopic constants of the very heavy rare gas dimers

The binding energy of the superheavy dimer Uuo₂ is considerably larger than that of its lighter homologues, despite a 40% reduction due to spin-other orbit interaction.

Evaluation of composite schemes for CCSDT(Q) calculations of interaction energies of noncovalent complexes

Evaluation of composite schemes for CCSDT(Q) calculations of interaction energies of noncovalent complexes.

Examination of tyrosine/adenine stacking interactions in protein complexes

The π-stacking interactions between tyrosine amino acid side chains and adenine-bearing ligands are examined. Crystalline protein structures from the protein data bank (PDB) exhibiting face-to-face tyrosine/adenine arrangements were used to construct 20 unique 4-methylphenol/N9-methyladenine (p-cresol/9MeA) model systems. Full geometry optimization of the 20 crystal structures with the M06-2X density functional theory method identified 11 unique low-energy conformations. CCSD(T) complete basis set (CBS) limit interaction energies were estimated for all of the structures to determine the magnitude of the interaction between the two ring systems. CCSD(T) computations with double-ζ basis sets (e.g., 6-31G*(0.25) and aug-cc-pVDZ) indicate that the MP2 method overbinds by as much as 3.07 kcal mol(-1) for the crystal structures and 3.90 kcal mol(-1) for the optimized structures. In the 20 crystal structures, the estimated CCSD(T) CBS limit interaction energy ranges from -4.00 to -6.83 kcal mol(-1), with an average interaction energy of -5.47 kcal mol(-1), values remarkably similar to the corresponding data for phenylalanine/adenine stacking interactions. Geometry optimization significantly increases the interaction energies of the p-cresol/9MeA model systems. The average estimated CCSD(T) CBS limit interaction energy of the 11 optimized structures is 3.23 kcal mol(-1) larger than that for the 20 crystal structures.

Achieving the CCSD(T) Basis-Set Limit in Sizable Molecular Clusters: Counterpoise Corrections for the Many-Body Expansion

An efficient procedure is introduced to obtain the basis-set limit in electronic structure calculations of large molecular and ionic clusters. This approach is based on a Boys-Bernardi-style counterpoise correction for clusters containing arbitrarily many monomer units, which is rendered computationally feasible by means of a truncated many-body expansion. This affords a tractable way to apply the sequence of correlation-consistent basis sets (aug-cc-pVXZ) to large systems and thereby obtain energies extrapolated to the complete basis set (CBS) limit. A three-body expansion with three-body counterpoise corrections is shown to afford errors of ≲0.1-0.2 kcal/mol with respect to traditional MP2/CBS results, even for challenging systems such as fluoride-water clusters. A triples correction, δCCSD(T) = ECCSD(T) - EMP2, can be estimated accurately and efficiently as well. Because the procedure is embarrassingly parallelizable and requires no electronic structure calculations in systems larger than trimers, it is extendible to very large clusters. As compared to traditional CBS extrapolations, computational time is dramatically reduced even without parallelization.

Basis set dependence of higher-order correlation effects in π-type interactions

The basis set dependence of higher-order correlation effects on π-type interaction energies was examined by scanning the potential energy surfaces of five dimer systems. The dimers of acetylene (H-C≡C-H), diacetylene (H-C≡C-C≡C-H), cyanogen (N≡C-C≡N), diphosphorous (P≡P), and 1,4-diphosphabutadiyne (P≡C-C≡P) were studied in three different configurations: cross, parallel-displaced, and t-shaped. More than 800 potential energy curves (PECs) were generated by computing the interaction energies for all 15 dimer configurations over a range of intermolecular distances with the MP2, coupled-cluster single double (CCSD), and coupled-cluster single double triple (CCSD(T)) methods in conjunction with 21 basis sets ranging from a small 6-31G*(0.25) split-valence basis set to a large aug-cc-pVQZ correlation consistent basis set. Standard extrapolation techniques were also used to construct MP2, CCSD, and CCSD(T) complete basis set (CBS) limit PECs as well as CBS limit higher-order correlation corrections based on the differences between CCSD(T) and MP2 interaction energies, denoted δ(MP2)(CCSD(T)), and the corresponding differences between CCSD(T) and CCSD interactions energies, denoted δ(CCSD)(CCSD(T)). Double-ζ basis sets struggled to reproduce the former but provided quite reasonable descriptions of the latter as long as diffuse functions were included. The aug-cc-pVDZ basis deviated from the δ(CCSD)(CCSD(T)) CBS limit by only 0.06 kcal mol(-1) on average and never by more than 0.24 kcal mol(-1), whereas the corresponding deviations were approximately twice that for the δ(MP2)(CCSD(T)) term. While triple-ζ basis sets typically improved results, only aug-cc-pVTZ provided appreciable improvement over utilizing the aug-cc-pVDZ basis set to compute δ(CCSD)(CCSD(T)). Counterpoise (CP) corrections were also applied to all double- and triple-ζ basis sets, but they rarely yielded a better description of these higher-order correlation effects. CP corrections only consistently improved results when the aug-cc-pVDZ basis set was used to compute δ(MP2)(CCSD(T)), yielding mean and maximum absolute deviations from the CBS values of 0.10 and 0.39 kcal mol(-1), respectively, for all five dimer systems.

Accurate Interaction Energies for Problematic Dispersion-Bound Complexes: Homogeneous Dimers of NCCN, P2, and PCCP

All intermolecular interactions involve London dispersion forces. The accurate treatment of dispersion is essential for the computation of realistic interaction potentials. In general, the most reliable method for computing intermolecular interactions is coupled-cluster singles and doubles with perturbative triples [CCSD(T)] in conjunction with a sufficiently flexible Gaussian atomic orbital basis set, a combination which is not routinely applicable due to its excessive computational demands (CPU time, memory, storage). Recently, many theoretical methods have been developed that attempt to account for dispersion in a more efficient manner. It is well-known that dispersion interactions are more difficult to compute in some systems than others; for example, π-π dispersion is notoriously difficult, while alkane-alkane dispersion is relatively simple to compute. In this work, numerous theoretical methods are tested for their ability to compute reliable interaction energies in particularly challenging systems, namely, the P2, PCCP, and NCCN dimers. Symmetry-adapted perturbation theory (SAPT) is applied to these dimers to demonstrate their sensitivity to the treatment of dispersion. Due to the small size of these systems, highly accurate CCSD(T) potential energy curves could be estimated at the complete basis set limit. Numerous theoretical methods are tested against the reliable CCSD(T) benchmarks. Methods using a treatment of dispersion that relies on time-dependent density functional theory (TDDFT) response functions are found to be the most reliable.

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

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