Volume 131, Issue 6, 14 August 2009

We present a planewave basis set implementation of charge constrained density functionalmolecular dynamics (CDFTMD) for simulation of electron transferreactions in condensed phase systems. Following the earlier work of Wu and Van Voorhis [Phys. Rev. A72, 024502 (2005)], the density functional is minimized under the constraint that the charge difference between donor and acceptor is equal to a given value. The classical ion dynamics is propagated on the Born–Oppenheimer surface of the charge constrained state. We investigate the dependence of the constrained energy and of the energy gap on the definition of the charge and present expressions for the constraint forces. The method is applied to the electron selfexchange reaction in aqueous solution. Sampling the vertical energy gap along CDFTMD trajectories and correcting for finite size effects, a reorganization free energy of 1.6 eV is obtained. This is 0.1−0.2 eV lower than a previous estimate based on a continuum model for solvation. The smaller value for the reorganization free energy can be explained by the fact that the Ru–O distances of the divalent and trivalent Ru hexahydrates are predicted to be more similar in the electron transfer complex than for the separated aqua ions.
 COMMUNICATIONS


Propagation lengths of surface plasmon polaritons on metal films with arrays of subwavelength holes by infrared imaging spectroscopy
View Description Hide DescriptionMetalfilms with arrays of subwavelength holes (mesh) exhibit extraordinary transmission resonances to which many attribute a role for surface plasmonpolaritons (SPPs); others debated this point. Experimental measurements of propagation lengths are presented under conditions that pertain to the use of SPPs for surface spectroscopy. The lateral extent of electromagnetic propagation along the mesh surface is measured by recording absorption spectra of a line of latex microspheres as a function of distance away from the line along the mesh. Measurements reveal an exponential functional form for decay of absorption signal laterally from the absorption source. Results at , which are closest to the strongest transmission resonance of the mesh, reveal a propagation distance along the surface of . This is 40% larger than the lattice spacing implicating the holes as the SPP damping mechanism, however, this is significantly shorter than smooth metal expectations.

Cistrans switchable metallosupramolecular polymers: Computer modeling
View Description Hide DescriptionUsing computer simulations we study linear oligomers end functionalized with ligands that can form trans or cis2:1 complexes with metal ions in a saltscreened good solvent. We show that transcisisomerization of ligandmetal complexes can significantly increase the average molecular weight as well as trigger formation of reversible metallosupramolecular network based on 3:1 ligandmetal complexes acting as crosslinkers. We predict the conditions under which the most dramatic changes in the properties of metallosupramolecular polymers, such as network formation or increase in elastic plateau modulus of the network, occur upon isomerization.
 Top

 ARTICLES

 Theoretical Methods and Algorithms

Charge constrained density functional molecular dynamics for simulation of condensed phase electron transfer reactions
View Description Hide DescriptionWe present a planewave basis set implementation of charge constrained density functionalmolecular dynamics (CDFTMD) for simulation of electron transferreactions in condensed phase systems. Following the earlier work of Wu and Van Voorhis [Phys. Rev. A72, 024502 (2005)], the density functional is minimized under the constraint that the charge difference between donor and acceptor is equal to a given value. The classical ion dynamics is propagated on the Born–Oppenheimer surface of the charge constrained state. We investigate the dependence of the constrained energy and of the energy gap on the definition of the charge and present expressions for the constraint forces. The method is applied to the electron selfexchange reaction in aqueous solution. Sampling the vertical energy gap along CDFTMD trajectories and correcting for finite size effects, a reorganization free energy of 1.6 eV is obtained. This is 0.1−0.2 eV lower than a previous estimate based on a continuum model for solvation. The smaller value for the reorganization free energy can be explained by the fact that the Ru–O distances of the divalent and trivalent Ru hexahydrates are predicted to be more similar in the electron transfer complex than for the separated aqua ions.

Timereversible molecular dynamics algorithms with bond constraints
View Description Hide DescriptionTimereversible molecular dynamics algorithms with bond constraints are derived. The algorithms are stable with and without a thermostat and in double precision as well as in singleprecision arithmetic. Time reversibility is achieved by applying a centraldifference expression for the velocities in the expression for Gauss’ principle of least constraint. The imposed time symmetry results in a quadratic expression for the Lagrange multiplier. For a system of complex molecules with connected constraints the corresponding set of coupled quadratic equations is easily solved by a consecutive iteration scheme. The algorithms were tested on two models. One is a dumbbell model of Toluene, the other system consists of molecules with four connected constraints forming a triangle and a branch point of constraints. The equilibrium particle distributions and the meansquare particle displacements for the dumbbell model were compared to the corresponding functions obtained by GROMACS. The agreement is perfect within statistical error.

Efficient and accurate approximations to the local coupled cluster singles doubles method using a truncated pair natural orbital basis
View Description Hide DescriptionA production level implementation of the closedshell local quadratic configuration interaction and coupled cluster methods with single and double excitations (QCISD and CCSD) based on the concept of pair natural orbitals [local pair natural orbital LPNOQCISD and LPNOCCSD) is reported, evaluated, and discussed. This work is an extension of the earlier developed LPNO coupledelectron pair approximation (LNPOCEPA) method [F. Neese et al., Chem. Phys.130, 114108 (2009)] and makes extended use of the resolution of the identity (RI) or density fitting (DF) approximation. Two variants of each method are compared. The less accurate approximations still recover 98.7%–99.3% of the correlation energy in the given basis and have modest disk space requirements. The more accurate variants typically recover 99.75%–99.95% of the correlation energy in the given basis but require the Coulomb and exchange operators with up to twoexternal indices to be stored on disk. Both variants have comparable computational efficiency. The convergence of the results with respect to the natural orbital truncation parameter has been studied. Extended numerical tests have been performed on absolute and relative correlation energies as function of basis set size and as well as on reaction energies, isomerization energies, and weak intermolecular interactions. The results indicate that the errors of the LPNO methods compared to the canonical QCISD and CCSD methods are below 1 kcal/mol with our default thresholds. Finally, some calculations on larger molecules are reported (ranging from 40–86 atoms) and it is shown that for medium sized molecules the total wall clock time required to complete the LPNOCCSD calculations is only two to four times that of the preceding selfconsistent field (SCF). Thus these methods are highly suitable for largescale computational chemistry applications. Since there are only three thresholds involved that have been given conservative default values, the methods can be confidentially used in a “blackbox” fashion in the same way as their canonical counterparts.

Construction of basis sets for timedependent studies
View Description Hide DescriptionThe common basis sets constructed for use in electronic structure calculations have been found inadequate for the representation of electrons participating in nonadiabatic timedependent dynamics calculations. In this paper we outline an approach to construct electronic bases better suited for dynamical processes such as energy deposition and charge transfer in binary collisions of ions, atoms, and molecules. Since electrons of manyatom systems commonly are represented by orbitals formed as linear combinations of atomic orbitals, the focus is on atomic basis sets. The main idea is to construct basis sets that adequately reproduce the first few excitation energies of neutral atoms. In this paper we outline a method for such basis set construction of various levels of accuracy for firstrow atoms and give a few illustrative examples.

Fullconfigurationinteraction calculation of threebody nonadditive contribution to helium interaction potential
View Description Hide DescriptionThe threebody nonadditive interactionenergy between helium atoms was calculated at 253 trimer configurations using the fullconfigurationinteraction (FCI) method. The analytic potential fitted to these energies is the best current representation of the threebody nonadditive interactions between helium atoms. At the equilateral triangle configuration with , near the minimum of the total potential, the nonadditive threebody energy calculated at the FCI level amounts to −88.5 mK, compared to −98.5 mK at the coupled cluster with single, double, and noniterative triple excitations [CCSD(T)] level. The uncertainty of the former result resulting from basis set incompleteness is estimated to be 1.5 mK. The relative uncertainty of our present complete threebody fit, including the uncertainties resulting from the fitting procedure, is estimated at 2%, a fivefold improvement over the previous best potential. Overall, the FCI contribution beyond CCSD(T) is rather important, being of the same order of magnitude as the uncertainty of the sum of twobody interactions. The inclusion of this contribution makes uncertainties of the total trimer interactionenergies dominated by the uncertainties of the twobody component.

Transient hydrodynamical behavior by dynamical nonequilibrium molecular dynamics: The formation of convective cells
View Description Hide DescriptionWe present a method based on dynamical nonequilibrium molecular dynamics (DNEMD) that allows one to produce rigorous ensemble averages for the transient regimes. We illustrate the method by describing the formation of convective cells within a twodimensional fluid system of soft disks in which a gravity field and a thermal gradient are present. We analyze two different physical settings, with the thermal gradient orthogonal or parallel to the gravity field. In both settings, we follow the formation of the convective flows from the initial time, when the perturbation is turned on, to the steady state. In the first setting (orthogonal fields) we investigate several different cases, varying the initial stationary ensemble and the perturbing field. We find that the final steadystate convective cell is independent of the specific sequence of perturbation fields, which only affects the transient behavior. In all cases, we find that the convective roll is formed through a sequence of damped oscillations of the local fields (density, temperature, and velocity), superimposed to an overall relaxation toward the local steadystate values. Then, we show how DNEMD can be applied to the Rayleigh–Bénard (RB) setting (parallel fields). In these conditions, the convective flow only establishes above a threshold, without a preferred verse of rotation. We analyze only the response to the ignition of the gravity field in a stationary system under the action of a vertical thermal gradient. Also in this case we characterize the transient response by following the evolution of the density, temperature, and velocity fields until the steadystate RB convective cell is formed. The observed transients are similar to those observed in the case of orthogonal fields. However, the final steady states are quite different. Finally, we briefly discuss the conditions for the general applicability of the DNEMD method.

On the use of shifted Jacobi polynomials in accurate evaluation of roots and weights of Rys polynomials
View Description Hide DescriptionIn this paper it is shown that shifted Jacobi polynomials can be used in connection with the Gaussian quadrature modified moment technique to greatly enhance the accuracy of evaluation of Rys roots and weights used in Gaussian integral evaluation in quantum chemistry. A general fourterm inhomogeneous recurrence relation is derived for the shifted Jacobi polynomial modified moments over the Rys weight function . It is shown that for this general fourterm inhomogeneous recurrence relation reduces to a threeterm dependent inhomogeneous recurrence relation. Adjusting to proper values depending on the Rys exponential parameter , the method is capable of delivering highly accurate results for large number of roots and weights in the most difficult to treat intermediate range. Examples are shown, and detailed formulas together with practical suggestions for their efficient implementation are also provided.

Efficient exact and skip methods for stochastic simulation of coupled chemical reactions
View Description Hide DescriptionGillespie’s direct method (DM) [D. Gillespie, J. Chem. Phys.81, 2340 (1977)] for exact stochastic simulation of chemical reactionsystems has been widely adopted. It is easy to implement but requires large computation for relatively large systems. Recently, two more efficient methods, next reaction method (NRM) [M. A. Gibson and J. Bruck, J. Phys. Chem. A105, 1876 (2000)] and optimized DM (ODM) [Y. Cao et al., J. Chem. Phys.121, 4059 (2004)], have been developed to improve simulation speed. It has been demonstrated that the ODM is the stateoftheart most efficient method for exact stochastic simulation of most practical reactionsystems. In this paper, we first develop an exact stochastic simulation algorithm named ODMK that is more efficient than the ODM. We then develop an approximate method named skip method to further accelerate simulation. Using two chemical reactionsystems, we demonstrate that our ODMK and skip method can save 20%–30% and 70%–80% simulation time, respectively, comparing to the ODM. We also show that our ODMK and skip method provide almost the same simulation accuracy as the ODM.

Analytic gradients for the statespecific multireference coupled cluster singles and doubles model
View Description Hide DescriptionThe general theory of analytic energy gradients is presented for the statespecific multireference coupled cluster method introduced by Mukherjee and coworkers [Mol. Phys.94, 157 (1998)], together with an implementation within the singles and doubles approximation, restricted to two closedshell determinants and Hartree–Fock orbitals. Expressions for the energy gradient are derived based on a Lagrangian formalism and cast in a densitymatrix notation suitable for implementation in standard quantumchemical program packages. In the present implementation, we exploit a decomposition of the multireference coupled cluster gradient expressions, i.e., lambda equations and the corresponding density matrices, into a socalled singlereference part for each reference determinant and a coupling term. Our implementation exhibits the proper scaling, i.e., with as the number of reference determinants and as the number of orbitals, and it is thus suitable for largescale applications. The applicability of our multireference coupled cluster gradients is illustrated by computations for the equilibrium geometry of the 2,6isomers of pyridyne and the pyridynium cation. The results are compared to those from singlereference coupled cluster calculations and are discussed with respect to the future perspectives of multireference coupled clustertheory.

Optimizing the hyperpolarizability tensor using external electromagnetic fields and nuclear placement
View Description Hide DescriptionWe investigate the effects of an external electric and magnetic field on the first hyperpolarizability tensor of a quantum system, such as a molecule or nanoparticle, whose nonlinear response is well below the fundamental limit. We find that the intrinsic hyperpolarizability is optimized when the applied electric and magnetic fields are comparable to the internal molecular fields. Indeed, the nonlinear response is just as large for an electron in the presence of the external field without the nuclei as it is for an electron bound to a molecule and in the presence of the applied field. We find that all combinations of fields and molecular structures that optimize the largest diagonal component of the intrinsic hyperpolarizability share the same universal properties: The threelevel ansatz is obeyed, the normalized transition moment to the dominant state is about 0.76, the ratio of the two dominant excited state energies is about 0.48, the electron density tends toward being onedimensional, and the intrinsic hyperpolarizability is less than 0.71. Thus, strategies for optimizing the hyperpolarizability should focus on ways to achieve these universal properties. On the other hand, when is optimized, the three level ansatz appears to hold for a pair of degenerate states. In this case, the energy ratio between the pairs of degenerate states is 0.42 and the normalized transition moment to the pair of dominant states is 0.87. Most importantly, the intrinsic hyperpolarizability is 0.9, the largest ever calculated for a system described by a potential energy function.
 Gas Phase Dynamics and Structure: Spectroscopy, Molecular Interactions, Scattering, and Photochemistry

The Stark effect in Rydberg states of a highly polar diatomic molecule: CaF
View Description Hide DescriptionThe Stark effect in molecular Rydberg states is qualitatively different from the Stark effect in atomic Rydberg states because of the anisotropy of the ion core and the existence of rotational and vibrational degrees of freedom. These uniquely molecular features cause the electricfieldinduced decoupling of the Rydberg electron from the body frame to proceed in several stages in a molecule. Because the transitiondipole moment among the sameRydberg states is much larger than the permanent dipole moment of the ion core, the decoupling of the Rydberg electron from the ion core proceeds gradually. In the first stage, analyzed in detail in this paper, and are mixed by the external electric field, while is conserved. In the further stages, as the external electric field increases, , , and are expected to undergo mixing. We have characterized these stages in , states of CaF. The large permanent dipole moment of makes CaF qualitatively different from the other molecules in which the Stark effect in Rydberg states has been described (, , , NO, and ) and makes it an ideal testbed for documenting the competition between the external and dipole electric fields. We use the weakfield Stark effect to gain access to the lowest rotational levels of , , and states and to assign their actual or nominal quantum number. Lowest rotational levels provide information needed to disentangle the shortrange and longrange interactions between the Rydberg electron and the ion core. We diagonalize an effective Hamiltonian matrix to determine the characters of the corenonpenetrating states and to characterize their mixing with the corepenetrating states. We conclude that the mixing of the , state with lower states is stronger than documented in our previous multichannel quantum defect theory and longrange fits to zerofield spectra.

Timedependent quantum wavepacket description of H and D atom tunneling in N–H and N–D photodissociation of methylamine and methylamine
View Description Hide DescriptionThe degree to which tunneling through a barrier in the N–H and N–D photodissociation channels of methylamine and its deuterated variant , respectively, plays a role was investigated by timedependent quantum wavepacket dynamics calculations. Two dimensional potential energy surfaces (PESs) of methylamine, presenting the N–H stretch and the HNC bend, were constructed employing multireference ab initio electronicstructure methods, allowing full description of the H motion on the plane. The timedependent Schrödinger equation was solved employing the Fourier method for calculating the Hamiltonian operation together with the Chebychev polynomial expansion of the evolution operator. The results show that tunneling and decay to vibrational resonant states on the first excited electronic PES are faster for the H atom than for the D. The decay into two of the resonant states found on the first PES strongly depends on the initially excited vibrational state on the ground electronic PES.

A model Hamiltonian to simulate the complex photochemistry of benzene II
View Description Hide DescriptionThe photophysics and photochemistry of benzene is a classic example of the richness of competing pathways available to a molecule after photoexcitation. Computer simulations are one way to provide a molecular picture for the dynamics behind the experimental observations. In this paper we develop a vibronic coupling Hamiltonian prepared in a previous paper [G. A. Worth, J. Photochem. Photobiol., A190, 190 (2007)]. Using CASPT2 we add dynamic correlation to the description of the excited states, improving their accuracy dramatically. Seven coupled states and all vibrational modes are included in the model and the parameters are obtained by fitting to points provided by the quantum chemistry calculations. The model is shown to be a good fit of the adiabatic surfaces and its accuracy is demonstrated by the calculation of three absorption bands, which compare favorably with the experimentally obtained spectra.

Photoelectron spectroscopy of small cluster anions
View Description Hide DescriptionWe report the photoelectron spectra of small cluster anions . The vibrational stateresolved spectrum of permits reliable identification of the origins of the excited and states of neutral IBr through a highquality Franck–Condon spectral simulation. As a result, we directly determine several important spectroscopic parameters: the adiabatic electron affinity (EA) of IBr, , the ground electronic state bond strength of , , its equilibrium bond length, , and its vibrational frequency, . These values represent a substantial improvement over existing experimental information and are in good agreement with recent theoretical studies. The photoelectron spectra of the first three cluster anions, , do not exhibit resolved vibrational structure, but the similarity to the photoelectron spectrum indicates minimal electron delocalization onto the solvent. The cluster anion spectra shift to progressively higher electron binding energies, providing information on the magnitude of the solvent perturbation and estimates of the EA of .

Gas phase electronic spectrum of Tshaped radical
View Description Hide DescriptionGas phase electronic transitions for the and band systems of Tshaped radical have been measured in the 345–475 nm range. Vibrational analyses of both band systems are reported. Simulation of several rotationally resolved bands confirms previously obtained rotational parameters for the state. The radical is produced by ablating an aluminum rod in the presence of acetylene gas. The resulting supersonic molecular beam is probed using both massselective resonant twocolor twophoton ionization and laser induced fluorescence.Ab initio calculations and vertical electronic excitation energies help the assignment. Vibrational frequencies for the , , and states have been determined. Rotational analysis of a number of bands yields spectroscopic constants for one vibronic state in the manifold and the origin band of the system.

A systematic search for minimum structures of small gold clusters and their electronic properties
View Description Hide DescriptionA systematic search for global and energetically lowlying minimum structures of neutral gold clusters is performed within a seeded genetic algorithm technique using density functional theory together with a relativistic pseudopotential. Choosing the energetically lowest lying structures we obtain electronic properties by applying a larger basis set within an energyconsistent relativistic smallcore pseudopotential approach. The possibility of extrapolating these properties to the bulk limit for such small cluster sizes is discussed. In contrast to previous calculations on cesium clusters [B. Assadollahzadeh et al., Phys. Rev. B78, 245423 (2008)] we find a rather slow convergence of any of the properties toward the bulk limit. As a result, we cannot predict the onset of metallic character with increasing cluster size, and much larger clusters need to be considered to obtain any useful information about the bulk limit. Our calculated properties show a large oddeven cluster size oscillation in agreement, for example, with experimental ionization potentials and electron affinities. For the calculated polarizabilities we find a clear transition to lower values at , the first cluster size where the predicted global minimum clearly shows a compact threedimensional (3D) structure. Hence, the measurement of cluster polarizabilities is ideal to identify the transition at low temperatures for gold. Our genetic algorithm confirms the pyramidal structure for .

and radicalmolecule van der Waals complex
View Description Hide DescriptionOH and SH radicals are important in atmospheric chemistry because of their high reactivity. We examine the Van der Waals radicalmolecule complexes formed by OH and SH radicals with molecular nitrogen. The van der Waals radicalmolecule complex between OH and CO, which is isoelectronic to OH and , is also examined as a calibration of the computational results to literature experimental findings. In this work, we employ high level ab initio methods to investigate the stability and spectroscopic properties of these complexes. Natural bond analysis is also performed in order to study their bonding features.

Influence of the carrier gas molar mass on the particle formation in a vapor phase
View Description Hide DescriptionThe influence of the molar mass of a carrier gas on the formation of nanoparticles in the vapor phase is investigated. The function of the carrier gas atmosphere is the regulation of the particle temperature by collisions with the cluster surface. The aim of this work is to optimize the carrier gas in a simulation in order to mimic a large amount of carrier gas atoms by few gas atoms with effective parameters. In this context the efficiency of the heat exchange with the carrier gas depending on its molar mass is analyzed. As a result one finds for varying molar masses and unchanged interaction parameters a competition between the efficiency and the number of the collisions. For too small molar masses the energy exchange per collision is too small while for too high masses the carrier gas atoms become very slow, decreasing the number of collisions.