text:
# From ab initio quantum chemistry to molecular dynamics: The delicate case of hydrogen bonding in ammonia

## I Introduction and Motivation

Wavefunction–based quantum chemistry methods as well as density functional theory (DFT) are now well established electronic structure techniques within Theoretical Chemistry TC . The former is traditionally used to investigate small finite systems, i.e. “gas phase problems”, whereas the latter had its first successes in solid state physics, i.e., in the framework of “periodic lattice problems”. There is also ample evidence that wavefunction–based quantum chemistry methods are able to predict properties of molecular systems up to virtually any desired accuracy as a result of systematically improving the correlation treatment together with the wavefunction representation. However, the numerical complexity, as e.g. measured by the computer time needed for such studies, explodes easily in praxi either before reaching a sufficient level of accuracy or before treating systems that are sufficiently large. DFT methods, on the other hand, are known to lead to astonishingly good predictions for a wide variety of systems at moderate cost – of course dramatic failures are also well documented in the literature. When describing molecular solids, associated liquids, or large biosystems, noncovalent interactions such as hydrogen bonds and van der Waals interactions might play an important role. The ideal approach would be to use post Hartree-Fock methods converged to the basis set limit, which turns out to be computationally prohibitive, especially when used in conjunction with molecular dynamics. While simple DFT methods are unable to correctly describe van der Waals interactions Rappe ; Gdanitz ; Hirao , hydrogen bonds are known to be well within the capabilities of DFT approaches once suitable gradient–corrected functionals are devised. A pioneering step in the early nineties was the investigation of hydrogen bonds in the water dimer kari1 ; sim92 ; cpwatern  and in small water clusters kari2  based on such functionals, which are of the generalised gradient approximation (GGA) type. Soon thereafter it was demonstrated that also the structure and dynamics of liquid bulk water can be described faithfully cpwatern . Now there is ample evidence that various properties of water can be computed quite reliably using GGA–type functionals, see e.g. Refs. cpwatern ; tuck ; Sprik ; SP ; SBP ; MTHP ; prr00 ; pbe-water ; pbe-water2 ; HCTH120 ; MTP ; krack ; IV  and Refs. Pais ; TBH ; Boyd ; tHCTH  for dimer calculations. But how about other hydrogen bonded systems?

Among hydrogen–bonded dimers the ammonia dimer is peculiar since it features a rather unusual hydrogen bond, see Refs. nelson87 ; havenith  for an early and a very recent overview. Instead of an essentially linear hydrogen bond, similar to the one of the water dimer having a H–O$`\mathrm{}`$O angle of only 5.5 klopper00 , the hydrogen bond in the ammonia dimer is significantly bent as depicted in Fig. 1c or d. Furthermore, close in energy is the so–called cyclic structure where an even more pronounced bending of the hydrogen bond produces two equivalent “hydrogen bonding contacts” according to Fig. 1e. As an interlude we stress that such strongly bent hydrogen bonds play an important role in many biological systems, in particular if the bending allows for more than one interaction of a donor with acceptor molecules which leads to so–called bifurcated or multifurcated hydrogen bonds steiner02 ; noteHB . In carbohydrates over 25% of the hydrogen bonds turn out to be multifurcated steiner02  and in proteins Saenger  about 40 and 90% of the hydrogen bonds in $`\beta `$–sheet and $`\alpha `$–helix structures, respectively, are of this type. In this context the ammonia dimer could serve as a prototype system in order to investigate nonlinear hydrogen bonding on the level of a minimal model, i.e. a free dimer in the gas phase.

Guided by analogy or chemical intuition, one might guess that the ammonia dimer possesses a “classical” quasi–linear hydrogen bond similar to other dimers such as those of water or hydrogen fluoride. Early calculations, see e.g. Refs. lin1 ; rlin1 ; rlin2  and references cited therein, actually assumed from the outset a perfectly linear bond similar to the one shown in Figs. 1a or b or predicted a quasi–linear hydrogen bond as depicted in Figs. 1c or d, while others C2h ; ahlrichs ; cyclic1  favoured strongly bent up to cyclic dimers similar to Figs. 1c–e. This controversy got fueled by microwave measurements MW  providing evidence for a quasi–rigid cyclic structure such as the one in Fig. 1e. The issue of a quasi–linear hydrogen bond vs. a cyclic structure remained controversial (see e.g. Refs. KT ; rlin3 ) up to early nineties where both experiment and theory started to converge. Thanks to sophisticated experiments hav91 ; hav92 ; loeser92 ; say-vrtrev ; DP ; heineking95 ; linnartz95 ; cotti96  and extensive computations Olthof ; Olthof2  there is now consensus havenith  that the ammonia dimer is a fluxional molecular complex with an equilibrium structure that is characterized by a bent hydrogen bond such as the one shown in Fig. 1c. Combining far–infrared spectra with static and dynamic calculations Olthof ; Olthof2 ; noteFIX , the H–N$`\mathrm{}`$N hydrogen bond angle (between the N–H bond of the proton donor and the N$`\mathrm{}`$N axis) was determined to be 26 in the equilibrium structure with a N$`\mathrm{}`$N distance of $`3.354`$ Å. Note, however, that this “experimental” structure of the ammonia dimer was deduced within the rigid monomer approximation Olthof ; noteFIX , which might introduce a systematic bias. In fact, a recent study saykally02  of the water dimer has shown that monomer flexibility is indeed important in determining its vibration–rotation–tunneling spectrum and thus “experimental” structures derived from it. It is expected that such monomer flexibility would also be important for the ammonia dimer. Furthermore, the global potential energy surface is found to be very flat Olthof ; Olthof2 , e.g. the cyclic C2h structure was estimated to be about only 7 cm-1 higher in energy, which is qualitatively consistent with the observation that the dimer is a very fluxional complex. Most interestingly, there is now experimental support nh3super  advocating a more cyclic equilibrium structure for (NH3)2 in superfluid Helium droplets at 0.4 K.

More recent wavefunction–based quantum chemical calculations LeePark ; Lindh  on the structure and energetics of the ammonia dimer support in some aspects the earlier conclusions Olthof ; Olthof2 ; havenith . Second–order Møller–Plesset perturbation calculations (MP2) with medium–sized basis sets (for geometry optimization followed by MP4 single–point energy calculations)LeePark  resulted however in a significantly smaller H–N$`\mathrm{}`$N hydrogen bond angle of 12 together with a reasonable low barrier of 7.6 cm-1. Most recently, another (valence–only) MP2 study Lindh , this time near the one–particle basis set limit, was reported. At the MP2 limit, the hydrogen bond angle of 23 was closer to the recommended value Olthof ; Olthof2 , but the barrier of 24 cm-1 seems too high and the dissociation energy $`D_\mathrm{e}`$ (into separated and optimized monomer fragments) was estimated to be $`13.5\pm \mathrm{0..3}`$ kJ/mol. This clearly hints to the importance of higher–order correlation effects and core–valence corrections, in particular since they account for 0.4 kJ/mol in the case of the water dimer Klopper2  (which already exceeds the error bar given in Ref. Lindh ). DFT calculations of the ammonia dimer are both rather scarce Zhu ; Suhai ; Jursic ; DMT  and somewhat inconclusive as to what level of accuracy DFT methods describe this intermolecular interaction. In the following we will clearly demonstrate that none of the 14 widely used GGA, meta–GGA, and hybrid functionals checked by us lead to a satisfactory description of both the structure and relative energies of the ammonia dimer (despite large basis sets being used and basis set superposition error being corrected for). Thus, even at the level of the ammonia dimer there is still room for improvement based on purely theoretical grounds !

Thus far we considered the behavior of ammonia in the limit of forming a dimer, e.g. in the very dilute gas phase. The question arises, if or how its peculiar behavior concerning hydrogen–bonding in the dimer is present in its condensed phases? Indeed, based on the apparent stability and rigidity of the cyclic dimer in the gas phase MW  it was concluded in quite general terms some time ago that NH3 might well be best described as a powerful hydrogen–bond acceptor with little propensity to donate hydrogen bonds nelson87 ; note that now it is however recognized that the ammonia dimer is non–rigid and non–cyclic havenith . The x–ray structure of crystalline ammonia at low temperatures CS1 ; CS2  features a staggered, directed hydrogen bonded network, with the monomers being arranged like the dimer in Fig. 1d. Furthermore, based on neutron powder diffraction experiments, it was proposed that solid ammonia has an “unusual shared–hydrogen bond geometry” nelmes96  and, in particular, the high–pressure phase ND3IV contains bifurcated hydrogen bonds with bending angles as large as about 30. A most recent theoretical study fortes03 , however, did not support the presence of bifurcated hydrogen bonds in phase IV of solid ammonia. In this latter study, it was concluded through calculations of electron density between intermolecular nitrogen and hydrogen atoms that the apparent bifurcated hydrogen bond geometry in solid ammonia IV is actually a single hydrogen bond perturbed by a neighboring interacting atom. The structure of liquid ammonia has been studied by x–ray AHN  and neutron CB  diffraction for a long time. Only the total structure factor and the total radial distribution function could be extracted in these early experiments and the presence of hydrogen bonds in liquid ammonia was inferred from a comparison of the experimental structure factor with the known structure of the solid. More recently, a set of neutron diffraction experiments was performed RNR  on liquid NH3, ND3, and an isomolar NH3/ND3 mixture. The powerful isotopic substitution technique, in conjunction with sophisticated data analysis, allowed extraction of all partial radial distribution functions at 213 and 273 K. Contrary to earlier conceptions (as presented in many standard textbooks on that subject) the spatial arrangement of nitrogen atoms showed that no extended hydrogen bonded network exists in liquid ammonia. Nevertheless, some degree of hydrogen bonding was inferred from the temperature dependence of the N-H and H-H radial distribution functions. However, the hydrogen bond interaction in liquid ammonia proved to be much weaker than that in water and no clear hydrogen bond peak was observed in either N-H or H-H correlations, unlike the case of water.

Due to the importance of liquid ammonia and in particular its metallic solutions, also many theoretical attempts have been made in past two decades to understand its liquid phase. Most of these studies are based on classical Monte Carlo (MC) or molecular dynamics (MD) simulations using empirical site–site WLJ ; KS1 ; IK ; heinzinger88 ; GXG ; KKH ; beu  or ab initio pair potentials hannongbua00  which all lead to essentially linear hydrogen bonds in the limit of treating the ammonia dimer. A recent pioneering study DMT  investigated the structure of liquid ammonia through Car–Parrinello ab initio molecular dynamics CP ; MH , i.e. without resorting to any pre–determined (pairwise) potentials but using the full many–body interactions as obtained by “on the fly” DFT (BLYP) instead. It was found that the results of the N-N, N-H and H-H radial distribution functions are in better agreement with experimental results than previous ones relying on model potentials. Some important differences, however, still remained. For example, the results of this study were found to overestimate the hydrogen bonded structure as revealed by the N-H correlations and predict a narrower first peak of the N-N distribution when compared with the corresponding experimental results. Also, a later such study LT  of liquid ammonia containing an ammonium or an amide solute reported somewhat different results for the N-N and N-H correlations between the bulk ammonia molecules compared to those of Ref. DMT . Besides, the nature of the hydrogen bonds in the liquid phase is not fully understood yet, in particular in view of the experimental situation. As discussed before, unlike water, the ammonia dimer has a bent hydrogen bond in its equilibrium geometry. If this bent character of the dimer hydrogen bond has the same significant effects on the liquid as in solid phases nelmes96  is still an unresolved issue, in particular since all pair potentials by construction yield a quasi–linear hydrogen bond as do existing density functionals.

The outline of this paper is as follows. In Section II, we describe the details of the computational methodologies that we employ for the calculations of the dimer and the liquid phase. In Section III, we present results for important dimer structures and their energies based on highly accurate ab initio methods and compare these data to those obtained from a wide variety of currently available density functional families. Section IV focuses in detail on the development of the novel HCTH/407+ functional, on its performance in various systems, and on the improvement achieved on the dimer level. In Section V, we compare the results of ab initio MD simulations of liquid ammonia using HCTH/407+ to those obtained from the widely used BLYP functional for reference purposes, and to experimental data wherever available. We summarise and conclude in Section VI.

## II Methods and Technical Details

### II.1 The HCTH Density Functional Family

In this contribution, we develop a new generalised gradient approximation (GGA) functional – which will be termed HCTH/407+ – that only involves the density $`\rho `$(r) and its gradient $`\rho `$(r). The latter has been shown to be a necessary ingredient for the description of hydrogen bonds within DFT, see for instance Refs. kari1 ; sim92 ; cpwatern ; kari2 ; TBH ; tHCTH ; landman . To the contrary of with the so–called hybrid functionals that involve “Hartree–Fock exchange” (such as e.g. B3LYPB3LYP ), this advantage is not impaired by a prohibitive computational overhead relative to the local (spin) density approximation if used in the framework of plane wave pseudopotential calculations, see e.g. the arguments presented in Sect. 3.3 of Ref. MH . Thus, we opted for reparametrising the Hamprecht-Cohen-Tozer-Handy (HCTH) form, which has been described in detail elsewhere HCTH93  and was originally used by Becke for a hybrid functional B97 . In general, HCTH is a post-local spin density approximation (post- LSDA) functional, meaning that it factorises the LSDA functional forms ($`F_{LSDA}`$)

$`E_{xc}={\displaystyle \underset{\gamma =x,c_{\sigma \sigma },c_{\alpha \beta }}{}}E_\gamma ={\displaystyle \underset{\gamma }{}}{\displaystyle \underset{q=0}{\overset{m}{}}}c_{q,\gamma }{\displaystyle F_{LSDA,\gamma }(\rho _\alpha ,\rho _\beta )u_\gamma ^q(x_\alpha ^2,x_\beta ^2)𝑑𝐫},`$ (1)
here, $`u_\gamma ^q`$ denotes the perturbation from the uniform electron gas if $`c_{0,\gamma }`$ is unity. The sum over the integers $`q`$ going from zero to $`m`$ implies the power series defined through $`u_\gamma `$, with $`\theta _{\gamma \sigma }`$ and $`\eta _{\gamma \sigma }`$ as fixed coefficients:

$`u_\gamma ^q(x_\alpha ^2,x_\beta ^2)=\left({\displaystyle \frac{\theta _{\gamma \sigma }x_\sigma ^2}{1+\eta _{\gamma \sigma }x_\sigma ^2}}\right)^q`$ (2)
The variable $`x_\sigma `$ is closely related to the reduced density gradient

$`x_\sigma ^2={\displaystyle \frac{(\rho _\sigma )^2}{\rho _\sigma ^{8/3}}}.`$ (3)
This functional form was previously employed for reparametrisations HCTH120 ; HCTH407  resulting in a highly accurate GGA functional, HCTH/407HCTH407 . Although the above form obeys none of the exact conditions for slowly varying densities Langreth , and violates part of the scaling relations Levy , it does obey the most crucial scaling relations and the Lieb-Oxford bound L-O ; L-O2  within the physically important region Shopping .

When employing the form in Eqn. 1 up to fourth order in $`m`$, we obtain 15 linear coefficients (because of having exchange, spin-like and spin-unlike correlation), which are easily parametrised by minimising an error function $`\mathrm{\Omega }`$, which is constructed out of the sum of three components

$`\mathrm{\Omega }`$ $`=`$ $`{\displaystyle \underset{m}{\overset{n_E}{}}}w_m(E_m^{exact}E_m^{KS})^2+{\displaystyle \underset{l,X}{\overset{n_G}{}}}w_{l,G}\left({\displaystyle \frac{E_l^{KS}}{X}}\right)^2`$ (4)
$`+{\displaystyle \underset{j,\sigma }{\overset{n_v}{}}}w_{j,v}{\displaystyle (v_{j,\sigma }^{ZMP}+k_{j,\sigma }v_{j,\sigma }^{KS})^2\rho _{j,\sigma }^{2/3}𝑑𝐫}.`$
The three sums represent the mean-square deviations from our reference data of energies, gradients and exchange-correlation potentials for each molecule, respectively, of the result of a Kohn–Sham density functional calculation (denoted by the superscript “K-S”). In the case of the total energies reference data are high–level quantum chemistry results that are denoted by “exact”. In the second sum the exact gradients (at equilibrium geometry) are zero by definition and thus do not appear in the formula. In the final term, we fit to the exchange-correlation potentials as determined by the Zhao-Morrison Parr(ZMP)-method ZMP , which are shifted by a constant $`k`$ because of the effects of the quantum-mechanical integer discontinuity. The ZMP method has been proven to be an important aspect of the fit potential . All these contributions need to be balanced using weights $`w`$, which have been determined and reported previously HCTH407 . The weights $`w`$ consist of a product of multiple partial weights making contributions for each molecule in order to ensure a balanced functional.

We also have to consider the molecular set for which the new functional was determined. In general, the new HCTH/407+ functional was determined by fitting it to the “407 set” of molecules as used for HCTH/407 www . This set is similar to the G3 set of molecules G3 , with added inorganic molecules, but a smaller proportion of large organic species which we considered to be over–represented. This original training set was supplemented with data from the ammonia dimer, including non-equilibrium structures, as described in detail in Sect. IV.

### II.2 Quantum Chemical Calculations

For the accurate determination of some points of the potential energy surface (PES) of the ammonia dimer, we used the W2 method W2 . This is basically an extrapolation towards the full CCSD(T) basis set limit including relativistic (but not Born-Oppenheimer) corrections. W2 energies are calculated at the CCSD(T)/APVQZ geometry, this notation meaning an aug-cc-pVQZ basis set for nitrogen and a cc-pVQZ basis set for hydrogen PVXZfirst . The $`MOLPRO`$ package MOLPRO  was used for these calculations. This method is one of the most accurate standard ab initio methods currently available, with an average error of less than 0.5 kcal/mol for the G2-1 setG2-1  of molecules W2eval . In particular, it is more reliable and accurate than the G1, G2 and G3 methods (although computationally much more expensive) W2eval ; Mark .

Two types of higher-order contributions were considered. Firstly, higher-order $`T_3`$ effects were assessed by means of full CCSDTCCSDT  calculations using the $`AcesII`$ program system Aces . Secondly, the importance of connected quadruple excitations was probed by means of CCSD(TQ) CCSDtq  and BD(TQ) BD ; BDt  calculations using $`AcesII`$ and $`Gaussian98`$ g98 , respectively.

For the DFT calculations, we employed the $`Cadpac`$ suite of programs Cadpac , and assessed a variety of GGA and meta–GGA functionals (namely BLYPBLYP , PBEPBE , HCTH/120HCTH120 , HCTH/407HCTH407 , $`\tau `$-HCTHtHCTH , and BP86P86 ; we also tested but do not report results from PW91PW91  and HCTH/93HCTH93 ) and hybrids (B3LYPB3LYP , B97-1HCTH93 , and $`\tau `$-HCTH hybridtHCTH ; not reported are the results obtained from B97-2B97-2  and PBE0PBE0 ). In case of the density functional calculations we used an APVTZ basis set, which is reasonably close to the DFT basis set limit.

For the DFT values, which were obtained from using finite basis sets without extrapolation as done in the W1/W2 methods, we employed the counterpoise (CP) correction BSSE  to account for the basis set superposition error (BSSE). In the fitting procedure, the TZ2P basis setTZ2P  was utilised in addition to the APVTZ basis set.

### II.3 Ab Initio Molecular Dynamics Simulations

The ab initio molecular dynamics simulations were performed by means of the Car-Parrinello method CP ; MH  and the CPMD code MH ; cpmd . A simple cubic box of 32 ammonia molecules with a box length of 11.229 Å was periodically replicated in three dimensions and the electronic structure of the extended system was represented by the Kohn-Sham formulation KS  of DFT within a plane wave basis. The core electrons were treated via the atomic pseudopotentials of Goedecker et al. GTH  and the plane wave expansion of the KS orbitals was truncated at a kinetic energy of 70 Ry. A fictitious mass of $`\mu `$=875 a.u. was assigned to the electronic degrees of freedom and the coupled equations of motion describing the system dynamics was integrated by using a time step of 5 a.u. As usually done, see e.g. Refs. cpwatern ; tuck ; Sprik , the hydrogen atoms were given the mass of deuterium in order to allow for a larger molecular dynamics timestep which also reduces the influence of (the neglected) quantum effects on the dynamical properties.

The ab initio MD simulations have been performed using the HCTH/407+ and BLYP BLYP  functionals. We used the BLYP functional in addition to the new functional because the former has been shown to provide a good description of hydrogen bonded liquids such as water cpwatern ; Sprik ; SP , methanol tsukada  and also ammonia DMT ; LT . Thus, it is worthwhile to compare the results of the two functionals for various structural and dynamical properties of liquid ammonia. The initial configuration of ammonia molecules was generated by carrying out a classical molecular dynamics simulation using the empirical multi-site interaction potential of Ref. GXG . Then, for simulation with the HCTH/407+ (BLYP) functional, we equilibrated the system for 8.75 (10.1) ps at 273 K in NVT ensemble using Nose-Hoover chain method and, thereafter, we continued the run in NVE ensemble for another 9.30 (9.34) ps for calculation of the various structural and dynamical quantities; the average temperatures of these microcanonical runs were about 275 (252) K. We note that the size of the simulation box used in the present simulations corresponds to the experimental density of liquid ammonia at 273 K which is $`2.26\times 10^2`$ molecules/ Å3 RNR .

## III The Ammonia Dimer: Structure and Energetics

In order to briefly validate the W2 method for some representative hydrogen bonded systems, (HF)2 and (H2O)2, we compared the W2 results to MP2 basis set limit corrections applied to coupled cluster calculationsKlopper2 ; Klopper1 , whereas the reference data was additionally empirically refined by scaling the calculations to certain quantum energy levels. Moreover, to get an estimate of the accuracy of the basis set extrapolation, we compared with the more cost-effective alternative of W2, which is W1 W2 . When calculated at the same geometry as W2, it usually yields a fair estimate as to how well the basis set limit is obtained in the extrapolation. In Table 1, the dissociation energies for the HF and H2O dimers are compared to best values from the literature, and displayed along with the results for the (NH3)(H2O) complex. All the W1 and W2 results are within the stated error margins of the reference data. While the dissociation energy of the hydrogen fluoride dimer is exactly the same, that of water dimer is slightly lower than predicted in reference Klopper2 . In the latter case, relativistic corrections were not included in their best estimate. This will reduce the energy somewhat; however we would expect a geometry relaxation to lower the energy. The authors of Ref.Klopper2  used a CCSD(T)/aug-cc-pVTZ geometry and estimated the geometry relaxation, whereas the W2 calculation is carried out at an CCSD(T)/A’PVQZ geometry. Nevertheless, the results show only 0.03 kJ/mol difference between the W1 and W2 method for the HF dimer and 0.09 kJ/mol for the water dimer, indicating the high accuracy of these methods for such interactions. Finally, the (H2O)(NH3) dissociation energy is predicted to be 26.8 kJ/mol by W2, compared to 27.4 kJ/mol determined by experiment Wurzberger ; unfortunately no high-level ab initio data were found in the literature for this mixed dimer.

Let us now focus on the ammonia dimer. Five different structures on the dimer PES can be considered important. The global minimum found on the CCSD(T)/APVQZ surface, as shown in Figure 1c, has an HNN angle of 19.9 with eclipsed hydrogens and Cs symmetry. Its staggered counterpart, displayed in Figure 1d, has an HNN angle of 13.1. The two completely linear structures with the HNN angle being 0 (Figures 1a and 1b) are salient points and the C2h structure (Figure 1e) is a transition state. In Table 2, W1 and W2 results for the five structures are summarised, with column two and four (with BSSE) corresponding to the regular W2 method without the cp-correction. Because of the low barrier between the Cs minimum and the C2h structure, we additionally calculated a counterpoise-corrected W2 barrier. As to the accuracy of the electron correlation methods underlying the W2 method, we can make the following observations:

* quasiperturbative connected triple excitations, i.e., the CCSD(T) – CCSD difference, only contribute 0.5 cm-1 to the C2h barrier with CP correction and 0.7 cm-1 without;
* this has also been verified by CCSD(TQ) \[BD(TQ)\] calculations, where the contribution of the perturbative connected quadruples lowers the barrier by 0.2 cm-1 \[0.3 cm-1\] with the APVDZ basis set;
* higher-order connected triples (i.e., CCSDT-CCSD(T)) additionally lower the barrier by 0.2 cm -1 when using APVDZ. This leads to the conclusion that not only has basis set convergence been achieved, but also that higher-order excitations beyond CCSD(T) will contribute less than 1 cm-1;
* further optimisation of the Cs geometry at the MP2/A’PV5Z level lowers the binding energy by only 0..3 cm-1 (0.004 kJ/mol) relative to MP2/APV5Z//MP2/APVQZ;
* we do not expect relaxation of the geometry at the core-valence level to make a significant contribution, since even for the monomer the core-valence correction to the geometry is very small core .

Still, the difference between the BSSE-corrected and uncorrected results is rather large. Based on these results, we suggest that the Cs to C2h energy difference is 3.5 $`\pm `$ 3 cm-1, and we would expect the barriers between the linear and minimum geometries to be equally well described by W2. For the slightly improved CCSD(T) geometry, by correcting CCSD(T)/APVQZ - MP2/APVQZ + MP2/APV5Z, we predict the HNN angle to be 20.7 degrees. As for the C2h barrier, the interaction energy is raised by 0.002 kJ/mol when going from CCSD(T) to CCSDT, showing that the quasiperturbative triples treatment (and the quadruples estimation as well) is quite accurate. Using an APVDZ basis set, the quasi-perturbative quadruples treatment lowers the binding energy by 0.092 kJ/mol for the coupled cluster and 0.087 kJ/mol when using the Brueckner-Doubles method. With the geometry relaxation contributing less than 0.005 kJ/mol at MP2 level and the estimate of the higher order excitations, we predict the dissociation energy to be 13.1 $`\pm `$ 0.2 kJ/mol, since the main error is likely to come from the quadruple excitations.

The results obtained using a selection of currently published functionals are shown in Table 3. For comparison, the CCSD(T), CCSD, MP2 and HF values are given. It should be noted, however, that BSSE for the coupled-cluster values is more than twice as high as for DFT (0.9 compared to 0.4 kJ/mol), and therefore we would expect the DFT calculations to be closer to their respective basis set limits than the ab initio methods. Nevertheless, their BSSE is still larger than that at the HF level. As expected, the CCSD(T) and W2 numbers are reasonably close to each other. Most density functional methods reproduce the dissociation energy as listed in the first column, reasonably well. In general, the GGA and meta-GGA functionals tested give a good description of this interaction underestimate and PW91 overestimating the dissociation energy. The tested hybrid functionals are even more accurate, yet all of them underestimating the dissociation energy.

However, if we consider the energy difference between the staggered and eclipsed conformation of the NH3-dimer (see the data in the second column of Table 3), it becomes obvious that available functionals are incapable of reproducing the effect of a bent hydrogen bond. It is found that density functionals typically underestimate the difference between the local and global minima by at least a factor of two. We note here that functionals other than those listed in our tables (such as e.g. PW91, B97-2, PBE0) yield similar results which are therefore not displayed. On the other hand, all other wavefunction-based ab initio methods considered are able to reproduce this effect. This becomes even more apparent when considering the C2h transition state. Here, the coupled-cluster differences are close to 3 cm-1, whereas HF renders 24 cm-1 and the DFT values are between 60 and 95 cm-1 with the exception of HCTH/407 (27 cm-1) and the $`\tau `$-HCTH functionals (100 and 140 cm-1). In turn, the other barrier is underestimated by 30-60% and both linear structures are stabilised relative to the global minimum. The same trends are seen in Table 4, where the HNN angles of the two minimum structures are shown. Here, available density functionals underestimate the bending of the hydrogen bond by almost a factor of two for the global minimum and by about 30% for the local (staggered) minimum, although the hydrogen bond distance itself is reasonably well reproduced. Both $`\tau `$-HCTH functionals, while they perform well with their large ’training’ fit set, are the worst performers here. On the other hand, HCTH/407 is the only functional which is somewhat capable of reproducing the effect, outperforming Hartree-Fock.

In summary, although standard density functionals do yield a correct dissociation energy and a slightly nonlinear hydrogen bond, they completely overestimate the barrier towards the symmetrical C2h structure, thereby preferring a much more linear hydrogen-bonding configuration than our accurate reference methods. Although the energy differences might be viewed as intrinsically too small to be accurately rendered by DFT, the fact that even HF gives better results is somewhat disturbing. Nevertheless, the observation that Hartree-Fock and HCTH/407 are somewhat capable of properly predicting the structure and energetics of the ammonia dimer system suggests that the problem at hand is not beyond the grasp of density functional methods as such.

## IV The HCTH/407+ Functional: Construction and Performance

Considering that standard density functionals were shown to be unable to describe the energy difference between linear and bent hydrogen bonds and do not capture the large bending angle of the hydrogen bond, the most obvious remedy would be to change the functional form itself. However, if all tested functionals give such a bad description, this is unlikely to solve the problem. The culprit in this case, therefore, has to be the set of molecules which the HCTH functionals are fit to, insofar as only linearly hydrogen-bonded systems have been included. Nevertheless, including the minimum structure and energetics of the ammonia dimer determined in the previous section into the set probably only partly solves this problem. The main problem, as established before, is the energy difference between the linear and bent hydrogen bonds. In addition, the TZ2P basis set commonly used in the fit set might not be appropriate. Even when employing the APVTZ basis set mentioned above, the basis set superposition error for the dimerisation energy is 29 cm-1, contributing 8 cm-1 to the energy difference between the minimum and symmetric structure when using the B97-1 functional. In comparison, using the TZ2P basis set, the BSSE is 92 cm-1, contributing 33 cm-1 to the energy difference. Clearly, the TZ2P basis set is inappropriate for this problem. The above numbers indicate that a counterpoise procedure has to be employed within the fit in order to get a good estimate. Note that even when using an extrapolation to the full basis set limit (in the case of W2), the difference between the counterpoise corrected and uncorrected values is already 3 cm-1.

All of the above issues were considered in the fit.. Thus, we fit to all points shown in Figure 1 employing an APVTZ basis set with counterpoise correction in the fit. This is the first time we are fitting to non-equilibrium structures, and we only fit to their energy differences and exchange-correlation potentials of these. Due to the energy difference of only a couple of wavenumbers, an additional weight has to be employed for any change in the fit to be noticeable. Thus, equation (3) now reads

$`\mathrm{\Omega }`$ $`=`$ $`{\displaystyle \underset{m}{\overset{n_E}{}}}w_m(E_m^{exact}E_m^{KS}E_m^{BSSE})^2+{\displaystyle \underset{l,X}{\overset{n_G}{}}}w_{l,G}\left({\displaystyle \frac{E_l^{KS}}{X}}\right)^2`$ (5)
$`+{\displaystyle \underset{j,\sigma }{\overset{n_v}{}}}w_{j,v}{\displaystyle (v_{j,\sigma }^{ZMP}+k_{j,\sigma }v_{j,\sigma }^{KS})^2\rho _{j,\sigma }^{2/3}𝑑𝐫}`$
with

$`E_m^{BSSE}=E_m^{ghost,donor}+E_m^{ghost,acceptor}E_m^{monomer,donor}E_m^{monomer,acceptor}.`$ (6)
Since we do not want to fit to a given value of basis set superposition error, it is included in $`E_m^{KS}`$. The weights $`w_{m,l,j}`$ now also include the additional weights for the potential energy surface with the other weights defined as before HCTH407 . Thus, $`w_m`$ is composed of

$`w_{m,PES}=750\times \mathrm{weight}_{confidence}\times \mathrm{weight}_{atom}\times \mathrm{weight}_{dimer}\times \mathrm{weight}_{PES}.`$ (7)
When determining the weights for the new members of the fit set, we need to consider both their APVTZ optimised energies as well as their APVTZ single-point energies at the CCSD(T)/APVQZ optimised geometries. The latter is necessary since we can only fit to single points. As a result of the small energy differences, the frozen-geometry single-point energies and optimised energies differ by a significant amount. The development procedure is similar to the one employed for the HCTH/407 functional, and the set includes numerous atomisation energies, proton affinities, ionisation potentials, electron affinities, and dissociation energies of hydrogen-bonded dimers and transition metal complexes HCTH407 . First, we fit to a subset of the 407 molecules in order to get close to the global minimum, then we fit to the full set. The latter subset (small set) is the G2-1 set plus two hydrogen-bonded systems (H2O)2 and (HF)2, and the additional structures on the potential energy surface, making 153 systems in all (the 147 of the HCTH/147 functional plus the five PES points, in addition to the NH3 dissociation energy at the TZ2P basis). For each point on the PES, the counterpoise correction was employed. The results using the newly determined functionals as a function of the weights have been determined in both the first fit (to 153 systems) and the second fit (to 413 systems). Thus, we obtained new functionals with weights ranging from 0/0 (which corresponds to the HCTH/407 functional) to 200/200. The first weight is the one used in the subset, while the second corresponds to the large fit. The importance of the weights in the first fit can be seen from the differences by using weight pairs of 30/80 and 80/80 \- they differ by about a factor of two in the energy difference between the staggered and eclipsed geometry. Hence, the minimum obtained by the second fit will of necessity be close to the global minimum obtained by the first one. When increasing these weights, we approach the reference values as expected, with the 200/200 functional rendering the closest values to the reference method. While not explicitly included in the fit, the BSSE also increased from 0.444 kJ/mol to 0.866 kJ/mol. This resembles the tendency of the functional to have a larger intermolecular distance than the reference methods when increasing the weights, hence entering a regime where the BSSE increases more rapidly than for the optimised structures. In comparison, the largest BSSE observed by the functionals tested in Table 3 is 0.48 kJ/mol. However, for all the calculations in the fit we have performed only DFT single-point energy calculations at the CCSD(T) optimised geometries.

Another, more appropriate, point of comparison consists of the optimised structures of all the newly obtained functionals, as done in Table 3 for the contemporary functionals. Since the energy differences in general are very small, these results differ significantly from the ones obtained by the single-point energies calculated at the coupled-cluster geometries. In addition, the C2h structure becomes the global minimum for the functionals with a weight of 80 or larger in the first fit. Based on these results, the 30/80 functional appears to be the most reliable for this interaction, and hence will be denoted as the new HCTH/407+ functional. This new functional (compare to Table 3) now renders a dissociation energy of 13.18 kJ/mol for the ammonia dimer, in much closer agreement to the reference values than any other method employed. The difference between the Cs minimum and the C2h structure is 4.0 cm-1 for HCTH/407+, compared to our best estimate of 3.5 cm-1 in Table 2. The staggered Cs structure lies 18 cm-1 above its eclipsed counterpart, which is about five wavenumbers higher than our best estimate but still within its error bars. The difference to the linear eclipsed and linear staggered structures are slightly underestimated at 51.5 and 53.5 cm-1 respectively, but this is still a vast improvement over the other functionals. In Figure 2, we compare the energies of the new HCTH/407+ functional to W2 and one of the most commonly used GGA functionals (BLYP). The BSSE for the relaxed Cs structures does not increase as rapidly as for their CCSD(T) optimised counterparts. It increases from 0.43 (HCTH/407) to 0.45 (HCTH/407+) to 0.61 kJ/mol (weight of 200), and does not yield as unreasonably large values as the single-point calculations at the CCSD(T) geometries. The hydrogen bond lengths of the Cs and C2h structures increase with larger weights, emphasising the importance of this analysis. Generally, all HCTH functionals are found to slightly overestimate this hydrogen bond length. The HNN angle of the HCTH/407+ functional is 17  for the minimum geometry. This is within the accuracy that can be expected from contemporary density functionals ($`\pm `$ 2), showing a significant improvement over all other density functionals in Table 4. The results for all functionals obtained with different weights (from which we determined the HCTH/407+ functional) are given in the supplementary material. supplementary

The parameters which define the HCTH/407+ functional, compared to HCTH/407, are given in Table 5. As might be expected, they differ only slightly from those for the standard HCTH/407 functional; only the correlation and higher-order coefficients have significantly been affected by the change.. Nevertheless, the new functional should now be able to describe non-directed hydrogen bonds better than HCTH/407. On the level of the ammonia dimer, the improvement with respect to other density functionals becomes clearly visible in Fig. 2; see Table 3 for the corresponding numbers. It is seen that HCTH/407+ describes the relative stability of all isomers much better with respect to the W2 reference data than BLYP, which was selected as a prominent representative of the GGA family. Most importantly, the dramatic failure of BLYP to capture the stability of the cyclic C2h isomer relative to all other isomers is cured. According to HCTH/407+, the two $`C_s`$ structures as well as the $`C_{2h}`$ structure are essentially degenerate, which is in accord with the reference data.. Apart from the energetics, also the structure of the ground state of the ammonia dimer, i.e. the eclipsed $`C_s`$ structure, is dramatically improved, see Table 4 for details. In particular the hydrogen bond angle HNN is now much larger and thus closer to the reference value than any other density functional method; note that also the functionals PW91, B97-2, and PBE0 were considered. The same is true for the HNN angle in the staggered variant of the $`C_s`$ symmetric structure. However, it is also clear from this table that the N$`\mathrm{}`$H distance is clearly overestimated: HCTH/407+ yields a distance of about 2.5 Å instead of around 2.3 Å as required. This seems to be typical for the HCTH family as both HCTH/120 and HCTH/407 yield similarly large values exceeding 2.4 Å (and HCTH/93 leads to 2.614 Å). In addition, the same trend of producing too long hydrogen bonds is found for the $`C_{2h}`$ structure. This deficiency is corrected by both $`\tau `$-HCTH and $`\tau `$-HCTH hybrid functionals. Here it should be noted that both BLYP and PBE describe the hydrogen bond length quite well, however at the expense of making it much too linear and grossly disfavoring the cyclic structure energetically.

Furthermore, this significant improvement in describing non-directed hydrogen bonds, however, comes at the expense of a slightly increased error for the other molecules in the fit set. Table 6 shows the errors of the HCTH/407 and HCTH/407+ functionals to the fit set, summarising these as RMS energy and geometry errors, together with the errors for the energies, geometry shifts and frequency shifts of nine hydrogen-bonded systems. The values for B3LYP, BLYP and BP86 are included for comparison. The hydrogen-bonded systems are the (HF)2, (HCl)2, $`(\mathrm{H}_2\mathrm{O})_2`$, (CO)(HF), (OC)(HF), (FH)(NH3), (ClH)(NH3), $`(\mathrm{H}_2\mathrm{O})(\mathrm{NH}_3)`$ and $`(\mathrm{H}_3\mathrm{O}^+)(\mathrm{H}_2\mathrm{O})`$ complexes, but not the NH3 dimer itself. The performance of other functionals over all the sets has been discussed elsewhere TBH ; HCTH407 ; tHCTH . Overall, errors for the HCTH/407+ functional are quite similar to HCTH/407, which in turn means that they are still excellent compared to first-generation hybrid functionals like B3LYP and especially GGA functionals like BLYP or BP86. Although we expect the geometries to be still slightly worse than B3LYP (although the gradient error is comparable), the energetic properties are better described, as are the hydrogen bonds. Whereas the errors of the HCTH/407+ functional for the fit set have barely changed compared to HCTH/407, the accuracy of the dissociation energies of the nine directed hydrogen bonds has actually decreased by about 20%. Nonetheless, the accuracy of the frequency shifts of the H-X bond stretches has increased, indicating that we seem to have improved the description of the potential energy surface of these complexes.

## V Liquid Ammonia: Structure, Dynamics, and Hydrogen Bonding

### V.1 Radial Distribution Functions

We have calculated the nitrogen-nitrogen, hydrogen-hydrogen and nitrogen-hydrogen radial distribution functions (RDF) from the atomic trajectories generated by ab initio molecular dynamics simulations; note again that we used the mass of the deuteron instead of that of the proton but we still use the symbol H for convenience. The results for both HCTH/407+ and BLYP functionals are shown in Fig. 3. In this Figure, we have also included the experimental results RNR ; Soper  of the RDFs of liquid ammonia at 273 K; note that the underlying experimental technique is based on the assumption that the structure of all isotopically substituted systems is identical. For the nitrogen-nitrogen RDF, it is seen that at the short distances the BLYP $`g_{NN}(r)`$ is in better agreement with experimental results. The HCTH/407+ peak is located at a distance of about 0.15 Å larger than the experimental peak position whereas the BLYP peak appears at a somewhat shorter distance. We note in this context that the HCTH/407+ functional gives an overestimated equilibrium N-N distance for the isolated dimer compared to the reference CCSD(T) result as can be inferred from Table 4. The larger N-N distance seems to be translated to the liquid phase configurations and we observe a somewhat larger most probable N-N distance for the HCTH/407+ functional than what is found in the experiments. Still, the overall shape of the first peak is represented reasonably well by both the functionals. When the RDFs are integrated up to their first minimum to obtain the coordination number in the first solvation shell, we obtained the values of 13.2 and 12.1 for HCTH/407+ and BLYP functionals compared to an experimental value of 12.75. Clearly, the agreement with experimental results is found to be rather good.

Considering the results of the H-H and N-H RDFs, the positions of the calculated intramolecular H-H and N-H peaks for both the functionals agree rather well with experimental results as expected. For example, the liquid phase N-H bond length and intramolecular H-H distance for the HCTH/407+ (BLYP) functionals are 1.03 (1.03) Å and 1.63 (1.65) Å which can be compared with the experimental values of 1.04 Å and 1.60 Å, respectively. Also, these intramolecular distances in the liquid phase are very close to their gas phase values. This was also noted earlier for the BLYP functional DMT . Because of the quantum disperson effects which are present in real liquid but not in the current simulated systems, the experimental intramolecular peaks are found to be broader than the theoretical peaks. The somewhat higher value of the experimental H-H correlation at the first minimum at around 2 Å can also be attributed to such quantum dispersion effects.

The intermolecular H-H peaks are found to be well represented by both functionals. The intermolecular part of the N-H RDF shows the presence of both H-bonded and non H-bonded H atoms in the first solvation shell. The shoulder up to about 2.7 Å can be attributed to hydrogen bonded H atoms, while the more pronounced peak at the larger distance at around 3.75 Å corresponds to the hydrogen atoms in the solvation shell which are not hydrogen bonded. Again, both functionals are found to describe well the non-hydrogen bonded first solvation shell peak when compared with the corresponding experimental result. For the hydrogen bonded part, however, the HCTH/407+ result is found to exhibit a less pronounced shoulder than what is found in the experiments and for the BLYP functional. This is probably a result of the slight overestimation of the hydrogen bond length by HCTH/407+ already on the level of the dimer, see Table 4. Considering both hydrogen bonded and non-hydrogen bonded H atoms, a sphere of radius 5.3 (5..2) Å around a central N atom is found to contain 34.1 (34) intermolecular H atoms for the HCTH/407+ (BLYP) functional. The radius of the sphere is set to 5.3 (5.2) Å because the first minimum of the intermolecular N-H RDF is located at this distance. The corresponding radius for the experimental N-H RDF is 5.0 Å, yielding the experimental value of 34.8 for the above coordination number. A more detailed analysis of the distribution of hydrogen bonds in liquid ammonia is given in the following subsection, and the dynamical aspects of hydrogen bonds are considered in Sect. V.4.

### V.2 Distribution of Hydrogen Bonds

The analysis of the hydrogen bond (HB) distribution among ammonia molecules is based on a calculation of the fractions $`P_n`$ of ammonia molecules that engage in $`n`$ hydrogen bonds and the average number of hydrogen bonds that an ammonia molecule donates and accepts noteHB . Following previous work on liquid ammonia, water and other hydrogen bonding liquids DMT ; GXG ; KKH ; GSR ; LC ; XB ; SNS ; AC ; BPB ; FHM , we have adopted a geometric definition of the hydrogen bonds where two ammonia molecules are assumed to be hydrogen bonded if they satisfy the following configurational criteria with respect to N-N and N-H distances

$`R^{(NN)}`$ $`<`$ $`R_c^{(NN)},`$
$`R^{(NH)}`$ $`<`$ $`R_c^{(NH)}`$ (8)
and the values of these distance cut-offs are usually determined from the positions of the first minimum of the corresponding (intermolecular) RDFs. We have used a similar procedure where the N-N RDF gives a value of R$`{}_{c}{}^{}{}_{}{}^{NN}`$ =5.25 (5.1) Å for the N-N cut-off for HCTH/407+ (BLYP) functionals. However, the procedure can not be applied unambiguously to determine the N-H cut-off because no clear minimum that separates the hydrogen bonded and nonbonded H atoms is found in the N-H RDF. In the latter case, the shoulder which is assigned to hydrogen bonded H atoms merges into the broad and more intense peak corresponding to the H atoms of the first solvation shell which are not hydrogen bonded. Since this shoulder – which corresponds to the hydrogen bonded H atoms – extends up to about 2.7 Å in the experimental and BLYP RDFs we have taken this distance of 2.7 Å as the cutoff N-H distance. For HCTH/407+ also, we have used the same cut-off of the N-H distance although the hydrogen bonded shoulder is less pronounced for this functional. Still, variation of $`R_c^{NH}`$ within reasonable bounds does not change qualitatively the result of this analysis.

In Fig. 4, we have shown the distribution of donor and acceptor hydrogen bonds for both HCTH/407+ and BLYP functionals. In this figure, panels (a) and (b) show the fraction of ammonia molecules that accept or donate $`n`$ hydrogen bonds and panel (c) shows the fraction of hydrogen atoms that donate $`n`$ hydrogen bonds. The majority of the ammonia molecules is found to accept one hydrogen bond and to donate one. Also, the fraction of hydrogen atoms donating two or more hydrogen bonds is found to be negligibly small, which means that the bifurcated or multifurcated hydrogen bonds are practically absent in liquid ammonia. The average number of donor hydrogen bonds per ammonia molecule, which is also equal to the average number of acceptor hydrogen bonds per molecule, is found to be 0.93 (1.34) for the HCTH/407+ (BLYP) functional. Given that there are as many as about 13 ammonia molecules in the first solvation shell, we conclude that non–hydrogen–bonded interactions, i.e. packing or steric effects, are crucial in determining the solvation behavior of ammonia in the liquid state. This is in sharp contrast to the case of liquid water where the solvation structure of a water molecule is determined primarily by hydrogen bonding interactions FHS .

The ammonia dimer has a strongly bent hydrogen bond in its equilibrium geometry and it would be interesting to investigate whether the hydrogen bonds remain bent in the liquid phase. We have carried out such an analysis by calculating the distribution function of the cosine of the hydrogen bond angle $`\theta `$, which is defined as the N$`\mathrm{}`$ N–H angle of a hydrogen bonded pair and the results are shown in Fig.. 5a for both the functionals. It is seen that, for both the functionals, the probability distribution $`P(\mathrm{cos}\theta )`$ is peaked at $`\mathrm{cos}\theta =1`$ indicating that the preferred hydrogen bond geometry is linear in the liquid phase, see Ref. cone-corr . The peak is somewhat more sharp and narrower for the BLYP functional, which means a stronger preference for the linear hydrogen bonds for this functional as compared to that for the HCTH/407+ functional. This preferential presence of linear hydrogen bonds in the liquid phase is indeed an interesting result given that the dimer has a bent hydrogen bond in the gas phase. We have also investigated the distribution of N$`\mathrm{}`$ N–H angle for H atoms that belong to nearest neighbours but are not hydrogen bonded, i.e. the N$`\mathrm{}`$H distance is greater than 2.7 Å. The results of these angular distributions are shown in Fig. 5b for 2.7 Å $`<`$ R(NH) $`<`$ 3.7 Å and in Fig. 5c for 3.7 Å $`<`$ R(NH) $`<`$ 4.7 Å. A wide range of values are found for the N$`\mathrm{}`$ N–H angle with very little or no preference for a particular orientation. These findings are indeed consistent with the above analysis that the first solvation shell in liquid ammonia is dominated by simple packing rather than by directional hydrogen bonding.

### V.3 Self Diffusion and Orientational Relaxation

The translational self diffusion coefficient of an ammonia molecule in the liquid state is calculated from the long-time limit of the mean-square displacement (MSD)

$$D=\underset{t\mathrm{}}{lim}\frac{<[r(t)r(0)]^2>}{6t},$$
(9)
where $`r(t)`$ is the position of the center of mass of a molecule at time $`t`$ and the average is carried out over the time origin for autocorrelation and over all the molecules as usual. The results of the diffusion coefficients, which are obtained from a least square linear fit of the simulation data excluding the initial ballistic regime up to 0.5 ps, are included in Table (VII). This Table also contains the experimental values of the diffusion coefficients of deuterated liquid ammonia ND3 at 275 and 252 K that are found from Eq.(2) of Ref.OPS , which was shown to fit the experimental data over a rather wide temperature range very well. The agreement of the calculated diffusion coefficients with the experimental result is found to be reasonably good for both functionals. In particular it is gratifying to see that the experimental value is underestimated because quantum effects on nuclear motion, which are neglected in the present study, tend to increase the diffusion in a liquid.

The self-orientational motion of ammonia molecules is analyzed by calculating the orientational time correlation function md-ana , $`C_l^\alpha (t)`$, defined by

$$C_l^\alpha (t)=\frac{<P_l(e^\alpha (t)e^\alpha (0))>}{<P_l(e^\alpha (0)e^\alpha (0))>},$$
(10)
where $`P_l`$ is the Legendre polynomial of order $`l`$ and $`e^\alpha `$ is the unit vector which points along the $`\alpha `$-axis in the molecular frame. Here, we have studied the time dependence of $`C_l^\alpha (t)`$ for $`l=1,2`$ and for three different $`e^\alpha `$, the unit vectors along the molecular dipole axis, an N-H bond and an intramolecular H-H axis. The geometric dipole vector of an ammonia molecule is calculated by assigning partial charges to N and H atoms as given by the classical model of Ref.GXG . For the orientational relaxation of unit vectors along the N-H and H-H axes, the results are averaged over three such intramolecular axes of each type. After an initial transient non-exponential decay, the relaxation becomes diffusional and $`C_l^\alpha (t)`$ decays exponentially. The orientational correlation time, $`\tau _l^\alpha `$, is defined as the time integral of the orientational correlation function

$$\tau _l^\alpha =_0^{\mathrm{}}𝑑tC_l^\alpha (t).$$
(11)
We note that the orientational relaxation of molecular vectors containing H atoms are usually measured by NMR relaxation experiments, which yield the Fourier transform of a correlation function for $`l=2`$ such that the integrated correlation time $`\tau _2^\alpha `$ is measured.

In Fig. 6, we depict the results of the orientational correlation function $`C_l^\alpha (t)`$ for $`l=1,2`$ and $`\alpha `$=dipole, N-H and H-H. The corresponding results of the orientational correlation times are included in Table (VII). We have also included the experimental result of $`\tau _2`$ of deuterated ammonia $`ND_3`$ as given by the fitted Arrhenius law of Ref. HZZ . Since this experimental relaxation time was obtained by NMR relaxation experiments, it is likely to correspond to the relaxation of N-H/H-H vectors rather than the dipole vector. The BLYP functional is found to predict a relatively slower orientational relaxation. Considering the quite small system size and short simulation trajectories, both functionals are found to describe the single-particle rotational dynamics of liquid ammonia reasonably well.

### V.4 Dynamics of Hydrogen Bonds

The key dynamical quantity in the context of hydrogen bond dynamics is the mean hydrogen bond life time, $`\tau _{HB}`$. The fast librational motion of ammonia molecules makes apparent breaking and reformation of a hydrogen bond a very fast process and, depending on how these fast transient events are taken into account, one can obtain different values for the hydrogen bond (HB) lifetime. The rotational motion of ammonia molecules is also a rather fast process and can also contribute to the short–time dynamics of hydrogen bonds. For a particular chosen definition of a breaking event, one can adopt a direct counting method FHM ; HFM ; RDM  or a time correlation function method LC ; XB ; SNS ; AC ; BPB ; RAP  to calculate the lifetime. In the present work, we have adopted the time correlation function approach which allows us to calculate the hydrogen bond lifetime in both cases: (a) The breaking of hydrogen bond occurs due to fast librational motion even though it may be a transient event and not a ‘true’ breaking event and (b) these fast librational breaking events are ignored and only rotational and translational diffusional breaking on a longer time scale is considered as the ‘true’ breaking of a hydrogen bond. Also, in the present context, it is implied that ‘fast librational breaking’ also includes the effects of fast rotational motion that might contribute to the short time dynamics of the hydrogen bonds.

In order to obtain the hydrogen bond lifetime for the two scenarios mentioned above, we calculate three time correlation functions: one continuous and two intermittent hydrogen bond correlation functions LC ; XB ; RAP . Before we define these time correlations, we first define two hydrogen bond population variables $`h(t)`$ and $`H(t)`$: $`h(t)`$ is unity when a particular tagged pair of ammonia molecules is hydrogen bonded at time $`t`$, according to the adopted definition from Sect. V.2 and zero otherwise. The function $`H(t)`$ is unity if the tagged pair of ammonia molecules remains continuously hydrogen bonded from $`t=0`$ to time $`t`$ given that the bond was formed for the last time at $`t=0`$, and it is zero otherwise. We define the continuous hydrogen bond time correlation function $`S_{HB}(t)`$ as

$$S_{HB}(t)=<h(0)H(t)>/<h>,$$
(12)
where $`<\mathrm{}>`$ denotes an average over all hydrogen bonds that are created at $`t=0`$. Clearly, $`S_{HB}(t)`$ describes the probability that a hydrogen bonded ammonia pair, which was created at $`t=0`$, remains continuously bonded upto time $`t`$. It differs from the continuous hydrogen bond correlation function of Rapaport RAP  in that we apply the condition of bond formation at time $`t=0`$. The mean hydrogen bond lifetime $`\tau _{HB}`$ can then be calculated from the time integral

$$\tau _{HB}=_0^{\mathrm{}}S_{HB}(t)𝑑t.$$
(13)
In Fig. 7a, we have shown the decay of $`S_{HB}(t)`$ for both functionals and the corresponding results of $`\tau _{HB}`$ are included in Table VII. The hydrogen bond lifetime is found to be about 0.1 ps for both functionals. We note that both fast librational and slower diffusional motion can contribute to the decay of $`S_{HB}(t)`$. Since the librational motion occurs on a faster time scale and since the correlation function $`S_{HB}(t)`$ does not allow any reformation event, the dynamics of $`S_{HB}(t)`$ primarily reveal the dynamics of hydrogen bond breaking due to fast librational motion and hence the lifetime that is obtained from Eq.(12) corresponds to the average time over which a hydrogen bond survives before it ‘breaks’ due to librations. Again, this short time breaking may also include contributions from fast rotational motion of ammonia molecules which also occurs at subpicosecond time scale.

A different way to analyze the hydrogen bond dynamics is to calculate the intermittent hydrogen bond correlation function LC ; RAP

$$C_{HB}(t)=<h(0)h(t)>/<h>,$$
(14)
where the average is now over all hydrogen bonds that were present at $`t=0`$. Note that the correlation function $`C_{HB}(t)`$ does not depend on the continuous presence of a hydrogen bond. It describes the probability that a hydrogen bond is intact at time $`t`$, given it was intact at time zero, independent of possible breaking in the interim time. Clearly, bonds which were briefly ‘broken’ by fast librational motions would continue to contribute to the correlation function at later times and this leads to a much slower decay of $`C_{HB}(t)`$ at longer times. The results of $`C_{HB}(t)`$ are shown in Fig. 7b which clearly shows a two-phase relaxation of this correlation function. The initial fast relaxation corresponds to rapid breaking of hydrogen bonds due to librational motion and the slower relaxation after this initial transient decay corresponds to the breaking and reformation of hydrogen bonds due to both rotational and translational diffusion of ammonia molecules.

After a hydrogen bond is broken, the two ammonia molecules can remain in the vicinity of each other for some time before either the bond is reformed or the molecules diffuse away from each other. We define $`N_{HB}(t)`$ as the time dependent probability that a hydrogen bond is broken at time zero but the two molecules remain in the vicinity of each other i.e. as nearest neighbors but not hydrogen bonded at time $`t`$.. Following previous workLC , we can then write a simple rate equation for the ‘reactive flux’ $`dC_{HB}/dt`$ in terms of $`C_{HB}(t)`$ and $`N_{HB}(t)`$

$$\frac{dC_{HB}(t)}{dt}=kC_{HB}(t)k^{}N_{HB}(t),$$
(15)
where $`k`$ and $`k^{}`$ are the forward and backward rate constants for hydrogen bond breaking. The inverse of $`k`$ can be interpreted as the average lifetime of a hydrogen bond.

The derivative of the intermittent hydrogen bond correlation or the ‘reactive flux’ of Eq.(14) is computed from the simulation results of $`C_{HB}(t)`$ that are presented in Fig. 7b. The probability function $`N_{HB}(t)`$ is also calculated from the simulation trajectories through the following correlation function approach XB

$$N_{HB}(t)=<h(0)[1h(t)]h^{}(t)>/<h>,$$
(16)
where $`h^{}(t)`$ is unity if the N–N distance of the pair of ammonia molecules is less than R$`{}_{}{}^{(NN)}{}_{c}{}^{}`$ at time $`t`$ and it is zero otherwise. The results of $`N_{HB}(t)`$ are shown in Fig.. 7c. We used a least-squares fit of Eq.(14) to the simulation results of the reactive flux, $`C_{HB}(t)`$ and $`N_{HB}(t)`$ to produce the forward and backward rate constants. We performed the fitting in the short time region $`0<t<0.15`$ ps to obtain the rate constants and the corresponding average hydrogen bond lifetime for the librational relaxation and we also carried out the fitting on the longer time region 0.25 ps $`<`$ t $`<`$ 3 ps to calculate these quantities for slower rotational and translational diffusional relaxation. The inverses of the corresponding forward rate constants, which correspond to the average hydrogen bond lifetimes and which we denote as $`1/k_{\mathrm{short}}`$ and $`1/k_{\mathrm{long}}`$, are included in Table VII. We note that the value of $`1/k_{\mathrm{short}}`$ is very close to average hydrogen bond lifetime $`\tau _{HB}`$ obtained from the continuous hydrogen bond correlation function $`S_{HB}(t)`$ which is not unexpected because both $`S_{HB}(t)`$ and the short-time part of the ‘reactive flux’ captures the hydrogen bond ‘breaking’ dynamics due to fast librational motion. The hydrogen bond lifetime as given by $`1/k_{long}`$ is significantly longer because it excludes the fast librational relaxation.

## VI Summary and Conclusions

We have investigated the properties of ammonia from the gas phase dimer to the liquid state by means of wavefunction based quantum chemistry techniques, density functional theory, and ab initio molecular dynamics. The ammonia dimer is a particularly interesting hydrogen–bonded system, as it is known to feature a strongly bent hydrogen bond. In order to investigate possible consequences of this nonlinear hydrogen bond, important points on the potential energy surface were first investigated by means of ab initio methods up to W2–extrapolated CCSD(T) theory, see Fig. 1 for the corresponding structures. The structure and energetics of important isomers are obtained up to unprecedented accuracy without resorting to experiment; see Fig. 1 and Table 4 for structures, and Fig. 2 and Table 3 for energies. It is confirmed that the structure of the global minimum (having eclipsed $`C_s`$ symmetry) possesses a substantially nonlinear hydrogen bond with an H–N$`\mathrm{}`$N angle that is predicted to be 20.7; note that the corresponding angle is predicted to be 5.5 in the water dimer. The energy difference between the global minimum and the C2h transition state (which is the famous “cyclic structure”) is about 3.5 cm-1, and the energy penalty to make the hydrogen bond artificially linear amounts to only 70 cm-1, i.e. 0.2 kcal/mol or 0.009 eV. This implies that the potential energy surface is very flat.

However, when investigating the ammonia dimer using a wide variety of available density functionals (including important representatives of the GGA, meta–GGA and hybrid functional families) it was found that none of the functionals checked by us describes the bent hydrogen bond in (NH3)2 satisfactorily. Typically, the hydrogen bond angle is too small by a factor of two and the energy of the cyclic $`C_{2h}`$ structure w.r.t. the eclipsed $`C_s`$ global minimum is overestimated by a factor of 10 to 30. Since ammonia is an important hydrogen bonded liquid, a density functional was developed for studying this subtle system, with particular focus on its applicability for condensed phase simulations in the framework of Car–Parrinello molecular dynamics. To this end, a novel GGA–type density functional, HCTH/407+, was developed, with special emphasis on the nonlinear hydrogen bond and potential energy surface of (NH3)2. This functional yields the bent hydrogen–bonded equilibrium structure of the ammonia dimer, as well as the correct energetics of the low–lying isomers including the energy barrier to a linear hydrogen–bond.

The performance of the new functional in describing the structural and dynamical properties of liquid ammonia was investigated by carrying out Car–Parrinello molecular dynamics simulations of the liquid phase. In addition, simulations were carried out using a GGA–type functional that is widely used for describing associated liquids, BLYP, and the results of the two functionals are compared with experiments wherever available. In particular, we focused on the atom-atom radial distribution functions, on structure and dynamics of hydrogen bonds, and on diffusion as well as orientational relaxation processes. Overall, both functionals are found to describe the structural and dynamical properties of the liquid phase reasonably well. Importantly, it is shown that the propensity to form a strongly bent hydrogen bond – which is characteristic for the equilibrium structure of the gas–phase ammonia dimer – is overwhelmed by steric packing effects that clearly dominate the solvation shell structure in the liquid state. Similarly, the propensity of ammonia molecules to form bifurcated and multifurcated hydrogen bonds in the liquid phase is found to be negligibly small. Thus, even functionals that lead to unreasonably linear hydrogen bonds in the limiting case of the in vacuo ammonia dimer, such as BLYP, yield a good description of liquid ammonia – albeit for the wrong reason!

###### Acknowledgements.

We are grateful to Professor Martina Havenith for pointing out the pecularities of the ammonia dimer and for stimulating discussions. The authors would also like to acknowledge Professor Nicholas Handy for enlightening discussions, his critical comments on an early draft of the manuscript, and his help in the very beginning of the project. A. D. Boese is grateful for financial support by the EPSRC, the Gottlieb Daimler- und Karl Benz-Stiftung, and the Feinberg Graduate School. A. Chandra wishes to thank Alexander von Humboldt –Foundation for a Research Fellowship. The simulations were carried out at Rechnerverbund–NRW and at Bovilab@RUB (Bochum). This research was also supported by the Lise Meitner Center for Computational Chemistry, the Helen and Martin A. Kimmel Center for Molecular Design, DFG, and FCI.  References

1 APVQZ basis set

a at 275 K, b at 252 K

Fig. 1. The five relevant optimized structures of the ammonia dimer; dotted lines are only guides to the eye..

Fig. 2. The potential energy profile obtained from W2 Theory (accurate reference calculations) in comparison to that obtained from HCTH/407+ and BLYP functionals at the five optimized structures, see insets, from Fig. 1. Note that the profiles obtained from other GGA functionals such as PBE are qualitatively similar to that of BLYP.

Fig.3. The atom-atom radial distribution functions. The solid, dashed and the dotted curves show the HCTH/407+, BLYP, and the experimental RNR ; Soper  results, respectively.

Fig.4. The fraction $`P(n)`$ of molecules having a certain number of accepted $`n_{\mathrm{acceptor}}`$ and donated $`n_{\mathrm{donor}}`$ hydrogen bonds, see panels (a) and (b), respectively. Panel (c) shows the fraction of hydrogen atoms that donates a certain number of hydrogen bonds. The squares and the triangles represent the results for the HCTH/407+ and BLYP functionals, respectively, and the dashed lines are drawn to guide the eye.

Fig.5. The distribution of the cosine of the N$`\mathrm{}`$N–H angle for H atoms that belong to the nearest neighbours. Panel (a) shows the distribution for hydrogen bonded H atoms (i.e.. R(NH) $`<`$ 2.7 Å) and panels (b) and (c) show, respectively, the results for non–hydrogen–bonded H atoms that appear in the region 2.7 Å $`<`$ R(NH) $`<`$ 3.7 Å and 3.7 Å $`<`$ R(NH) $`<`$ 4.7 Å. The solid and the dashed curves show the results for the HCTH/407+ and BLYP functionals, respectively.

Fig.6. The time dependence various orientational correlation functions, see text for definitions and further details. The solid and the dashed curves show the results for the HCTH/407+ and BLYP functionals, respectively.

Fig.7. The time dependence of various hydrogen bond correlation functions, see text for definitions and further details. The different curves are labelled as in Fig.6.


format: text

text:
# Dynamical Processes in the Central Kpc and Active Galactic Nuclei

## 1. Introduction: Emerging Global Connection

Departures from axial symmetry are destined to shorten the timescales of secular and dynamical evolution in disk galaxies. Stellar bars, triaxial halos and tidal interactions play the important roles of driving such evolution on larger spatial scales, but their effect diminishes sharply within the central kpc. Is there any comparable non-axisymmetric morphology within the circumnuclear region which can impose gravitational torques, trigger bursts of star formation, and fuel the nonstellar activity of supermassive black holes (SBHs) — ubiquitous, as recent observations confirm? What are the relevant processes which maintain this morphology, and how much they affect the galactic evolution?

Disk galaxies as a “norm” are barred or ovally distorted in the near-infrared (NIR). High-resolution ground-based instruments and the availability of the HST has allowed for the first time a meaningful analysis of central morphology and kinematics. Although our knowledge of the inner regions of disk galaxies is clearly incomplete, certain patterns in their dynamical evolution and their relationship to larger and smaller spatial scales have emerged.

At least dynamically, the inner parts of disk galaxies can be defined by the positions of the inner Lindblad resonances (ILRs), typically at about 1 kpc from the center. These resonances play an important role in filtering density waves, propagating between the bar corotation radius (CR) and the center. They are usually delineated by elevated star formation rates and the concentration of molecular gas in nuclear rings, which can serve as reservoirs for fueling the central activity.

One can distinguish the disk and spheroidal components within the central kpc. Surprisingly, the mass of the central SBH has been claimed to correlate with the bulge properties (e.g., reviews by Kormendy & Gebhardt 2001; Merritt & Ferrarese 2001), despite the common wisdom that black holes are fueled by disk accretion. However, formation of bulges can be tied directly to the properties of triaxially-shaped halos, thus strongly hinting about the ‘global connection’ between smallest and largest spatial scales within the forming and evolving galaxy (Section 3 and El-Zant et al. 2002).

Resolved plane morphology of central kpc in barred galaxies has revealed so-far grand-design (mini) spirals and bars in addition to “traditional” disks and bulges. We refer to large, kpc-scale bars as “primary,” and to the sub-kpc bars as “secondary.” The theoretical rationale behind these definitions is that secondary bars are believed to form as a result of radial gas inflow along the large-scale bars and, therefore, are expected to be confined within their ILRs (Shlosman, Frank & Begelman 1989). Below we review morphologies of the circumnuclear region, analyze dynamical properties of nested bars, their dynamical coupling to larger scales, and the resulting gas flows on scales of $`10`$ pc$`10`$ kpc, as well as discuss the possible origin of SBH-bulge correlation within the context of triaxial halos.

## 2. Nuclear Bars in Nested Bars: Gas, Stars and ‘Cocktails’

In the powerful display of nonlinear dynamics, many galaxies exhibit double stellar bars, which can tumble with different or identical pattern speeds. Although examples have been known for nearly three decades, the dynamical importance of nested bars has been first pointed out much later. The largest sample of (112) disk galaxies analyzed so far for this purpose reveals a substantial fraction of nested bars, probably in excess of $`2025`$%. Even more interestingly, about 1/3 of barred galaxies host a second bar (Laine et al. 2002). The gas contents of nuclear bars vary. In some cases, the cold gas can be dynamically important, as evident from the interferometric 2.6 mm CO emission and NIR lines of H2. Depending on the gas fraction contributing to the gravitational potential, one can distinguish stellar-, gas-dominated (i.e., gaseous) and mixed nuclear bars. The nuclear gaseous bars have no large-scale counterparts.

Nuclear bars extend the action of gravitational torques to smaller spatial scales. Their relevance for the AGN depends of course on the availability of the ‘fuel’ (i.e., gas) which loses angular momentum and falls toward the center. The gas-dominated nuclear bars can be especially important for this process, and therefore we discuss them in some detail below. Probably the most intriguing property of nested bars is their theoretically anticipated stage of a dynamical decoupling, when each bar exhibits a different pattern speed (Shlosman et al. 1989). Several aspects of this problem are analyzed below.

Gas Flow in Nested Bars. Here we assume that nested bars tumble with different pattern speeds, $`\mathrm{\Omega }_s>\mathrm{\Omega }_p`$, where ‘p’ stands for primary and ‘s’ — for secondary (sub-kpc) bar. Arguments dealing with chaos minimization in such time-dependent system will lead to certain limits, specifically that CR of the secondary bar must lie in the vicinity of the primary bar ILR, constraining $`\mathrm{\Omega }_s`$ (e.g., Pfenniger & Norman 1990). Such dynamical configuration, in principle, poses a problem for uninterrupting gas inflow towards smaller radii. According to this argument, the gas flow is repelled at the bar CR, inwards or outwards, because of the rim formed by the effective potential there, and hence may not cross the ILR/CR (i.e., bar-bar interface). Shlosman et al. (1989) have argued, in essence, that it is the gas self-gravity that overcomes repulsion by modifying the underlying potential. But even in the limit of neglibigle gravity in the gas, the flow is capable of crossing the bar-bar interface, although not in a steady manner and only for a restricted range of azimuthal angles (Shlosman & Heller 2002).

The pattern of shock dissipation in nested bars can be inferred from Fig. 1. It allows one to separate the incoming large-scale shocks from those driven by the secondary bar. Note that two systems of spiral shocks occur, each associated with the corresponding bar. The interaction between these shock systems shows detachment when the bars are perpendicular and attachment when they are aligned with each other. The shapes of the shocks depend on the angle between the bars.

Several factors characterize the gas dynamics in the decoupled bars: $`(i)`$ the time-dependent nature of the potential; $`(ii)`$ the nonsteady gas injection into the secondary bar which proceeds through the primary shocks penetrating the bar-bar interface. This phenomenon is absent at the CR of the primary bars. Unstable orbits in the interface region preclude the secondary bars from extending to their CR (El-Zant & Shlosman 2002, in preparation); $`(iii)`$ a fast-tumbling secondary bar which prevents the secondary ILRs from forming. Even in the case of a long-lived decoupled phase secondary bars are not expected to slow down. In fact, the gas inflow across the interface and the resulting central concentration can speed-up the bar (Heller, Noguchi, & Shlosman 1993, unpublished). The low-Mach-number gas flow is well organized and capable of following orbits within the bar with little dissipation. Non-linear orbit analysis reveals that, in the deep interior of the secondary bar, the $`x_1`$ orbits have a mild ellipticity and no end-loops. This result is robust. No offset large-scale shocks form under these conditions.

Knapen et al. (1995) have analyzed the shock dissipation in a self-consistent potential of ‘live’ stars and gas before the onset of decoupling, when both bars tumble with the same pattern speeds, and when the gas gravity is accounted for. No offset shocks have been found in this configuration either. These results indicate that gas inflows stagnate within the inner parts of fast rotating nuclear bars, forming nuclear disks if gas gravity is neglected — a condition similar to having a low surface density gas.

We conclude that compelling arguments show that no large-scale shocks and consequently no offset dust lanes will form inside secondary nuclear bars either when they are dynamically coupled and spin with the same pattern speeds as the primary bars, or dynamically decoupled, spinning much faster, if gas gravity is neglected. However, the fate of the gas settling inside the nuclear bars cannot be decided without invoking global gravitational effects in the gas which will completely change the nature of the flow.

Dynamical Decoupling of Nested Bars. Decoupling in nested bars is indirectly supported by the observed random orientation of primary and secondary bars. To complicate the matter, both bars can corotate, being completely synchronized. This configuration of nearly orthogonal bars may be a precursor to the future decoupled phase. The gas responding to the gravitational torques from the primary bar flows towards the center and encounters the $`x_2`$ orbits, which it populates. The forming secondary bar may be further strengthened by the gas gravity, or the amount of gas accumulating in the ILR resonance region may be insufficient to cause this runaway. The computational effort has so far gone into analyzing self-gravitating systems (e.g., Friedli 1999; Shlosman 1999).

If the secondary bar forms via gravitational instability (in stellar or gaseous disks), it will spin in the direction of the primary bar with $`\mathrm{\Omega }_s>\mathrm{\Omega }_p`$ (Shlosman et al. 1989; Friedli & Martinet 1993; Heller & Shlosman 1994). The presence of gas appears to be imperative for this to occur. Both bars are dynamically decoupled and the angle between them becomes arbitrary. The secondary bar can also rotate in the opposite sense to the primary bar, resulting from merging — a non-recurrent configuration.

Decoupling of Non-Self-Gravitating Bars. The actual degree of viscosity in the ISM is largely unknown. Heller, Shlosman & Englmaier (2001) have investigated the effect of viscosity on the gas settling in the ILR region of a single large-scale stellar bar. The most spectacular evolution occurred in the low-viscosity model, although all the models showed similar initial evolution, during which the gas accumulated in a double ring, corresponding roughly to two ILRs. In all the models the rings interacted hydrodynamically and merged (Fig. 2), forming a single oval-shaped ring corotating with the primary bar and leading it by $`\varphi _{\mathrm{dec}}`$ — depending on the gas viscosity. The remaining ring becomes increasingly oval and barlike (Fig. 2), its pattern speed changes abruptly, and it swings towards the primary bar, against the direction of rotation. In the inertial frame, this forming gaseous bar spins in the same sense as the main bar, albeit with $`\mathrm{\Omega }_s<\mathrm{\Omega }_p`$, but in the primary-bar frame it tumbles retrograde (!) for about 60 dynamical times, $``$ 2–3 $`\times 10^9`$ yr, until it is captured again by the primary. The shape of the decoupled bar and $`\mathrm{\Omega }_s`$ depend on bar orientation. The eccentricity, $`ϵ`$, reaches a maximum when bars are aligned. In standard and high-viscosity models, the secondary bar only librates about the primary.

The key to understanding this behavior lies in the distribution of gas particles with Jacobi energy, $`E_\mathrm{J}`$. After merging the ring is positioned close to the energy where the transition from $`x_2`$ to $`x_1`$ (at the inner ILR) occurs. The exact value of this transition energy is model-dependent, but this is of no importance to the essence of the decoupling. The crucial difference between the models comes from $`(i)`$ the position of the forming gaseous bar on the $`E_\mathrm{J}`$ axis after ring merging and $`(ii)`$ the value of $`\varphi _{\mathrm{dec}}`$.

As the gaseous bar forms at $`\varphi _{\mathrm{dec}}`$ to the primary, gravitational torques act to align the bars. In the low-viscosity model, gaseous bar resides on purely $`x_2`$ orbits and, therefore, responds to the torque by speeding up its precession while being pulled backwards, until it is almost at right angles to the bar potential valley. The decoupling happens abruptly when $`1/2`$ of the gas finds itself at $`E_\mathrm{J}`$ below the inner ILR. The absence of $`x_2`$ orbits at these $`E_\mathrm{J}`$ means that the gas loses its stable orientation along the primary-bar minor axis. The gaseous bar has a much smaller $`ϵ`$ in the fourth quadrant. Such an asymmetry with respect to the primary ensures that the torques from the primary bar are smaller in the fourth quadrant. But only for the least viscous model this becomes crucial, and the torques are unable to confine the bar oscillation, which continues for a full swing of $`2\pi `$. The nuclear bar is trapped again, after few rotations.

Correlation between $`ϵ`$ of the gaseous bar and its orientation can be tested observationally. Two additional effects should have observational consequence. First, the gas will cross the inner ILR on a dynamical timescale. Shlosman et al. (1989) pointed out that ILRs present a problem for radial gas inflow because the gas can stagnate there. A solution was suggested in the form of a global gravitational instability in the nuclear ring or disk, which will generate gravitational torques in the gas, driving it further in. Recently Sellwood & Moore (1999) resurrected the idea that ILRs would “choke” the gas inflow. Fig. 2 shows that even non-self-gravitating nuclear rings are prone to dynamical instability which drives the gas inward. Accounting for the gas gravity dramatically affects the stability of gas at the inner ILR, as discussed below (see also Fig. 3). This inflow is expected to be accompanied by star formation along the molecular bar, having a quasi-periodic, bursting character.

Decoupling of Self-Gravitating Bars. Gas gravity plays a crucial role in the formation and decoupling of bars, securing $`\mathrm{\Omega }_\mathrm{s}>\mathrm{\Omega }_\mathrm{p}`$. Simulations which tackle this issue have been very limited so far. Friedli & Martinet (1993, and priv. communication) and Combes (1994) experienced difficulties in decoupling pure stellar or gaseous bars. The inner bar in this case is very transient, which may be a result of an insuffient number of particles. For a mixed system the decoupled stage is more prolonged. Simulations explicitly demonstrate the necessity for gas to be present, and confirm that an increased central-mass concentration is important for the dynamical separation of the outer and inner parts. This can be achieved by moving the gas along the large-scale bar, accumulating it inside the ILR on the $`x_2`$ orbits, modifying the local potential and forming a double ILR. The gravity of the gas settling onto these orbits is sufficient to “drag” the stars along, but such a configuration still corotates with the primary. Settling onto the $`x_2`$ orbits depends upon gas sound speed and viscosity. When the gas is too viscous or hot, it will avoid the ILRs completely and remain on the $`x_1`$ orbits aligned with the bar. Moderately viscous gas will settle onto the innermost $`x_2`$ orbits. The evolution of nuclear regions in disks, therefore, depends on the (unknown) equation of state of the ISM. The inclusion of star formation in nuclear rings has already demonstrated how the resulting increase in viscosity leads to the mass transfer across the ILRs (Knapen et al. 1995).

Shlosman (2001) have described simulations which have accounted for gravitational effects in the gas. Ring evolution was similar to that without gas gravity, including the swing towards the primary bar. But this was followed by the rapid growth of prograde, $`\mathrm{\Omega }_s>\mathrm{\Omega }_p`$, self-gravitating modes with $`m=2`$ and 4 with the pattern speed much larger than $`\mathrm{\Omega }_\mathrm{p}`$, resulting in an avalanche-type inflow (Fig. 3). A model with the seed SBH revealed that inflowing gas feeds the SBH at peak rates, increasing its mass tenfold. This model shows that, at least theoretically, gaseous bars can dump the fuel at the very vicinity of the SBH.

## 3. Central Black Holes and Bulge-Halo Connection

Recent observations of galactic centers show clear correlations between the masses of the SBHs, $`M_{}`$, and masses and dispersion velocities of surrounding bulges, $`M_B`$ and $`\sigma _B`$, respectively. So far this is the most conclusive evidence that the very centers of galaxies are coupled in their evolution with the main bodies of their hosts. Analysis of the available data has led to the conclusion that the log $`M_{}\mathrm{log}\sigma _B`$ relation has a slope of 4 (Tremaine et al. 2002). El-Zant et al. (2002) have attempted to provide a physical explanation to this tantalizing phenomenon, based on the interaction between dark halos and baryonic components settling in their midst. In this formulation, the halo properties determine those of the bulge and the SBH.

The basic assumption is that of flat core triaxial halos subject to baryon infall (e.g., El-Zant, Shlosman & Hoffman 2001). The flat cores do not support loop orbits and stable gas motions. It is the baryon influx which perturbs the background harmonic core and leads to the appearance of progressively less oval loops. At the same time, the contracting baryonic component becomes self-gravitating. Interestingly, these two phenomena correlate (Fig. 4), so one expects that the onset of massive star formation coincides with termination of the dynamical infall, forming a bulge-like configuration. Unless the bulge is cuspy, the loop orbits are still absent deep inside, where the infall proceeds until the formation of a central SBH.

This settling of a baryonic matter in the mildly triaxial halos can in principle lead to the simultaneous formation of the SBH-bulge system with a characteristic bulge density, set by the density in the background halo. While the formation of round loop orbits in the vicinity of the SBH will choke the dynamical infall, the bulge will still grow, until a similar process will terminate this. The inner ‘bottleneck’ depends on the bulge-to-halo density ratio, $`\rho _B/\rho _H`$, while the outer bottleneck is governed by the bulge-to-halo mass ratio — limiting the SBH-to-bulge mass ratio $`K`$ to largely between $`10^410^2`$ (Fig. 5). Such a large spread in $`K`$ is indeed observed. The spread in $`M_{}\sigma _B`$ relation can be minimized if the halo follows the Faber-Jackson scaling (Faber & Jackson 1976).

The most important prediction from modeling this nonlinear process is that the bulge and SBH parameters are determined by the halo properties. Within this framework, larger and more massive halo cores are expected to produce, on average, larger and more massive bulges. Disks will form outside the halo core, or later, when the baryons destroy the core. If the core is very large, most of the baryonic material is consumed in the first phase and no significant disk forms. This effect is expected to be prominent in larger mass cores, since if these follow the Faber-Jackson relation, more massive halos should have proportionally larger cores. Other predictions include the requirement that the average density of the bulge in the central region should be close to that of the halo core. There is also a minimal bulge mass associated with a given core, although this varies significantly with the critical loop eccentricity assumed. A number of additional consequences for galaxy formation and evolution will be discussed elsewhere.

## 4. Conclusions

The subjects which have been emphasized in this review are nested bars and the SBH-bulge correlation. The nested bars are interesting for two main reasons, the possiblity to find them in the dynamically decoupled stage, and being a prime channel to fuel the AGN. Clearly, the ability of nested bars to trigger the gas inflow toward the very center depends on the gas fraction within the central kpc. Gaseous bars can dump the fuel at the very vicinity of the SBH. Neglecting the gravitational effects in this gas artificially terminates its subsequent evolution.

The origin of the SBH-bulge correlations can be the manisfestation of a global connection between the dark halos and the baryon-dominated luminous parts of galaxies. The flat-core triaxial halos provide the possible physical mechanism which can shape these inner parts, depending on the properties of dark matter component, in particular whether the latter exhibits the Faber-Jackson scaling. The halo shapes will have also a profound effect on the formation and subsequent growth of disk component and its stability to bar formation at higher redshifts.

### Acknowledgments.

I am grateful to my collaborators on the subjects discussed here, and to numerous colleagues for exciting discussions. Partially supported by NAG 5-10823, HST-AR-09546.01-A and AST-0206251 grants.

## References

Combes, F. 1994, Mass-Transfer Induced Activity in Galaxies, I. Shlosman, Ed. (Cambridge Univ. Press) 170

El-Zant, A., Shlosman, Begelman, M.C., Frank, J. 2002, ApJ, submitted

El-Zant, A., Shlosman, I., Hoffman, Y. 2001, ApJ, 760, 636

Faber, S.M., Jackson, R.E. 1976, ApJ, 204, 668

Friedli, D. 1999, Evolution of Galaxies on Cosmological Timescales, J.E. Beckman et al., Eds. (San Francisco: ASP), 88

Friedli, D., Martinet, L. 1993, A&A, 277, 2

Heller, C.H., Shlosman, I., Englmaier, P. 2001, ApJ, 553, 661

Heller, C.H., Shlosman, I. 1994, ApJ, 424, 84

Knapen, J.H., Beckman, J.E., Heller, C.H., Shlosman, I., de Jong, R.S. 1995, ApJ, 454, 623

Kormendy, J., Gebhardt, K. 2001, 20th Texas Symp. on Relativistic Astrophys., J.C. Wheeler & H. Martel, Eds. (AIP Conf. Proc.), 363

Laine, S., Shlosman, I., Knapen, J.H., Peletier, R.F. 2002, ApJ, 567, 97

Merritt, D., Ferrarese, L. 2001, The Central Kpc of Starbursts & AGN, J.H. Knapen et al., Eds. (ASP: San Francisco), 335

Pfenniger, D., Norman, C.A. 1990, ApJ, 363, 391

Sellwood, J.A., Moore, E.M. 1999, ApJ, 510, 125

Shlosman, I., Heller, C.H. 2002, ApJ, 565, 921

Shlosman, I. 2001, The Central Kpc of Starbursts & AGN, J.H. Knapen et al., Eds. (San Francisco), 55

Shlosman, I. 1999, Evolutiuon of Galaxies on Cosmological Timescales, J.E. Beckman et al., Eds. (San Francisco: ASP), 100

Shlosman, I., Frank, J., Begelman, M.C. 1989, Nat, 338, 45

Tremaine, S. et al. 2002, astro-ph/0203468