Volume 141, Issue 8, 28 August 2014

Despite its importance in atmospheric science, much remains unknown about the microscopic mechanism of heterogeneous ice nucleation. In this work, we perform hybrid Monte Carlo simulations of the heterogeneous nucleation of ice on a range of generic surfaces, both flat and structured, in order to probe the underlying factors affecting the nucleation process. The structured surfaces we study comprise one basal plane bilayer of ice with varying lattice parameters and interaction strengths. We show that what determines the propensity for nucleation is not just the surface attraction, but also the orientational ordering imposed on liquid water near a surface. In particular, varying the ratio of the surface's attraction and orientational ordering can change the mechanism by which nucleation occurs: ice can nucleate on the structured surface even when the orientational ordering imposed by the surface is weak, as the water molecules that interact strongly with the surface are themselves a good template for further growth. We also show that lattice matching is important for heterogeneous nucleation on the structured surface we study. We rationalise these bruteforce simulation results by explicitly calculating the interfacial free energies of ice and liquid water in contact with the nucleating surface and their variation with surface interaction parameters.
 ARTICLES

 Theoretical Methods and Algorithms

Concentration effects on the rates of irreversible diffusioninfluenced reactions
View Description Hide DescriptionWe formulate a new theory of the effects of likeparticle interactions on the irreversible diffusioninfluenced bimolecular reactions of the type A + B → P + B by considering the evolution equation of the triplet ABB number density field explicitly. The solution to the evolution equation is aided by a recently proposed method for solving the Fredholm integral equation of the second kind. We evaluate the theory by comparing its predictions with the results of extensive computer simulations. The present theory provides a reasonable explanation of the simulation results.

Computing thermal Wigner densities with the phase integration method
View Description Hide DescriptionWe discuss how the Phase Integration Method (PIM), recently developed to compute symmetrized time correlation functions[M. Monteferrante, S. Bonella, and G. Ciccotti, Mol. Phys.109, 3015 (2011)], can be adapted to sampling/generating the thermal Wigner density, a key ingredient, for example, in many approximate schemes for simulating quantum time dependent properties. PIM combines a path integral representation of the density with a cumulant expansion to represent the Wigner function in a form calculable via existing Monte Carlo algorithms for sampling noisy probability densities. The method is able to capture highly nonclassical effects such as correlation among the momenta and coordinates parts of the density, or correlations among the momenta themselves. By using alternatives to cumulants, it can also indicate the presence of negative parts of the Wigner density. Both properties are demonstrated by comparing PIM results to those of reference quantum calculations on a set of model problems.

Validity conditions for moment closure approximations in stochastic chemical kinetics
View Description Hide DescriptionApproximations based on momentclosure (MA) are commonly used to obtain estimates of the mean molecule numbers and of the variance of fluctuations in the number of molecules of chemical systems. The advantage of this approach is that it can be far less computationally expensive than exact stochastic simulations of the chemical master equation. Here, we numerically study the conditions under which the MA equations yield results reflecting the true stochastic dynamics of the system. We show that for bistable and oscillatory chemical systems with deterministic initial conditions, the solution of the MA equations can be interpreted as a valid approximation to the true moments of the chemical master equation, only when the steadystate mean molecule numbers obtained from the chemical master equation fall within a certain finite range. The same validity criterion for monostable systems implies that the steadystate mean molecule numbers obtained from the chemical master equation must be above a certain threshold. For mean molecule numbers outside of this range of validity, the MA equations lead to either qualitatively wrong oscillatory dynamics or to unphysical predictions such as negative variances in the molecule numbers or multiple steadystate moments of the stationary distribution as the initial conditions are varied. Our results clarify the range of validity of the MA approach and show that pitfalls in the interpretation of the results can only be overcome through the systematic comparison of the solutions of the MA equations of a certain order with those of higher orders.

Symmetrical windowing for quantum states in quasiclassical trajectory simulations: Application to electron transfer
View Description Hide DescriptionIt has recently been shown [S. J. Cotton and W. H. Miller, J. Chem. Phys. 139, 234112 (2013)] that a symmetrical windowing quasiclassical (SQC) approach [S. J. Cotton and W. H. Miller, J. Phys. Chem. A 117, 7190 (2013)] applied to the MeyerMiller model [H.D. Meyer and W. H. Miller, J. Chem. Phys. 70, 3214 (1979)] for the electronic degrees of freedom in electronically nonadiabatic dynamics is capable of quantitatively reproducing quantum mechanical results for a variety of test applications, including cases where “quantum” coherence effects are significant. Here we apply this same SQC methodology, within a fluxside correlation function framework, to calculate thermal rate constants corresponding to several proposed models of electron transfer processes [P. Huo, T. F. Miller III, and D. F. Coker, J. Chem. Phys. 139, 151103 (2013); A. R. Menzeleev, N. Ananth, and T. F. Miller III, J. Chem. Phys. 135, 074106 (2011)]. Good quantitative agreement with Marcus Theory is obtained over several orders of magnitude variation in nonadiabatic coupling. Moreover, the “inverted regime” in thermal rate constants (with increasing bias) known from Marcus Theory is also reproduced with good accuracy by this very simple classical approach. The SQC treatment is also applied to a recent model of photoinduced proton coupled electron transfer [C. Venkataraman, A. V. Soudackov, and S. HammesSchiffer, J. Chem. Phys. 131, 154502 (2009)] and population decay of the photoexcited donor state is found to be in reasonable agreement with results calculated via reduced density matrix theory.

Stochastic manybody perturbation theory for anharmonic molecular vibrations
View Description Hide DescriptionA new quantum Monte Carlo (QMC) method for anharmonic vibrational zeropoint energies and transition frequencies is developed, which combines the diagrammatic vibrational manybody perturbation theory based on the Dyson equation with Monte Carlo integration. The infinite sums of the diagrammatic and thus sizeconsistent first and secondorder anharmonic corrections to the energy and selfenergy are expressed as sums of a few m or 2mdimensional integrals of wave functions and a potential energy surface (PES) (m is the vibrational degrees of freedom). Each of these integrals is computed as the integrand (including the value of the PES) divided by the value of a judiciously chosen weight function evaluated on demand at geometries distributed randomly but according to the weight function via the Metropolis algorithm. In this way, the method completely avoids cumbersome evaluation and storage of highorder force constants necessary in the original formulation of the vibrational perturbation theory; it furthermore allows even higherorder force constants essentially up to an infinite order to be taken into account in a scalable, memoryefficient algorithm. The diagrammatic contributions to the frequencydependent selfenergies that are stochastically evaluated at discrete frequencies can be reliably interpolated, allowing the selfconsistent solutions to the Dyson equation to be obtained. This method, therefore, can compute directly and stochastically the transition frequencies of fundamentals and overtones as well as their relative intensities as pole strengths, without fixednode errors that plague some QMC. It is shown that, for an identical PES, the new method reproduces the correct deterministic values of the energies and frequencies within a few cm^{−1} and pole strengths within a few thousandths. With the values of a PES evaluated on the fly at random geometries, the new method captures a noticeably greater proportion of anharmonic effects.

FranckCondon factors perturbed by damped harmonic oscillators: Solvent enhanced X ^{1}A_{g} ↔ A^{1}B_{1u} absorption and fluorescence spectra of perylene
View Description Hide DescriptionDamped harmonic oscillators are utilized to calculate FranckCondon factors within displaced harmonic oscillator approximation. This is practically done by scaling unperturbed Hessian matrix that represents local modes of force constants for molecule in gaseous phase, and then by diagonalizing perturbed Hessian matrix it results in direct modification of Huang–Rhys factors which represent normal modes of solute molecule perturbed by solvent environment. Scaling parameters are empirically introduced for simulating absorption and fluorescence spectra of an isolated solute molecule in solution. The present method is especially useful for simulating vibronic spectra of polycyclic aromatic hydrocarbon molecules in which hydrogen atom vibrations in solution can be scaled equally, namely the same scaling factor being applied to all hydrogen atoms in polycyclic aromatic hydrocarbons. The present method is demonstrated in simulating solvent enhanced X ^{1}Ag ↔ A^{1}B1u absorption and fluorescence spectra of perylene (mediumsized polycyclic aromatic hydrocarbon) in benzene solution. It is found that one of six active normal modes v 10 is actually responsible to the solvent enhancement of spectra observed in experiment. Simulations from all functionals (TD) B3LYP, (TD) B3LYP35, (TD) B3LYP50, and (TD) B3LYP100 draw the same conclusion. Hence, the present method is able to adequately reproduce experimental absorption and fluorescence spectra in both gas and solution phases.

Average local ionization energy generalized to correlated wavefunctions
View Description Hide DescriptionThe average local ionization energy function introduced by Politzer and coworkers [Can. J. Chem.68, 1440 (1990)] as a descriptor of chemical reactivity has a limited utility because it is defined only for onedeterminantal selfconsistentfield methods such as the Hartree–Fock theory and the Kohn–Sham densityfunctional scheme. We reinterpret the negative of the average local ionization energy as the average total energy of an electron at a given point and, by rewriting this quantity in terms of reduced density matrices, arrive at its natural generalization to correlated wavefunctions. The generalized average local electron energy turns out to be the diagonal part of the coordinate representation of the generalized Fock operator divided by the electron density; it reduces to the original definition in terms of canonical orbitals and their eigenvalues for onedeterminantal wavefunctions. The discussion is illustrated with calculations on selected atoms and molecules at various levels of theory.

Validity of virial theorem in allelectron mixed basis density functional, Hartree–Fock, and GW calculations
View Description Hide DescriptionIn this paper, we calculate kinetic and potential energy contributions to the electronic groundstate total energy of several isolated atoms (He, Be, Ne, Mg, Ar, and Ca) by using the local density approximation (LDA) in density functional theory, the Hartree–Fock approximation (HFA), and the selfconsistent GW approximation (GWA). To this end, we have implemented selfconsistent HFA and GWA routines in our allelectron mixed basis code, TOMBO. We confirm that virial theorem is fairly well satisfied in all of these approximations, although the resulting eigenvalue of the highest occupied molecular orbital level, i.e., the negative of the ionization potential, is in excellent agreement only in the case of the GWA. We find that the wave function of the lowest unoccupied molecular orbital level of noble gas atoms is a resonating virtual bound state, and that of the GWA spreads wider than that of the LDA and thinner than that of the HFA.

Reference hypernetted chain theory for ferrofluid bilayer: Distribution functions compared with Monte Carlo
View Description Hide DescriptionProperties of ferrofluid bilayer (modeled as a system of two planar layers separated by a distance h and each layer carrying a soft sphere dipolar liquid) are calculated in the framework of inhomogeneous OrnsteinZernike equations with reference hypernetted chain closure (RHNC). The bridge functions are taken from a soft sphere (1/r ^{12}) reference system in the pressureconsistent closure approximation. In order to make the RHNC problem tractable, the angular dependence of the correlation functions is expanded into special orthogonal polynomials according to Lado. The resulting equations are solved using the NewtonGRMES algorithm as implemented in the publicdomain solver NITSOL. Orientational densities and pair distribution functions of dipoles are compared with Monte Carlo simulation results. A numerical algorithm for the FourierHankel transform of any positive integer order on a uniform grid is presented.
 Atoms, Molecules, and Clusters

Relativistic coupled cluster study of diatomic compounds of Hg, Cn, and Fl
View Description Hide DescriptionThe structure and energetics of eight diatomic heavyatom molecules are presented. These include the species MAu, M2, and MHg, with M standing for the Hg, Cn (element 112), and Fl (element 114) atoms. The infiniteorder relativistic 2component Hamiltonian, known to closely reproduce 4component results at lower computational cost, is used as framework. Highaccuracy treatment of correlation is achieved by using the coupled cluster scheme with single, double, and perturbative triple excitations in large converged basis sets. The calculated interatomic separation and bond energy of Hg 2, the only compound with known experimental data, are in good agreement with measurements. The binding of Fl to Au is stronger than that of Cn, predicting stronger adsorption on gold surfaces. The bond in the M2 species is strongest for Fl2, being of chemical nature; weaker bonds appear in Cn2 and Hg 2, which are bound by van der Waals interactions, with the former bound more strongly due to the smaller van der Waals radius. The same set of calculations was also performed using the relativistic density functional theory approach, in order to test the performance of the latter for these weakly bound systems with respect to the more accurate coupled cluster calculations. It was found that for the MAu species the B3LYP functional provides better agreement with the coupled cluster results than the B88/P86 functional. However, for the M2 and the MHg molecules, B3LYP tends to underestimate the binding energies.

Alloying and oxidation of in situ produced coreshell Al@Yb nanoalloy particles—An “onthefly” study
View Description Hide DescriptionCoreshellstructured nanoalloy particles with an Aldominated interior covered by few Yb monolayers have been fabricated using a vaporaggregation method involving magnetron sputtering. The radially segregated structure of the YbAl nanoparticles has been disclosed by “onthefly” photoelectron spectroscopy monitoring of the nanoparticle beam in Yb 4f and Al 2p electron binding energy regions. Both, the binding energy values and the electron microscopy images taken on the deposited nanoparticles, allow estimating their dimensions to be in the 5–10 nm range. The photoelectron spectroscopy results suggest that in these nanoparticles no trivalent Yb – the typical case for the macroscopic YbAl alloy – is present. The oxidation of preformed YbAl nanoparticles was successfully attempted, leading to the appearance of divalent Yb surface oxide – in contrast to the bulk macroscopic Yb which is trivalent in the oxide. Our results suggest that at intermediate oxygen exposures “sandwichlike” nanoparticles of YbO/Yb/Al were synthesized. At higher O2 exposures, the oxygen seems to penetrate all the way to the YbAl interface. The results of the present study have to be considered when photonic applications of Ybdoped garnet nanoparticles are planned.

Chargetransfer excitations in lowgap systems under the influence of solvation and conformational disorder: Exploring rangeseparation tuning
View Description Hide DescriptionCharge transfer excitations play a prominent role in the fields of molecular electronics and light harvesting. At the same time they have developed a reputation for being hard to predict with timedependent density functional theory, which is the otherwise predominant method for calculating molecular structure and excitations. Recently, it has been demonstrated that rangeseparated hybrid functionals, in particular with an “optimally tuned” range separation parameter, describe chargetransfer excitations reliably for different molecules. Many of these studies focused on molecules in vacuum. Here we investigate the influence of solvation on the electronic excitations of thiophene oligomers, i.e., paradigm low gap systems. We take into account bulk solvation using a continuum solvation model and geometrical distortions from molecular dynamics. From our study, three main findings emerge. First, geometrical distortions increase absorption energies by about 0.5 eV for the longer thiophene oligomers. Second, combining optimal tuning of the range separation parameter with a continuum solvation method is not straightforward and has to be approached with great care. Third, optimally tuned rangeseparated hybrids without a shortrange exchange component tend to inherit undesirable characteristics of semilocal functionals: with increasing system size the range separation parameter takes a smaller value, leading to a functional of effectively more semilocal nature and thus not accurately capturing, e.g., the saturation of the optical gap with increasing system size.

High resolution electronic spectroscopy of the A ^{2}Σ^{−} − X ^{2}Π_{1/2} transition of PtN
View Description Hide DescriptionThe (2,0) vibrational band of the A ^{2}Σ^{−} − X ^{2}Π1/2 transition of platinum nitride, PtN, was recorded at Dopplerlimited resolution using intracavity laser absorption spectroscopy (ILS) and at subDoppler resolution using molecular beam laser induced fluorescence (LIF) spectroscopy. Isotopologue structure for ^{194}PtN, ^{195}PtN, and ^{196}PtN, magnetic hyperfine splitting due to ^{195}Pt (I = ½), and nuclear quadrupole splitting due to ^{14}N (I = 1) were observed in the spectrum. Molecular constants for the ground and excited states are derived. The hyperfine interactions are used to illuminate the nature of the A ^{2}Σ^{−} excited electronic state.

Nonadiabatic photofragmentation dynamics of BrCN^{−}
View Description Hide DescriptionThe photofragmentation dynamics of BrCN^{−} in the 270–355 nm and the 430–600 nm wavelength regions is explored both experimentally and theoretically. In the case of excitation between 430 nm and 600 nm, it is found that the molecular ion accesses two dissociation channels with a measured 60:40 branching ratio that is nearly constant over this range of photon energies. The dominant product channel corresponds to Br^{−} + CN, while the second channel correlates to spinorbit excited Br^{*} with CN^{−}. A larger wavelength dependence of the branching ratio is observed at shorter wavelengths, where the fraction of Br^{−} based products ranges from 80% to 95% at 355 nm and 270 nm, respectively. These branching ratios are reproduced and the mechanisms are explored by quantum dynamics calculations based on ground and excited state potential energy surfaces for BrCN^{−}, evaluated at the SOMRCISD level of theory. It is found that the electronic states that correlate to the two observed product channels are coupled through the spinorbit terms in the electronic Hamiltonian. The strength of this coupling displays a strong dependence on the BrCN angle. Specifically, after promotion to the excited state that is energetically accessible with 430–600 nm photons, it is found that when the wave packet accesses BrCN separations of between 4 Å and 6 Å, predominantly the Br^{−} + CN products are formed when the BrCN angle is smaller than 120°. For larger values of the BrCN angle, the Br^{*} + CN^{−} channel dominates. At the shorter wavelength excitation, the dynamics is complicated by a pair of states that correlate to electronically excited CN^{*} + Br^{−} products that borrow oscillator strength from the bright state, leading to an increase in the amount of Br^{−} relative to CN^{−}. The implications of these findings are discussed and compared to the experimentally measured product branching ratios for the photodissociation of BrCN^{−}.

Pressureinduced oligomerization of benzene at room temperature as a precursory reaction of amorphization
View Description Hide DescriptionOligomerization of benzene at high pressures up to 16 GPa was investigated at room temperature using an opposedanvil type pressure apparatus. The recovered samples were analyzed using GCMS to identify and quantify the products after the highpressure experiments. Some structural isomers of benzene dimer as well as biphenyl, naphthalene, and terphenyl isomers were detected at pressures higher than 13 GPa. The molar yield of the polycyclic aromatic hydrocarbons increased concomitantly with increasing pressure, although benzene still remained. The oligomerization is likely to occur when the neighbor distance of the benzene molecules exceeds the threshold of the reaction distance. The oligomerization is regarded as a precursory phenomenon of the amorphization that occurs at higher pressure.

Accurate double manybody expansion potential energy surface for the 2^{1} A′ state of
View Description Hide DescriptionAn accurate double manybody expansion potential energy surface is reported for the 2^{1} A′ state of . The new double manybody expansion (DMBE) form has been fitted to a wealth of ab initio points that have been calculated at the multireference configuration interaction level using the fullvalencecompleteactivespace wave function as reference and the ccpVQZ basis set, and subsequently corrected semiempirically via double manybody expansionscaled external correlation method to extrapolate the calculated energies to the limit of a complete basis set and, most importantly, the limit of an infinite configuration interaction expansion. The topographical features of the novel potential energy surface are then examined in detail and compared with corresponding attributes of other potential functions available in the literature. Exploratory trajectories have also been run on this DMBE form with the quasiclassical trajectory method, with the thermal rate constant so determined at room temperature significantly enhancing agreement with experimental data.
 Liquids, Glasses, and Crystals

Effects of surface interactions on heterogeneous ice nucleation for a monatomic water model
View Description Hide DescriptionDespite its importance in atmospheric science, much remains unknown about the microscopic mechanism of heterogeneous ice nucleation. In this work, we perform hybrid Monte Carlo simulations of the heterogeneous nucleation of ice on a range of generic surfaces, both flat and structured, in order to probe the underlying factors affecting the nucleation process. The structured surfaces we study comprise one basal plane bilayer of ice with varying lattice parameters and interaction strengths. We show that what determines the propensity for nucleation is not just the surface attraction, but also the orientational ordering imposed on liquid water near a surface. In particular, varying the ratio of the surface's attraction and orientational ordering can change the mechanism by which nucleation occurs: ice can nucleate on the structured surface even when the orientational ordering imposed by the surface is weak, as the water molecules that interact strongly with the surface are themselves a good template for further growth. We also show that lattice matching is important for heterogeneous nucleation on the structured surface we study. We rationalise these bruteforce simulation results by explicitly calculating the interfacial free energies of ice and liquid water in contact with the nucleating surface and their variation with surface interaction parameters.

The individual and collective effects of exact exchange and dispersion interactions on the ab initio structure of liquid water
View Description Hide DescriptionIn this work, we report the results of a series of density functional theory (DFT) based ab initio molecular dynamics (AIMD) simulations of ambient liquid water using a hierarchy of exchangecorrelation (XC) functionals to investigate the individual and collective effects of exact exchange (Exx), via the PBE0 hybrid functional, nonlocal van der Waals/dispersion (vdW) interactions, via a fully selfconsistent densitydependent dispersion correction, and an approximate treatment of nuclear quantum effects, via a 30 K increase in the simulation temperature, on the microscopic structure of liquid water. Based on these AIMD simulations, we found that the collective inclusion of Exx and vdW as resulting from a largescale AIMD simulation of (H2O)128 significantly softens the structure of ambient liquid water and yields an oxygenoxygen structure factor, S OO(Q), and corresponding oxygenoxygen radial distribution function, g OO(r), that are now in quantitative agreement with the best available experimental data. This level of agreement between simulation and experiment demonstrated herein originates from an increase in the relative population of water molecules in the interstitial region between the first and second coordination shells, a collective reorganization in the liquid phase which is facilitated by a weakening of the hydrogen bond strength by the use of a hybrid XC functional, coupled with a relative stabilization of the resultant disordered liquid water configurations by the inclusion of nonlocal vdW/dispersion interactions. This increasingly more accurate description of the underlying hydrogen bond network in liquid water also yields higherorder correlation functions, such as the oxygenoxygenoxygen triplet angular distribution, P OOO(θ), and therefore the degree of local tetrahedrality, as well as electrostatic properties, such as the effective molecular dipole moment, that are in much better agreement with experiment.

Molecular response of liquid nitrogen multiply shocked to 40 GPa
View Description Hide DescriptionLiquid nitrogen was subjected to multiple shock compression to examine its response to pressures (1540 GPa) and temperatures (18004000 K) previously unexplored in static and shock compression studies. Raman spectroscopy measurements were used to characterize the molecular bond response and to determine temperatures in the peak state. By extending our analysis to include other Raman spectroscopy measurements, an empirical relation was developed that describes the pressure and temperature dependence of the Raman shift (of the 2330 cm^{−1} mode) for both shock and static compression. Based on the PT dependence of the Raman shifts, the liquid nitrogen molecular response is best understood by considering three temperature regimes: below 1500 K, 15004000 K, and above 4000 K. For the pressures and temperatures accessed in the present work, liquid nitrogen remains a molecular fluid, and becomes a greybody emitter at the highest pressures.

Dipolar correlations in liquid water
View Description Hide DescriptionWe present an analysis of the dipolar correlations in water as a function of temperature and density and in the presence of simple ionic solutes, carried out using molecular dynamics simulations and empirical potentials. We show that the dipoledipole correlation function of the liquid exhibits sizable oscillations over nanodomains of about 1.5 nm radius, with several isosbestic points as a function of temperature; the size of the nanodomains is nearly independent on temperature and density, between 240 and 400 K and 0.9 and 1.3 g/cm^{3}, but it is substantially affected by the presence of solvated ions. In the same range of thermodynamic conditions, the decay time (τ) of the system dipole moment varies by a factor of about 30 and 1.5, as a function of temperature and density, respectively. At 300 K, we observed a maximum in τ as a function of density, and a corresponding shallow maximum in the tetrahedral order parameter, in a range where the diffusion coefficient, the pressure and the dielectric constant increase monotonically.