Volume 132, Issue 6, 14 February 2010

The Poisson–Boltzmann (PB) formalism is among the most popular approaches to modeling the solvation of molecules. It assumes a continuum model for water, leading to a dielectric permittivity that only depends on position in space. In contrast, the dipolar Poisson–Boltzmann–Langevin (DPBL) formalism represents the solvent as a collection of orientable dipoles with nonuniform concentration; this leads to a nonlinear permittivity function that depends both on the position and on the local electric field at that position. The differences in the assumptions underlying these two models lead to significant differences in the equations they generate. The PB equation is a second order, elliptic, nonlinear partial differential equation(PDE). Its response coefficients correspond to the dielectric permittivity and are therefore constant within each subdomain of the system considered (i.e., inside and outside of the molecules considered). While the DPBL equation is also a second order, elliptic, nonlinear PDE, its response coefficients are nonlinear functions of the electrostatic potential. Many solvers have been developed for the PB equation; to our knowledge, none of these can be directly applied to the DPBL equation. The methods they use may adapt to the difference; their implementations however are PBE specific. We adapted the PBE solver originally developed by Holst and Saied [J. Comput. Chem.16, 337 (1995)] to the problem of solving the DPBL equation. This solver uses a truncated Newton method with a multigrid preconditioner. Numerical evidences suggest that it converges for the DPBL equation and that the convergence is superlinear. It is found however to be slow and greedy in memory requirement for problems commonly encountered in computational biology and computational chemistry. To circumvent these problems, we propose two variants, a quasiNewton solver based on a simplified, inexact Jacobian and an iterative selfconsistent solver that is based directly on the PBE solver. While both methods are not guaranteed to converge, numerical evidences suggest that they do and that their convergence is also superlinear. Both variants are significantly faster than the solver based on the exact Jacobian, with a much smaller memory footprint. All three methods have been implemented in a new code named AQUASOL, which is freely available.
 ARTICLES

 Theoretical Methods and Algorithms

AQUASOL: An efficient solver for the dipolar Poisson–Boltzmann–Langevin equation
View Description Hide DescriptionThe Poisson–Boltzmann (PB) formalism is among the most popular approaches to modeling the solvation of molecules. It assumes a continuum model for water, leading to a dielectric permittivity that only depends on position in space. In contrast, the dipolar Poisson–Boltzmann–Langevin (DPBL) formalism represents the solvent as a collection of orientable dipoles with nonuniform concentration; this leads to a nonlinear permittivity function that depends both on the position and on the local electric field at that position. The differences in the assumptions underlying these two models lead to significant differences in the equations they generate. The PB equation is a second order, elliptic, nonlinear partial differential equation(PDE). Its response coefficients correspond to the dielectric permittivity and are therefore constant within each subdomain of the system considered (i.e., inside and outside of the molecules considered). While the DPBL equation is also a second order, elliptic, nonlinear PDE, its response coefficients are nonlinear functions of the electrostatic potential. Many solvers have been developed for the PB equation; to our knowledge, none of these can be directly applied to the DPBL equation. The methods they use may adapt to the difference; their implementations however are PBE specific. We adapted the PBE solver originally developed by Holst and Saied [J. Comput. Chem.16, 337 (1995)] to the problem of solving the DPBL equation. This solver uses a truncated Newton method with a multigrid preconditioner. Numerical evidences suggest that it converges for the DPBL equation and that the convergence is superlinear. It is found however to be slow and greedy in memory requirement for problems commonly encountered in computational biology and computational chemistry. To circumvent these problems, we propose two variants, a quasiNewton solver based on a simplified, inexact Jacobian and an iterative selfconsistent solver that is based directly on the PBE solver. While both methods are not guaranteed to converge, numerical evidences suggest that they do and that their convergence is also superlinear. Both variants are significantly faster than the solver based on the exact Jacobian, with a much smaller memory footprint. All three methods have been implemented in a new code named AQUASOL, which is freely available.

On the linear response and scattering of an interacting moleculemetal system
View Description Hide DescriptionA manybody Green’s function approach to the microscopic theory of plasmonenhanced spectroscopy is presented within the context of localized surfaceplasmon resonancespectroscopy and applied to investigate the coupling between quantummolecular and classicalplasmonic resonances in monolayercoated silver nanoparticles. Electronic propagators or Green’s functions, accounting for the repeated polarizationinteraction between a single molecule and its image in a nearby nanoscale metal, are explicitly computed and used to construct the linearresponse properties of the combined moleculemetal system to an external electromagnetic perturbation. Shifting and finite lifetime of states appear rigorously and automatically within our approach and reveal an intricate coupling between molecule and metal not fully described by previous theories. Selfconsistent incorporation of this quantummolecular response into the continuumelectromagnetic scattering of the moleculemetal target is exploited to compute the localized surfaceplasmon resonance wavelength shift with respect to the bare metal from first principles.

Statistical model for small clusters transforming from one isomer to another
View Description Hide DescriptionBased on the fact that the kinetic energy of one atom in small cluster still obeys the Boltzmann distribution, a statistical model is developed to predict the time consumed by a small cluster transforming from one isomer to another and is tested by vast molecular dynamics simulations of isomers transformation in helium gas at high temperatures (2000–3500 K). Extrapolating the model to lower temperatures, we found that the time for the most probable isomer of formed at 2500 K turning into the most stable one is more than at room temperature.

Assessment of an analytical density matrix derived from a modified Colle–Salvetti approach to the electron gas
View Description Hide DescriptionThe Ragot–Cortona model of local correlation energy [S. Ragot and P. Cortona, J. Chem. Phys.121, 7671 (2004)] revisits the initial approach of Colle and Salvetti [Theor. Chim. Acta37, 329 (1975)] in order to reinstate the kinetic contribution to the total correlation energy . In this work, the oneelectron reduced density matrix underlying the amended model is fully derived in closed form. By construction, the said density matrix is parameterfree but not representable, owing to approximations used in the Ragot–Cortona approach. However, the resulting density matrix is shown to have formally correct short and longrange expansions. Furthermore, its momentumspace counterpart qualitatively agrees with known parametrized momentum distributions except at small momenta, where the disagreement reflects the nonrepresentability of the model and restricts to a small fraction of the slowest electrons only.

A Chebychev propagator with iterative time ordering for explicitly timedependent Hamiltonians
View Description Hide DescriptionA propagation method for timedependent Schrödinger equations with an explicitly timedependent Hamiltonian is developed where time ordering is achieved iteratively. The explicit time dependence of the timedependent Schrödinger equation is rewritten as an inhomogeneous term. At each step of the iteration, the resulting inhomogeneous Schrödinger equation is solved with the Chebychev propagation scheme presented in the work of M. Ndong et al. [J. Chem. Phys.130, 124108 (2009)]. The iteratively timeordering Chebychev propagator is shown to be robust, efficient, and accurate and compares very favorably with all other available propagation schemes.

Protein solvation from theory and simulation: Exact treatment of Coulomb interactions in threedimensional theories
View Description Hide DescriptionSolvation forces dominate protein structure and dynamics. Integral equation theories allow a rapid and accurate evaluation of the effect of solvent around a complex solute, without the sampling issues associated with simulations of explicit solvent molecules. Advances in integral equation theories make it possible to calculate the angle dependent average solvent structure around an irregular molecular solution. We consider two methodological problems here: the treatment of longranged forces without the use of artificial periodicity or truncations and the effect of closures. We derive a method for calculating the longranged Coulomb interaction contributions to threedimensional distribution functions involving only a rewriting of the system of integral equations and introducing no new formal approximations. We show the comparison of the exact forms with those implied by the supercell method. The supercell method is shown to be a good approximation for neutral solutes whereas the new method does not exhibit the known problems of the supercell method for charged solutes. Our method appears more numerically stable with respect to thermodynamic starting state. We also compare closures including the Kovalenko–Hirata closure, the hypernettedchain (HNC) and an approximate threedimensional bridge function combined with the HNC closure. Comparisons to molecular dynamics results are made for water as well as for the proteinsolute bovine pancreatic trypsin inhibitor. The proposed equations have less severe approximations and often provide results which compare favorably to molecular dynamics simulation where other methods fail.

Alcohol solubility in a lipid bilayer: Efficient grandcanonical simulation of an interfacially active molecule
View Description Hide DescriptionWe derive a new densitybiased Monte Carlo technique which preserves detailed balance and improves the convergence of grandcanonical simulations of a species with a strong preference for an interfacial region as compared to the bulk. This densitybiasing technique is applied to the solubility of “alcohol” molecules in a mesoscopic model of the lipid bilayer, a system which has anesthetic implications but is poorly understood.

Genetic algorithm optimization of laser pulses for molecular quantum state excitation
View Description Hide DescriptionConventionally optimal control theory has been used in the theoreticaldesign of laser pulses through the direct variation in the electric field of the laser pulse as a function of time. This often leads to designed laser pulses which contain a broad and seemingly arbitrary frequency structure that varies in time in a manner which may be difficult to realize experimentally. In contrast, the experimental design of laser pulses has used a genetic algorithm (GA) approach, varying only those laser parameters actually available to the experimentalist. We investigate in this paper the possibility of using GA optimization methods in the theoreticaldesign of laser pulses to bring about quantum state transitions in molecules. This allows us to select only a small limited number of parameters to vary and to choose these parameters so that they correspond to those available to the experimentalist. In the paper we apply our methods to the vibrationalrotational excitation of the HF molecule. We choose a small limited number of frequencies and vary only the associated electric field amplitudes and pulse envelopes. We show that laser pulses designed in this way can lead to very high transition probabilities.

Quantitative prediction of gasphase and nuclear magnetic shielding constants
View Description Hide DescriptionHighlevel ab initio benchmark calculations of the and NMR chemical shielding constants for a representative set of molecules are presented. The computations have been carried out at the Hartree–Fock selfconsistent field (HFSCF), density functional theory(DFT) (BP86 and B3LYP), secondorder Møller–Plesset perturbation theory (MP2), coupled cluster singles and doubles (CCSD), and CCSD augmented by a perturbative treatment of triple excitations [CCSD(T)] level of theory using basis sets of triple zeta quality or better. The influence of the geometry, the treatment of electron correlation, as well as basis set and zeropoint vibrational effects on the shielding constants are discussed and the results are compared to gasphase experimental shifts. As for the first time a study using highlevel postHF methods is carried out for a secondrow element, we also propose a family of basis sets suitable for the computation of shielding constants. The mean deviations observed for and are 0.9 [CCSD(T)/13s9p4d3f] and −3.3 ppm [CCSD(T)/15s12p4d3f2g], respectively, when corrected for zeropoint vibrational effects. Results obtained at the DFT level of theory are of comparable accuracy to MP2 for and of comparable accuracy to HFSCF for . However, they are not improved by inclusion of zeropoint vibrational effects. The PN molecule is an especially interesting case with exceptionally large electron correlation effects on shielding constants beyond MP2 which, therefore, represents an excellent example for further benchmark studies.

Analytical theory of finitesize effects in mechanical desorption of a polymer chain
View Description Hide DescriptionWe discuss a unique system that allows exact analytical investigation of first and secondorder transitions with finitesize effects: mechanical desorption of an ideal lattice polymer chain grafted with one end to a solid substrate with a pulling force applied to the other end. We exploit the analogy with a continuum model and use accurate mapping between the parameters in continuum and lattice descriptions, which leads to a fully analytical partition function as a function of chain length, temperature (or adsorption strength), and pulling force. The adsorptiondesorption phase diagram, which gives the critical force as a function of temperature, is nonmonotonic and gives rise to reentrance. We analyze the chain length dependence of several chain properties (bound fraction, chain extension, and heat capacity) for different cross sections of the phase diagram. Close to the transition a single parameter (the product of the chain length and the deviation from the transition point) describes all thermodynamic properties. We discuss finitesize effects at the secondorder transition(adsorption without force) and at the firstorder transition (mechanical desorption). The firstorder transition has some unusual features: The heat capacity in the transition region increases anomalously with temperature as a power law, metastable states are completely absent, and instead of a bimodal distribution there is a flat region that becomes more pronounced with increasing chain length. The reason for this anomaly is the absence of an excess surface energy for the boundary between adsorbed and stretched coexisting phases (this boundary is one segment only): The two states strongly fluctuate in the transition point. The relation between mechanical desorption and mechanical unzipping of DNA is discussed.
 Gas Phase Dynamics and Structure: Spectroscopy, Molecular Interactions, Scattering, and Photochemistry

Theoretical study for the reaction of with
View Description Hide DescriptionThe lowlying triplet and singlet potential energy surfaces of the reaction have been studied at the level. On the triplet surface, six kinds of pathways are revealed, namely, direct hydrogen abstraction, Caddition/elimination, Naddition/elimination, substitution, insertion, and Hmigration. Multichannel Rice–Ramsperger–Kassel–Marcus theory and transitionstate theory are employed to calculate the overall and individual rate constants over a wide range of temperatures and pressures. It is predicted that the direct hydrogen abstraction and Caddition/elimination on triplet potential energy surface are dominant pathways. Major predicted end products include and . At atmospheric pressure with Ar and as bath gases, (IM1) formed by collisional stabilization is dominated at , whereas and NCO produced by Caddition/elimination pathway are the major products at the temperatures between 800 and 1500 K; the direct hydrogen abstraction leading to plays an important role at higher temperatures in hydrocarbon combustion chemistry and flames, with estimated contribution of 64% at 2000 K. Furthermore, the calculated rate constants are in good agreement with available experimental data over the temperature range 300–600 K. The kinetic isotope effect has also been calculated for the triplet reaction. On the singlet surface, the atomic oxygen can easily insert into C–H or C–C bonds of , forming the insertion intermediates and or add to the carbon atom of CN group in , forming the addition intermediate ; both approaches were found to be barrierless. It is indicated that the singlet reaction exhibits a marked difference from the triplet reaction. This calculation is useful to simulate experimental investigations of the reaction in the singlet state surface.

Disparate product distributions observed in reactions with and
View Description Hide DescriptionResults of gas phase reactivity studies on group six transition metal suboxide clusters,, , , and (, ; ) with both and are reported. Sequential oxidation for the more reduced species, , and dissociative addition for certain species, , is evident in the product distributions observed in mass spectrometric measurements. Reactions with proceed at a rate that is on the order of higher than for . The pattern of reaction products reveals compositiondependent chemical properties of these group six unary and binary clusters. At the core of this variation is the difference in Mo–O and W–O bond energies, the latter of which is significantly higher. This results in a larger thermodynamic drive to higher oxidation states in clusters with more tungsten atoms. However, addition products for more oxidized Wrich clusters are not observed, while they are observed for the more Morich clusters. This is attributed to the following: In the higher oxides (e.g., ), addition reactions require distortion of local metaloxygen bonding, and will necessarily have higher activation barriers for W–O bonds, since the vibrational potentials will be narrower. The binary clusters generally show sequential oxidation to higher values of . This again is attributed to higher W–O bond energy, the result being that stable binary structures have W atoms in higher oxidation states, and Mo centers both in more reduced states and sterically unhindered. The reduced Mo center provides a locus of higher reactivity. An unusual result that is not readily explained is the chemically inert behavior of .

Cooperative and diminutive hydrogen bonding in Y⋯HCN⋯HCN and NCH⋯Y⋯HCN trimers
View Description Hide DescriptionA computational study of the cooperative effect of hydrogen bonding in Y⋯HCN⋯HCN and its diminutive effect in NCH⋯Y⋯HCN linear complexes relative to the Y⋯HCN dimer was undertaken at the level of theory. It was found that the additional hydrogen bond in Y⋯HCN⋯HCN leads to an enhanced Y⋯HCN dissociation energy, extended H–C bond length, and larger redshift of the H–C stretch relative to Y⋯HCN, while opposite features are observed in NCH⋯Y⋯HCN. The cooperativity is diminished as the hardness of the Y atom directly bonded to the HCN molecule increases. A particularly interesting result is that the small bond contraction and blueshift associated with the H–C bond in BF⋯HCN is converted to a small bond extension and redshift on the formation of the BF⋯HCN⋯HCN trimer.

Argon nucleation in a cryogenic supersonic nozzle
View Description Hide DescriptionWe have measuredpressures and temperatures corresponding to the maximum nucleation rate of argon in a cryogenic supersonic nozzle apparatus where the estimated nucleation rates are . As increases from 34 to 53 K, increases from 0.47 to 8 kPa. Under these conditions, classical nucleationtheory predicts nucleation rates of 11–13 orders of magnitude lower than the observed rates while mean field kinetic nucleationtheory predicts the observed rates within 1 order of magnitude. The current data set appears consistent with the measurements of Iland et al. [J. Chem. Phys.127, 154506 (2007)] in the cryogenicnucleation pulse chamber. Combining the two data sets suggests that classical nucleationtheory fails because it overestimates both the critical cluster size and the excess internal energy of the critical clusters.

Exploring the mechanisms of H atom loss in simple azoles: Ultraviolet photolysis of pyrazole and triazole
View Description Hide DescriptionThe photophysics of gas phase pyrazole and 2H1,2,3triazole molecules following excitation at wavelengths in the range has been investigated using the experimental technique of H (Rydberg) atom photofragment translational spectroscopy. The findings are compared with previous studies of pyrrole and imidazole , providing a guide to H atom loss dynamics in simple Ncontaining heterocycles. CASPT2 theoretical methods have been employed to validate these findings. Photoexcitation of pyrazole at the longest wavelengths studied is deduced to involve excitation, but photolysis at is characterized by rapid N–H bond fission on a potential energy surface. The eventual pyrazolyl radical products are formed in a range of vibrational levels associated with both the ground and first excited electronic states as a result of nonadiabatic coupling at large N–H bond lengths. The excitation energy of the lowest state of pyrazole is found to be significantly higher in energy than that of pyrrole and imidazole. Similar studies of 2H1,2,3triazole reveal that the lowest state is yet higher in energy and not accessible following excitation at . The N–H bond strength of pyrazole is determined as , significantly greater than that of the N–H bonds in pyrrole and imidazole. The correlation between the photochemistry of azoles and the number and position of nitrogen atoms within the ring framework is discussed in terms of molecular symmetry and orbital electron density. A photodissociation channel yielding H atoms with low kinetic energies is also clearly evident in both pyrazole and 2H1,2,3triazole. Companion studies of pyrazole suggest that these slow H atoms arise primarily from the N–H site, following excitation, and subsequent internal conversion and/or unintended multiphoton absorption processes.

Probing the structural evolution of , , through a comparison of computed electron removal energies and experimental photoelectron spectra
View Description Hide DescriptionComputed electron removal energies for clusters, , are presented for the three lowestenergy isomers obtained from extensive, unbiased searches for the minimum energy structure at each size. The density functional theory(DFT) computations make use of a scheme introduced by Jellinek and Acioli (JA) [J. Chem. Phys.118, 7783 (2003)] that obtains electron removal energies from DFT orbital energies using corrections based on DFT total energies. The computed removal energies are compared with the measured photoelectron spectra (PES) for . The patterns of computed removal energies are shown to be isomer specific for clusters in this size range. By matching the computed removal energies to the observed PES, the isomers responsible for the PES are identified. The results of the JA scheme are compared to those obtained using other DFTbased methods.

Ab initio characterization of the Ca–HCl van der Waals complex
View Description Hide DescriptionThe equilibrium structure and threedimensional potential energy surface of the Ca–HCl van der Waals complex in its ground electronic state have been determined from accurate ab initio calculations using the coupledcluster method, CCSD(T), in conjunction with basis sets of quadruple and quintuplezeta quality. The coreelectron correlation, highorder valenceelectron correlation, and scalar relativistic effects were investigated. The Ca–HCl complex was confirmed to be linear at equilibrium, with the vibrationless dissociation energy (into Ca and HCl) of . The vibrationrotation energy levels of various Ca–HCl isotopomers were predicted using the variational method. The predicted spectroscopic constants can be useful in a further analysis of highresolution vibrationrotation spectra of the Ca–HCl complex.

The barrier height, unimolecular rate constant, and lifetime for the dissociation of
View Description Hide DescriptionAlthough never spectroscopically identified in the laboratory, hydrogenated nitrogen is thought to be an important species in combustion chemistry. The classical barrier height and exothermicity for the reaction are predicted by high level ab initio quantum mechanical methods [up to CCSDT(Q)]. Total energies are extrapolated to the complete basis set limit applying the focal point analysis. Zeropoint vibrational energies are computed using fundamental (anharmonic) frequencies obtained from a quartic force field. Relativistic and diagonal Born–Oppenheimer corrections are also taken into account. The quantum mechanical barrier with these corrections is predicted to be and the reaction exothermicity to be . The importance of these parameters for the thermal decomposition process is discussed. The unimolecular rate constant for dissociation of the molecule and its lifetime are estimated by canonical transitionstate theory and Rice–Ramsperger–Kassel–Marcus theory. The lifetime of the molecule is here estimated to be at room temperature. Our result is in marginal agreement with the latest experimental kinetic modeling studies , albeit consistent with the very rough experimental upper limit . For the dissociationreaction,kinetic isotope effects are investigated. Our analysis demonstrates that the molecule has a longer lifetime than the molecule. Thus, might be more readily identified experimentally. The ionization potential of the molecule is determined by analogous high level ab initio methods and focal point analysis. The adiabatic IP of is predicted to be , in only fair agreement with the experimental upper limit of 7.92 eV deduced from sychrothonradiationbased photoionization mass spectrometry.
 Condensed Phase Dynamics, Structure, and Thermodynamics: Spectroscopy, Reactions, and Relaxation

Spectroscopic investigation of OCS complexes inside helium droplets: Evidence for superfluid behavior
View Description Hide DescriptionUp to 16 parahydrogen and orthodeuterium molecules have been assembled around an OCS carbonyl sulfide chromophore molecule inside the pure and mixed droplets at temperatures of 0.38 and , respectively. The infrared spectra of the resulting complexes exhibit a sequence of rotationally resolved vibrational bands in the vicinity of , which are sufficiently separated to assign them to clusters with specific numbers of attached molecules for . The present article contains the first complete analysis of the spectra for and a full documentation of the results for briefly described in a short report [Europhys. Lett.83, 66008 (2008)]. Distinct rotational branches are observed for all OCS clusters at the He droplet temperatures of and , indicating that the shell rotates nearly freely about the molecular OCS axis. In the case of OCS at , the branch is seen for most , with the exception of , 6 and . At , the branch has disappeared for all , indicating that the axial rotations are no longer active. Previously, the absence of a branch for and 6 was explained by the high group symmetry of the bosonic rigid (donut) rings around the OCS molecule. This model, however, fails in explaining the disappearance of the branch for . In essential agreement with recent pathintegral Monte Carlo calculations, the observed phenomenon is attributed to the onset of superfluidity in the multiring shell and the related permutations of bosonic molecules. A floppy shell model, which accounts for the effect of tunneling and exchange of molecules within the clusters, is able to explain the postulated superfluid behavior of the shell at low temperatures. Within this model the activation of states of low axial symmetry is responsible for the appearance of the branch at higher temperatures.

Quantum instanton evaluations of surface diffusion, interior migration, and surfacesubsurface transport for H/Ni
View Description Hide DescriptionThe quantum instanton approximation is extended to investigate dynamical processes of hydrogen on surface, from surface to subsurface, and between interior sites in nickel lattice. The path integral Monte Carlo and adaptive umbrella sampling techniques are employed to manipulate the quantum instanton formula. The free energy profiles along reaction paths, temperature dependence of free energies, and rates as well as diffusion coefficients are calculated for each process. The results manifest that the motions of nickel atoms beneath the surface have little effect on the hydrogen diffusion on Ni(111), and the hydrogen at the fcc binding site is much easier to get into bulk nickel than the one at the hcp site. The temperature dependence of free energy profiles also reveals that the hydrogen in the subsurface octahedral vacancy and interior tetrahedral vacancy becomes unstable at low temperatures, which proposes a temperature dependence of reaction mechanism. In addition, the relaxations of the lattices dramatically lower the free energy barriers except for the process of the hydrogen diffusion on Ni(111). The quantum motions of the lattice atoms affect the free energies little at 300 K, but they hinder the rates by 20%–40% compared with the classical motions of lattice atoms.