they provide rigid building blocks for anchoring 
 substituents in defined geometric positions , thereby determining the 
 overall molecular shape . their unique 
 electronic structure with a distinct -cloud located parallel 
 above and below the ring plane allows for a variety of interaction 
 patterns . the distribution as well as 
 the overall density of electrons within the conjugated -system 
 of a prototypical benzene molecule can either be modulated by substituents 
 in various geometric arrangements or by exchanging 
 ring carbon atoms by nitrogen , oxygen , or sulfur creating heterocyclic 
 systems . their unique interactions are 
 a reason for their frequent occurrence in natural as well as synthetic 
 bioactive molecules . interaction patterns exhibited by aromatic 
 heterocycles comprise 
 hydrophobic , polar , hydrogen bonding , cation- , amid- , halogen- , and -stacking interactions . -stacking interactions have been investigated 
 in detail by a wide range of experimental and computational methods : 
 for benzene dimers , in dna or in proteins . however , it should be noted that the study of noncovalent interactions 
 remains challenging for computational as well as experimental approaches . in addition to the intrinsic 
 relevance of this specific system , these calculations are often used 
 to evaluate the performance of quantum chemical and force field methods 
 in describing aromatic systems . previous studies used 
 coupled cluster with singles , doubles , and perturbative triples ( ccsd(t ) ) up to complete basis set ( cbs ) extrapolation to evaluate accurate interaction energies . plesset perturbation theory was found inadequate 
 to describe the interaction energies of -stacking complexes , 
 as it tends to significantly overestimate interaction energies and 
 underestimate equilibrium distances at the cbs limit . compared to the benzene dimer , studies of interactions 
 of aromatic heterocycles are far less numerous . some attention has 
 been paid to the interaction of benzene with pyridine and the pyridine 
 dimer . as the main focus of this study 
 is to characterize the potential energy surfaces ( pes ) for a diverse 
 set of 13 aromatic heterocycles , the coupled cluster approach would 
 be prohibitively expensive . nevertheless , we employed ccsd(t ) calculations 
 using cbs extrapolation and the cc - pvtz basis set as a reference to identify a suitable density functional 
 method for our calculations . density functional theory ( dft ) has 
 been the most widely used approach in electronic structure calculations , 
 providing a reasonable trade - off between accuracy and computational 
 efficiency . since long - range dispersion effects are known to be an 
 important contribution to the interaction energy of the benzene dimer , 
 computationally cheap hartree recently , several density functionals including a correction 
 term for long - range dispersion interactions have been developed . the most common functional form for a correction 
 term is an r potential using a nucleus - specific 
 dispersion contribution for each atom . these parameters are combined for each pair of atoms to give a dispersion 
 parameter for the respective interaction . furthermore , an exponential 
 damping function is usually applied to avoid unphysical attraction 
 at short - range . some studies suggest an overestimation of the stability 
 of large  complexes in pairwise long - range corrected 
 dft methods . they propose the inclusion of 
 an empirical 3-body dispersion term . however , our 
 reference calculations at the ccsd(t)/cbs level of theory suggest 
 that this effect is not relevant for the size of the systems investigated 
 in this study . interestingly , only a limited set of aromatic 
 heterocycles frequently 
 appears in drug molecules . we performed extensive dft calculations 
 to characterize -stacking interactions between benzene and 
 a variety of aromatic heterocycles relevant in drug design as identified 
 by ertl and co - workers . ligand 
 recognition processes and provide guidelines for ligand optimization 
 along the lines of bissantz et al . a variety of interaction geometries deviating from the optimal 
 relative orientation for heterocyclic - aromatic interactions can be 
 observed in crystal structures ( figure 1 ) . observed 
 interactions of aromatic heterocycles in ligands with 
 aromatic protein side chains in the pdb extracted with relibase+ . the two distinct maxima represent a parallel - displaced ( low interplanar 
 angle , left peak ) as well as a t - shaped interaction geometry ( high 
 interplanar angle , right peak ) . for the purposes of this study , furthermore , 
 the geometric requirements of a certain binding mode do not allow 
 for arbitrary placement of interaction partners . we provide fundamental 
 data useful in structure - based lead optimization by suggesting possible 
 ring substitutions that provide a gain in affinity through strengthened 
 -stacking interactions . a coordinate system was 
 defined for all complexes as outlined in 
 figure 2 . we set the origin to the centroid 
 of the benzene molecule , the z - axis as a normal on 
 the ring plane of the benzene molecule and the x - axis 
 from the origin through one of the c a second benzene molecule was placed in the exact same orientation 
 at an interplanar distance of 3.0  . the interplanar distance 
 of this complex was subsequently varied in steps of 0.1  . the 
 most favorable interaction energy along this path was defined as the 
  stacked  geometry ( point a , figure 2 ) and used as a starting point for all further geometric manipulations . procedure 
 for scanning parallel - displaced pes : after initial identification 
 of the stacked ( point a ) and parallel - displaced ( point b ) minima ( ring 
 centers denoted by + ) , a path x ( red ) for lateral 
 displacement is chosen through these points . starting from point a , 
 the aromatic heterocycle is moved along x in increments of 
 1/10 the distance ab . while moving the aromatic heterocycle 
 through 20 increments of x displacement beyond point b , it 
 is rotated through 360 in 30 intervals around a ring - plane 
 normal axis in its center of geometry z ( blue ) at each point . interaction energies were calculated 
 by subtracting the energy 
 of the two monomers , as suggested by the supermolecular approach : in order to identify 
 a suitable dft functional for this study , 
 the stacked pes of the benzene dimer was scanned by performing ab initio calculations using gaussian09 rev . a.02 at the ccsd(t)/cbs level of theory by using 
 aug - cc - pvdz , aug - cc - pvtz , and aug - cc - pvqz to obtain the interaction 
 energy at the mp2/cbs limit through three - point extrapolation as described 
 by helgaker et al . the difference between 
 mp2/aug - cc - pvdz and ccsd(t)/aug - cc - pvdz was added to the extrapolated 
 mp2/cbs energy to obtain the ccsd(t)/cbs value . the monomer geometry 
 of benzene was obtained by using the default geometry in gaussview 
 4.1 . for further validation , lateral displacement along dx was varied 
 in order to sample interaction energies of parallel - displaced geometries 
 at the ccsd(t)/cc - pvtz level of theory . at each interplanar distance 
 dz and lateral displacement fock 
 method 
 and several density functional methods comprising b3lyp and the dispersion corrected methods m06 - 2x and b97xd using the cc - pvdz , aug - cc - pvdz , cc - pvtz and aug - cc - pvtz basis sets 
 from dunning and co - workers . results 
 obtained using b97xd / cc - pvtz were consistently within less 
 than 0.5 kcal / mol compared to ccsd(t)/cbs calculations for stacked 
 geometries ( see supporting information figure 
 1 ) and ccsd(t)/cc - pvtz calculations of parallel - displaced geometries 
 ( see supporting information figure 2 ) . 
 therefore , all subsequent interaction energy and dipole moment calculations 
 were performed at the b97xd / cc - pvtz level of theory . a set of seven 5-membered and six 6-membered heteroaromatic systems 
 ( see figure 3 ) was prepared by fully optimizing 
 each aromatic heterocycle as a monomer at the b97xd / cc - pvtz 
 level of theory . a complex with benzene was prepared from each resulting 
 geometry in the relative orientation outlined in figures 2 and 3 . to scan the pes of 
 these complexes , rotational and translational degrees of freedom initial ring 
 orientations : all investigated aromatic heterocycles , 
 in their z = 0 orientation and the direction of rotation 
 along z are indicated in this plot . the parallel - displaced minimum ( point b , figure 2 ) was identified by a sequence of coarse manual 
 optimization 
 followed by free geometry optimization . manual optimization was performed 
 by displacing the aromatic heterocycle from the stacked minimum along 
 the x direction in steps of 0.1  while keeping the benzene at 
 the origin stationary . at each point , the aromatic heterocycle was 
 rotated to sequentially align each ring atom with the x - axis in order to identify the preferred orientation . subsequently , 
 a free geometry optimization was performed starting from the energetically 
 most favorable orientation identified manually . the aromatic heterocycle 
 centroid of the resulting geometry was defined as point b. after 
 identifying both the respective stacked and parallel - displaced 
 minimum geometries , a straight transition path x through these 
 points was defined ( figure 2 ) . each aromatic 
 heterocycle was moved along this path in increments of 1/10 of the 
 distance ab for a total of 20 increments while keeping the 
 ring planes parallel . at each increment of x displacement , 
 the aromatic heterocycle was rotated in 12 steps of 30 in the 
 counterclockwise direction . single point energies were calculated 
 for each geometry in order to sample the pes starting from the stacked 
 configuration in point a through the parallel - displaced minimum in 
 point b and beyond . location of points a and b as well as the slope 
 for x for all aromatic heterocycles are given in supporting information table 2 . as dx denotes lateral 
 displacement , it is 0 for the stacked geometry and is given only for 
 the parallel - displaced minimum at point b. the slope denotes the change 
 in dz per 1  of lateral displacement along dx when following 
 the transition path x. using this approach , we can scan the 
 most relevant region of the respective pes of parallel - displaced arrangements 
 using just 2 parameters . all calculated single points for each 
 system were combined to form 
 a pes for parallel - displaced geometries of the respective system contour plots of the pes for 
 all systems near their respective parallel - displaced minima are depicted 
 in figure 4 . experimental data was extracted 
 from the pdb using relibase+ ( version 
 3.1 ) . for the analysis of geometries , 
 all aromatic 5- and 6-rings within 
 a distance of less than 6.5  of any aromatic side chain atom 
 of a receptor were selected . heme and the purine and pyrimidine dna 
 and rna bases were excluded from this selection in order to focus 
 the analysis on drug - like molecules . furthermore , all nmr structures 
 as well as all x - ray structures with a resolution worse than 2.5  
 were not considered in this analysis . in order to validate our theoretical 
 findings the tested 
 density functionals show significant differences in 
 their ability to accurately describe -stacking interactions : 
 b3lyp did not show any discernible favorable interaction between two 
 benzene molecules . the m06 - 2x functional exhibited oscillations in 
 interaction energy in areas away from the minimum . this was deemed 
 undesirable , since the study aims at characterizing large areas of 
 the pes in order to estimate preferred as well as disfavored interaction 
 geometries for the respective aromatic heterocycles . whereas 
 b97xd with both double- basis sets gave unsatisfactory 
 results with energy deviations exceeding 1 kcal / mol compared to ccsd(t)/cbs , 
 the triple- basis sets were able to yield consistent interaction 
 energies to ccsd(t)/cbs with deviations smaller than 0.5 kcal / mol . usage of aug - cc - pvtz did not yield significant improvements in accuracy 
 over the cc - pvtz basis set . therefore , cc - pvtz was used for all subsequent 
 calculations in order to optimize computational efficiency . pes were 
 calculated as outlined in the methods section 
 for all 12 heteroaromatic complexes and the benzene dimer . the resulting 
 energy landscapes at a distance of up to 2.5  in the x 
 direction are given in figure 4 . 3d depictions 
 of the parallel - displaced minima as identified through the pes scans 
 are given in figure 5 . energy surfaces 
 of parallel - displaced geometries for 6-membered ( a benzene , b pyridine , 
 c pyridazine , d pyrimidine , e pyrazine , f s - triazine ) and 5-membered 
 ( g furan , h thiophene , i oxazole , j thiazole , k isoxazole , l isothiazole , 
 m 1,3,4-thiadiazole ) rings investigated in this study are given here . a common coloring scheme across all plots gives an indication of the 
 relative energetics of the interactions . 3d depictions of the minimum geometries identified on the pes for 
 all investigated systems . ( a benzene , b pyridine , c pyridazine , d 
 pyrimidine , e pyrazine , f s - triazine , g furan , h thiophene , i oxazole , 
 j thiazole , k isoxazole , l isothiazole , m 1,3,4-thiadiazole ) . placing an electronegative atom 
 over the benzene ring surface is 
 energetically disfavored in all cases . this also holds true for multiple 
 heteroatoms in a single ring , e.g. in the case of pyrazine , where 
 two distinct saddle points in the rotational dimension can be observed 
 at z angles of 0 and 180 , respectively ( figure 4e ) . the preferred rotational orientation is highly 
 dependent on the type of aromatic heterocycle as shown by the pes 
 in figure 4 . a consistent interaction pattern 
 can be observed for all investigated aromatic heterocycles : all systems 
 exhibit a pronounced minimum at a parallel - displaced geometry at dx 
 of approximately 1.5  as outlined in figure 4 . interaction energies and minimum geometries 
 resulting from unrestricted optimization are given in figure 6 . interaction energies for all investigated heterocycles 
 at the parallel - displaced 
 minimum : the x - axis corresponds to the centroid distance 
 at the energy minimum . rotation of heterocycles is consistent with 
 the minimum geometry analogous to the benzenepyridine example . the interaction minima for aromatic 
 heterocycles in parallel - displaced 
 geometries are up to 1.5 kcal / mol more favorable than for the benzene 
 dimer . unfavorable orientations carry an energy penalty of up to 2 
 kcal / mol over benzene benzene even in close proximity ( displacement 
 below 0.5  ) to the minimum geometry . nitrogen is generally oriented 
 away from the ring center of the interaction partner in the respective 
 minimum geometry . a closer look at the series of 1,2- , 1,3- , 
 and 1,4-diazabenzenes 
 ( pyridazine , pyrimidine , and pyrazine ) reveals an interesting trend : 
 in the case of pyrazine , it is geometrically impossible for both nitrogens 
 to point away from the ring center of the interacting benzene due 
 to the 1,4 positioning of the nitrogen atoms . it also has the weakest 
 interaction energy at the minimum geometry of this series . pyrimidine 
 has a strong preferential orientation with both nitrogen atoms pointing 
 away from the benzene ring center . furthermore , 
 pyridazine disfavors an orientation with the nitrogen atoms located 
 directly above the benzene ring most strongly of all molecules in 
 this series . in the series oxazole to thiazole or isoxazole 
 to isothiazole no 
 significant changes of the interaction energies of these systems with 
 benzene are observed . however , due to the larger size of the sulfur 
 atom , the interplanar distance at the minimum is increased by approximately 
 0.1  in the case of isothiazole and by 0.05  in the case 
 of thiazole compared to isoxazole and oxazole , respectively . still , 
 replacing furan by thiophene increases the interaction energy slightly 
 and also increases the equilibrium distance by approximately 0.05 
  . if we consider only systems with nonzero dipole moment 
 ( excluding 
 benzene , pyrazine , and s - triazine ) , the resulting correlation of the 
 magnitude of the dipole moment with the minimum interaction energy 
 in parallel - displaced geometry is striking . the squared correlation 
 coefficient is 0.78 as outlined in figure 7 . in general , a significant portion 
 of the observed variance in interaction energy at the parallel - displaced 
 minimum geometry can be explained by the magnitude of the dipole moment . all systems with a dipole moment of zero do have nonzero quadrupole 
 moments by which they can show electrostatic interactions as well . the observed deviation of pyrimidine from this trend is likely caused 
 by the partial compensation of dipole moment vectors induced by the 
 1,3-positioning of the nitrogen atoms . this will lead to higher electrostatic 
 moments being more pronounced , explaining the high interaction energy 
 of pyrimidine . correlation of interaction energy of 10 aromatic heterocycles 
 and 
 benzene at the minimum geometry with dipole moment : all systems without 
 dipole moment were excluded from this graph ( benzene , pyrazine and 
 s - triazine ) . ( left to right : thiophene , furan , thiazole , oxazole , pyridine , pyrimidine , 
 isothiazole , isoxazole , 1,3,4-thiadiazole ) almost 80% of the variance 
 in interaction energy can be explained by the magnitude of the dipole 
 moment alone . a higher dipole moment (  , debye ) yields a more 
 favorable interaction energy ( ie , kcal / mol ) at the minimum . as expected , the benzene dimer 
 ( figure 4a ) shows no preferred orientation . 
 while this is also true for s - triazine 
 ( figure 4f ) , s - triazine is the only investigated 
 system featuring a low - energy stacked geometry . pyridine ( figure 4b ) prefers any parallel - displaced orientation locating 
 the nitrogen atom away from the benzene -cloud . this is also 
 true for pyridazine ( figure 4c ) , showing a 
 stronger gain in interaction energy than pyridine for favored configurations . 
 whereas pyrimidine ( figure 4d ) also prefers 
 to orient its nitrogen atoms away from the benzene -cloud , 
 pyrazine ( figure 4e ) does not allow this configuration 
 due to the 1,4 positioning of its nitrogen atoms . pyrazine prefers 
 to have its nitrogen atoms oriented perpendicular to the axis of displacement 
 over one nitrogen atom placed above and the other outside the -cloud 
 of benzene . furan ( figure 4 g ) and thiophene 
 ( figure 4h ) show very similar interaction patterns 
 in terms 
 of energetics as well as orientational preference . both systems favor 
 placing the heteroatom above the hydrogen atoms of benzene located 
 at a 6090 angle of the axis of displacement . due to 
 the larger size of sulfur , the location of the minimum in thiophene 
 is at a slightly ( 0.2  dx ) larger distance than for 
 furan . the similarity of oxo- and thio - substituted aromatic heterocycles 
 holds also true for oxazole ( figure 4i ) and 
 thiazole ( figure 4j ) as well as isoxazole ( figure 4k ) and isothiazole ( figure 4l ) . in all cases , the positioning of the nitrogen atom is again the 
 crucial determinant of preferred orientation . placing either oxygen 
 or sulfur directly above the center of the benzene molecule does not 
 lower the interaction energy . the case of 1,3,4-thiadiazole ( figure 4 m ) is especially interesting when compared to isothiazole 
 or isoxazole , as 1,3,4-thiadiazole exhibits the most favorable interaction 
 energy at an orientation that is an energetic maximum for both isoxazole 
 and isothiazole . from the presented data we deduce two main 
 contributing factors 
 to the energy and geometric preference of heteroaromatic stacking 
 interactions . this is exemplified by the benzene dimer , 
 which has no dipole moment , exhibiting only weakly polarized c h 
 bonds , yet showing a favorable interaction of larger than 2 kcal / mol 
 at the parallel - displaced minimum . dispersion interactions do not 
 show a significant directional preference and almost exclusively depend 
 on the approximately constant contact area of the involved species . in addition to the dispersion component , we find an equally important 
 electrostatic contribution . we presume that the geometric preference 
 of specific systems is mainly driven by these electrostatic interactions . this is consistent with the findings of sherrill in his recent study 
 of stacking interactions of benzene and pyridine homo- and heterodimers . we find that the minimum interaction energy 
 is strongly dependent 
 on the magnitude of the dipole moment . correlating the minimum interaction 
 energy with the magnitude of the dipole moment of a heterocyclic monomer 
 revealed an unexpectedly strong relationship between these factors . 
 around 80% of the observed variance in minimum interaction energy this finding is especially 
 relevant , since the most commonly occurring natural aromatic heterocycles , 
 i.e. pyrimidine and purine in the nucleobases , have a nonzero dipole 
 moment . exploiting electrostatic interaction 
 patterns might facilitate the design of environment - specific stacking 
 interactions by rational choice of specific aromatic heterocycles . consistent series of molecules can be identified wherein variation 
 of the placement of heteroatoms causes a systematic shift in the preferred 
 interaction geometry and strength : for one , this is exemplified by 
 a series of nitrogen - containing aromatic heterocycles : pyridine - pyridazine - pyrimidine - pyrazine . 
 throughout the series , interaction energy increases steadily . second , 
 substitution of ring atoms with atoms from the same group of the periodic 
 table cause differences in centroid - centroid distances consistent 
 with the size differences of these atoms while not affecting interaction 
 energies significantly . such substitutions are also demonstrated for 
 the cases of oxazole / thiazole and isoxazole / isothiazole , respectively . most investigated systems show broad energy minima 
 that enable a certain degree of flexibility in relative orientation . however , certain orientations should be avoided , as improper positioning 
 of electronegative atoms is associated with a significant energy penalty . pyridine for example allows a rotation of around 240 without 
 any significant deviation form the minimum interaction energy . the 
 remaining 120 of possible relative orientations with nitrogen 
 directly above the benzene ring can reduce the interaction energy 
 by more than 1 kcal / mol , thereby making the interaction less favorable 
 than an interaction with benzene . this effect is even more pronounced 
 in pyridazine with two adjacent nitrogen atoms , where the energy penalty 
 for improper relative rotation is more than 2 kcal / mol . interaction 
 geometries of pyridyl substrucures with phenylalanine 
 were identified using relibase+ as outlined in the methods section and subsequently mapped to our coordinate 
 system by calculating a best - fit plane through the centroid of the 
 phenyl substructure of phenylalanine , calculating the projection of 
 the ring centroid of the pyridyl substructure of the ligand to this 
 plane and the respective ring rotation z . it can be seen , that 
 the majority of interactions occur in regions identified as energetically 
 highly favorable ( complexes 1tka , 1lc8 , 1ioe , 1iqj , 2yri , 1iqg ) . however , 
 certain complexes present the pyridyl substructure to the receptor 
 in a nonoptimal way ( complexes 2r5c , 2fxd ) . these examples illustrate 
 the necessity for a more detailed picture of interaction energies 
 in nonminimum regions . our calculations suggest that reorientation 
 or replacement of the pyridyl moiety in 2fxd would provide a gain 
 in affinity . 1iqj is a co - crystal of coagulation factor xa with the inhibitor m55124 ( figure 8) . it is well - known that the hydrophobic s4 pocket of factor 
 xa interacts 
 with aromatic substructures . the depicted 
 example shows excellent -stacking with pyridine . according 
 to our analysis , the interaction geometry of the pyridyl moiety in 
 m55124 is in the optimal orientation relative to the side chain phe-174 
 in terms of its stacking interaction ( ie = 3.92 kcal / mol ) . however , we can not exclude the possibility for the pyridine to be 
 protonated , which would add a cation- component to this interaction . 2yri is an x - ray structure of alanine - pyruvate aminotransferase 
 in complex with an inhibitor that exhibits a substituted pyridine 
 as its central structural feature ( figure 8) . based on our analysis of interaction energies , the core pyridine 
 is ideally positioned relative to phe-85 in terms of -stacking 
 energy . the broad interaction energy minimum allows for this moiety 
 to be rotated at around 60 along z relative to the factor 
 xa - m55124 interaction described previously and still be optimal in 
 terms of energy ( ie = 3.97 kcal / mol ) . 2fxd is a crystal structure of hiv-1 
 protease in complex with the inhibitor atazanavir ( figure 8) . atazanavir contains a benzene moiety , which is 
 directly linked to the analyzed pyridine ring . observed stacking interactions 
 of the pyridyl moiety with phe-82 show a distinct deviation from identified 
 optimal relative geometries . according to our calculations , this interaction 
 geometry is more than 1 kcal / mol unfavorable over optimal phe - pyridine 
 stacking ( ie = 2.84 kcal / mol ) . however , the neighboring benzene 
 in atazanavir is also involved in stacking with phe-82 . we argue that 
 rotating the pyridine by 120 for the nitrogen to be in para - position 
 would improve stacking interactions . however , this might further reduce 
 coplanarity of the benzene and pyridine rings in atazanavir , as the 
 added hydrogen in the ortho - position would induce additional strain 
 energy in the ligand . this reduction in coplanarity other marketed 
 drugs that target hiv-1 protease frequently contain a single benzene 
 ring in an analogous position ( e.g. , ritonavir , lopinavir ) . occurrences of parallel - displaced phe - pyridine interactions 
 in 
 protein  ligand complexes identified in the pdb are overlaid 
 with the corresponding benzenepyridine pes . most interaction 
 geometries occur in a relative orientation identified as favorable , 
 as exemplified by the depicted interactions in 1iqj and 2yri . however , some interaction 
 geometries are clearly in a nonpreferred region as exemplified by 
 2fxd . we expect substituent effects to have a significant impact , both in 
 terms of geometric preferences and interaction energies . substituent 
 effects are known to change electrostatic properties of aromatic ring 
 systems ( e.g. , the inversion of electrostatics for hexafluorobenzene 
 at the extreme ) . although a detailed 
 investigation of substituent effects is beyond the scope of this study , 
 we assume the general workflow presented here to be applicable to 
 the investigation of substituent effects as well . moreover as the ligand and the protein are solvated in aqueous solution 
 prior to complex formation , ( partial ) desolvation of both interaction 
 partners occurs during binding . we assume that different solvation 
 energies of the investigated aromatic heterocycles have a significant 
 impact on the overall binding energetics . however , we do not expect 
 solvation to be crucial for the identified geometric preferences . in the course of this study we evaluated dispersion corrected density 
 functionals for their capacity to reproduce ccsd(t ) energy surfaces . we considered the b97xd and m06 - 2x functionals developed by 
 the head - gordon and truhlar groups , respectively . we found that both functionals 
 are well suited to reproduce energies and positions of the minima 
 of high - level post - hf calculations . we found b97xd to be better 
 suited for our purposes , as m06 - 2x exhibited oscillating energies 
 in geometries away from the local energy minima . energetic differences 
 between different interacting aromatic heterocycles 
 are not immediately obvious even to the trained chemist . therefore , 
 our study can provide a reference for structure - based drug design 
 by comparing ring orientations and their associated energies in a 
 comprehensive manner . the calculated pes allow to immediately identify 
 which ring replacements might produce a gain in affinity in cases 
 where -stacking interactions are relevant . the b97xd 
 functional employing high- basis sets is 
 well suited to investigate aromatic stacking interactions . we performed 
 extensive pes scans for 13 complexes of 5- and 6-membered aromatic 
 heterocycles with benzene . furthermore , we present 
 pes of parallel - displaced heteroaromatic systems as a guide to identify 
 especially favorable or unfavorable interaction geometries . moreover , 
 we identified a strong dependence on electrostatic effects for the 
 preferred interaction geometry as well as the minimum energy . we conclude that dispersion is a significant contributor to the 
 overall interaction energy , which is augmented to a varying degree 
 by electrostatic effects . the presented pes as well as the identified 
 optimal interaction geometries can provide guidance for structure - based 
 lead optimization . this information can assist a medicinal chemist 
 in identifying an appropriate aromatic heterocycle for a given interaction 
 geometry ( compare bissantz et al . ) .
in 
 this study we investigate -stacking interactions of a 
 variety of aromatic heterocycles with benzene using dispersion corrected 
 density functional theory . 
 we calculate extensive potential energy 
 surfaces for parallel - displaced interaction geometries . 
 we find that 
 dispersion contributes significantly to the interaction energy and 
 is complemented by a varying degree of electrostatic interactions . 
 we identify geometric preferences and minimum interaction energies 
 for a set of 13 5- and 6-membered aromatic heterocycles frequently 
 encountered in small drug - like molecules . 
 we demonstrate that the 
 electrostatic properties of these systems are a key determinant for 
 their orientational preferences . 
 the results of this study can be 
 applied in lead optimization for the improvement of stacking interactions , 
 as it provides detailed energy landscapes for a wide range of coplanar 
 heteroaromatic geometries . 
 these energy landscapes can serve as a 
 guide for ring replacement in structure - based drug design .