Volume 139, Issue 5, 07 August 2013
Index of content:

We measure the desorption of anions stimulated by the impact of 0–20 eV electrons on highly uniform thin films of plasmid DNAdiaminopropane. The results are accurately correlated with film thickness and composition by AFM and XPS measurements, respectively. Resonant structures in the H^{−}, O^{−}, and OH^{−} yield functions are attributed to the decay of transient anions into the dissociative electron attachment (DEA) channel. The diamine induces ammoniumphosphate bridges along the DNA backbone, which suppresses the DEA O^{−} channel and in counterpart increases considerably the desorption of OH^{−}. The close environment of the phosphate groups may therefore play an important role in modulating the rate and type of DNA damages induced by low energy electrons.
 COMMUNICATIONS


Communication: Favorable dimensionality scaling of rectangular collocation with adaptable basis functions up to 7 dimensions
View Description Hide DescriptionWe show that by using a rectangular collocation method with a small basis of parameterized functions, it is possible to compute a vibrational spectrum by solving the Schrödinger equation in 7D from a small number of ab initio calculations without a potential surface. The method is ideal for spectra of molecules adsorbed on a surface. In this paper, it is applied to calculate experimentally relevant energy levels of acetic acid adsorbed on the (101) surface of anatase TiO2. In this case, to obtain levels of experimental accuracy, increasing the number of dimensions from 4 to 7 increases the number of required potential points from about 1000 to about 10 000 and the number of basis functions from 126 to 792: the scaling is very attractive.

 ARTICLES

 Theoretical Methods and Algorithms

A convergent reactiondiffusion master equation
View Description Hide DescriptionThe reactiondiffusion master equation (RDME) is a lattice stochastic reactiondiffusion model that has been used to study spatially distributed cellular processes. The RDME is often interpreted as an approximation to spatially continuous models in which molecules move by Brownian motion and react by one of several mechanisms when sufficiently close. In the limit that the lattice spacing approaches zero, in two or more dimensions, the RDME has been shown to lose bimolecular reactions. The RDME is therefore not a convergent approximation to any spatially continuous model that incorporates bimolecular reactions. In this work we derive a new convergent RDME (CRDME) by finite volume discretization of a spatially continuous stochastic reactiondiffusion model popularized by Doi. We demonstrate the numerical convergence of reaction time statistics associated with the CRDME. For sufficiently large lattice spacings or slow bimolecular reaction rates, we also show that the reaction time statistics of the CRDME may be approximated by those from the RDME. The original RDME may therefore be interpreted as an approximation to the CRDME in several asymptotic limits.

Calculation of the “absolute” free energy of a βhairpin in an allatom force field
View Description Hide DescriptionWe propose a new approach to calculate the conformational free energy of a macromolecule in a compact stable state in implicit solvent. The free energy is evaluated with respect to an artificial reference system without internal interactions that is confined to a small welldefined multidimensional volume of a regular shape occupying approximately the same part of the conformational space as the macromolecule of interest. We present a practical implementation of our method, successfully apply it to a βhairpin in allatom representation, verify the results with direct parallel tempering simulations, and discuss the possibilities of further improvements.

A qubit coupled with confined phonons: The interplay between true and fake decoherence
View Description Hide DescriptionThe decoherence of a qubit coupled with the phonons of a finitesize lattice is investigated. The confined phonons no longer behave as a reservoir. They remain sensitive to the qubit so that the origin of the decoherence is twofold. First, a qubitphonon entanglement yields an incomplete true decoherence. Second, the qubit renormalizes the phonon frequency resulting in fake decoherence when a thermal average is performed. To account for the initial thermalization of the lattice, the qua ntum Langevin theory is applied so that the phonons are viewed as an open system coupled with a thermal bath of harmonic oscillators. Consequently, it is shown that the finite lifetime of the phonons does not modify fake decoherence but strongly affects true decoherence. Depending on the values of the model parameters, the interplay between fake and true decoherence yields a very rich dynamics with various regimes.

Orbitaloptimized coupledelectron pair theory and its analytic gradients: Accurate equilibrium geometries, harmonic vibrational frequencies, and hydrogen transfer reactions
View Description Hide DescriptionOrbitaloptimized coupledelectron pair theory [or simply “optimized CEPA(0),” OCEPA(0), for short] and its analytic energy gradients are presented. For variational optimization of the molecular orbitals for the OCEPA(0) method, a Lagrangianbased approach is used along with an orbital direct inversion of the iterative subspace algorithm. The cost of the method is comparable to that of CCSD [O(N 6) scaling] for energy computations. However, for analytic gradient computations the OCEPA(0) method is only half as expensive as CCSD since there is no need to solve the λ2amplitude equation for OCEPA(0). The performance of the OCEPA(0) method is compared with that of the canonical MP2, CEPA(0), CCSD, and CCSD(T) methods, for equilibrium geometries, harmonic vibrational frequencies, and hydrogen transfer reactions between radicals. For bond lengths of both closed and openshell molecules, the OCEPA(0) method improves upon CEPA(0) and CCSD by 25%–43% and 38%–53%, respectively, with Dunning's ccpCVQZ basis set. Especially for the openshell test set, the performance of OCEPA(0) is comparable with that of CCSD(T) (ΔR is 0.0003 Å on average). For harmonic vibrational frequencies of closedshell molecules, the OCEPA(0) method again outperforms CEPA(0) and CCSD by 33%–79% and 53%–79%, respectively. For harmonic vibrational frequencies of openshell molecules, the mean absolute error (MAE) of the OCEPA(0) method (39 cm−1) is fortuitously even better than that of CCSD(T) (50 cm−1), while the MAEs of CEPA(0) (184 cm−1) and CCSD (84 cm−1) are considerably higher. For complete basis set estimates of hydrogen transfer reaction energies, the OCEPA(0) method again exhibits a substantially better performance than CEPA(0), providing a mean absolute error of 0.7 kcal mol−1, which is more than 6 times lower than that of CEPA(0) (4.6 kcal mol−1), and comparing to MP2 (7.7 kcal mol−1) there is a more than 10fold reduction in errors. Whereas the MAE for the CCSD method is only 0.1 kcal mol−1 lower than that of OCEPA(0). Overall, the present application results indicate that the OCEPA(0) method is very promising not only for challenging openshell systems but also for closedshell molecules.

Treatment of scalarrelativistic effects on nuclear magnetic shieldings using a spinfree exacttwocomponent approach
View Description Hide DescriptionA costeffective treatment of scalarrelativistic effects on nuclear magnetic shieldings based on the spinfree exacttwocomponent theory in its oneelectron variant (SFX2C1e) is presented. The SFX2C1e scheme gains its computational efficiency, in comparison to the fourcomponent approach, from a focus on spinfree contributions and from the elimination of the small component. For the calculation of nuclear magnetic shieldings, the separation of spinfree and spindependent terms in the parent fourcomponent theory is carried out here for the matrix representation of the Dirac equation in terms of a restrictedmagnetically balanced gaugeincluding atomic orbital basis. The resulting spinfree fourcomponent matrix elements required to calculate nuclear magnetic shieldings are then used to construct the corresponding SFX2C1e Hamiltonian and its perturbed counterpart in the context of SFX2C1e analytic derivative theory. To demonstrate the applicability of the approach, we report coupledcluster calculations for prototypical problems such as the 17O shieldings of transitionmetal oxo complexes ( , M = Cr, Mo, and W) and the 129Xe shieldings of xenon fluorides (XeF2, XeF4, and XeF6).

Van der Waals interactions in density functional theory by combining the quantum harmonic oscillatormodel with localized Wannier functions
View Description Hide DescriptionWe present a new scheme to include the van der Waals (vdW) interactions in approximated Density Functional Theory (DFT) by combining the quantum harmonic oscillator model with the maximally localized Wannier function technique. With respect to the recently developed DFT/vdWWF2 method, also based on Wannier Functions, the new approach is more general, being no longer restricted to the case of well separated interacting fragments. Moreover, it includes higher than pairwise energy contributions, coming from the dipole–dipole coupling among quantum oscillators. The method is successfully applied to the popular S22 molecular database, and also to extended systems, namely graphite and H2 adsorbed on the Cu(111) metal surface (in this case metal screening effects are taken into account). The results are also compared with those obtained by other vdWcorrected DFT schemes.

A variational method for density functional theory calculations on metallic systems with thousands of atoms
View Description Hide DescriptionA new method for finitetemperature density functional theory calculations which significantly increases the number of atoms that can be simulated in metallic systems is presented. A selfconsistent, direct minimization technique is used to obtain the Helmholtz free energy of the electronic system, described in terms of a set of nonorthogonal, localized functions which are optimized in situusing a periodicsinc basis set, equivalent to plane waves. Most parts of the calculation, including the demanding operation of building the Hamiltonian matrix, have a computational cost that scales linearly with the number of atoms in the system. Also, this approach ensures that the Hamiltonian matrix has a minimal size, which reduces the computational overhead due to diagonalization, a cubicscaling operation that is still required. Large basis set accuracy is retained via the optimization of the localized functions. This method allows accurate simulations of entire metallic nanostructures, demonstrated with calculations on a supercell of bulk copper with 500 atoms and on gold nanoparticles with up to 2057 atoms.

Ad hoc methods for accurate determination of Bader's atomic boundary
View Description Hide DescriptionIn addition to the recently published triangulation method [P. M. Polestshuk, J. Comput. Chem.34, 206 (Year: 2013)]10.1002/jcc.23121, two new highly accurate approaches, ZFSX and SINTY, for the integration over an atomic region covered by a zeroflux surface (zfs) were developed and efficiently interfaced into the TWOE program. ZFSX method was realized as three independent modules (ZFSX1, ZFSX3, and ZFSX5) handling interatomic surfaces of a different complexity. Details of algorithmic implementation of ZFSX and SINTY are discussed. A special attention to an extended analysis of errors in calculations of atomic properties is paid. It was shown that uncertainties in zfs determination caused by ZFSX and SINTY approaches contribute negligibly (less than 10−6 a.u.) to the total atomic integration errors. Moreover, the new methods are able to evaluate atomic integrals with a reasonable time and can be universally applied for the systems of any complexity. It is suggested, therefore, that ZFSX and SINTY can be regarded as benchmark methods for the computation of any Quantum Theory of Atoms in Molecules atomic property.

Improving the accuracy and efficiency of timeresolved electronic spectra calculations: Cellular dephasing representation with a prefactor
View Description Hide DescriptionTimeresolved electronic spectra can be obtained as the Fourier transform of a special type of time correlation function known as fidelity amplitude, which, in turn, can be evaluated approximately and efficiently with the dephasing representation. Here we improve both the accuracy of this approximation—with an amplitude correction derived from the phasespace propagator—and its efficiency—with an improved cellular scheme employing inverse Weierstrass transform and optimal scaling of the cell size. We demonstrate the advantages of the new methodology by computing dispersed timeresolved stimulated emission spectra in the harmonic potential, pyrazine, and the NCO molecule. In contrast, we show that in strongly chaotic systems such as the quartic oscillator the original dephasing representation is more appropriate than either the cellular or prefactorcorrected methods.

The tensor hypercontracted parametric reduced density matrix algorithm: Coupledcluster accuracy with O(r^{4}) scaling
View Description Hide DescriptionTensor hypercontraction is a method that allows the representation of a highrank tensor as a product of lowerrank tensors. In this paper, we show how tensor hypercontraction can be applied to both the electron repulsion integral tensor and the twoparticle excitation amplitudes used in the parametric 2electron reduced density matrix (p2RDM) algorithm. Because only O(r) auxiliary functions are needed in both of these approximations, our overall algorithm can be shown to scale as O(r 4), where r is the number of singleparticle basis functions. We apply our algorithm to several small molecules, hydrogen chains, and alkanes to demonstrate its low formal scaling and practical utility. Provided we use enough auxiliary functions, we obtain accuracy similar to that of the standard p2RDM algorithm, somewhere between that of CCSD and CCSD(T).

Domain decomposition for implicit solvation models
View Description Hide DescriptionThis article is the first of a series of papers dealing with domain decomposition algorithms for implicit solvent models. We show that, in the framework of the COSMO model, with van der Waals molecular cavities and classical charge distributions, the electrostatic energy contribution to the solvation energy, usually computed by solving an integral equation on the whole surface of the molecular cavity, can be computed more efficiently by using an integral equation formulation of Schwarz's domain decomposition method for boundary value problems. In addition, the soobtained potential energy surface is smooth, which is a critical property to perform geometry optimization and molecular dynamics simulations. The purpose of this first article is to detail the methodology, set up the theoretical foundations of the approach, and study the accuracies and convergence rates of the resulting algorithms. The full efficiency of the method and its applicability to large molecular systems of biological interest is demonstrated elsewhere.

Permutation invariant polynomial neural network approach to fitting potential energy surfaces
View Description Hide DescriptionA simple, general, and rigorous scheme for adapting permutation symmetry in molecular systems is proposed and tested for fitting global potential energy surfaces using neural networks (NNs). The symmetry adaptation is realized by using loworder permutation invariant polynomials (PIPs) as inputs for the NNs. This socalled PIPNN approach is applied to the H + H2 and Cl + H2 systems and the analytical potential energy surfaces for these two systems were accurately reproduced by PIPNN. The accuracy of the NN potential energy surfaces was confirmed by quantum scattering calculations.

Interacting hard rods on a lattice: Distribution of microstates and density functionals
View Description Hide DescriptionWe derive exact density functionals for systems of hard rods with firstneighbor interactions of arbitrary shape but limited range on a onedimensional lattice. The size of all rods is the same integer unit of the lattice constant. The derivation, constructed from conditional probabilities in a Markov chain approach, yields the exact joint probability distribution for the positions of the rods as a functional of their density profile. For contact interaction (“sticky core model”) between rods, we give a lattice fundamental measure form of the density functional and present explicit results for contact correlators, entropy, free energy, and chemical potential. Our treatment includes inhomogeneous couplings and external potentials.

The accuracy of the GaussianandfiniteelementCoulomb (GFC) method for the calculation of Coulomb integrals
View Description Hide DescriptionWe analyze the accuracy of the Coulomb energy calculated using the GaussianandfiniteelementCoulomb (GFC) method. In this approach, the electrostatic potential associated with the molecular electronic density is obtained by solving the Poisson equation and then used to calculate matrix elements of the Coulomb operator. The molecular electrostatic potential is expanded in a mixed Gaussianfiniteelement (GF) basis set consisting of Gaussian functions of s symmetry centered on the nuclei (with exponents obtained from a full optimization of the atomic potentials generated by the atomic densities from symmetryaveraged restricted openshell Hartree–Fock theory) and shape functions defined on uniform finite elements. The quality of the GF basis is controlled by means of a small set of parameters; for a given width of the finite elements d, the highest accuracy is achieved at smallest computational cost when tricubic (n = 3) elements are used in combination with two (γH = 2) and eight (γ1st = 8) Gaussians on hydrogen and firstrow atoms, respectively, with exponents greater than a given threshold ( ). The error in the calculated Coulomb energy divided by the number of atoms in the system depends on the system type but is independent of the system size or the orbital basis set, vanishing approximately like d 4 with decreasing d. If the boundary conditions for the Poisson equation are calculated in an approximate way, the GFC method may lose its variational character when the finite elements are too small; with larger elements, it is less sensitive to inaccuracies in the boundary values. As it is possible to obtain accurate boundary conditions in linear time, the overall scaling of the GFC method for large systems is governed by another computational step—namely, the generation of the threecenter overlap integrals with three Gaussian orbitals. The most unfavorable (nearly quadratic) scaling is observed for compact, truly threedimensional systems; however, this scaling can be reduced to linear by introducing more effective techniques for recognizing significant threecenter overlap distributions.

Polarization as a field variable from molecular dynamics simulations
View Description Hide DescriptionA theoretical and computational framework for systematically calculating the macroscopic polarization density as a field variable from molecular dynamics simulations is presented. This is done by extending the celebrated Irving and Kirkwood [J. Chem. Phys.18, 817 (Year: 1950)10.1063/1.1747782] procedure, which expresses macroscopic stresses and heat fluxes in terms of the atomic variables, to the case of electrostatics. The resultant macroscopic polarization density contains molecular dipole, quadrupole, and higherorder moments, and can be calculated to a desired accuracy depending on the degree of the coarsegraining function used to connect the molecular and continuum scales. The theoretical and computational framework is verified by recovering the dielectric constant of bulk water. Finally, the theory is applied to calculate the spatial variation of the polarization vector in the electrical double layer of a 1:1 electrolyte solution. Here, an intermediate asymptotic length scale is revealed in a specific region, which validates the application of mean field PoissonBoltzmann theory to describe this region. Also, using the existence of this asymptotic length scale, the lengths of the diffuse and condensed/Stern layers are identified accurately, demonstrating that this framework may be used to characterize electrical double layers over a wide range of concentrations of solutions and surface charges.

When is the next extending of FickJacobs equation necessary?
View Description Hide DescriptionApplicability of the effective onedimensional equations, such as FickJacobs equation and its extensions, describing diffusion of particles in 2D or 3D channels with varying cross section A(x) along the longitudinal coordinate x, is studied. The leading nonstationary correction to ZwanzigRegueraRubí equation[R. Zwanzig, J. Phys. Chem.96, 3926 (Year: 1992)10.1021/j100189a004; D. Reguera and J. M. Rubí, Phys. Rev. E64, 061106 (Year: 2001)10.1103/PhysRevE.64.061106] is derived and tested on the exactly solvable model, diffusion in a 2D linear cone. The effects of such correction are demonstrated and discussed on elementary nonstationary processes, a time dependent perturbation of the stationary flow and calculation of the mean first passage time.

Multiscale enhanced path sampling based on the OnsagerMachlup action: Application to a model polymer
View Description Hide DescriptionWe propose a novel path sampling method based on the OnsagerMachlup (OM) action by generalizing the multiscale enhanced sampling technique suggested by Moritsugu and coworkers [J. Chem. Phys.133, 224105 (Year: 2010)10.1063/1.3510519]. The basic idea of this method is that the system we want to study (for example, some molecular system described by molecular mechanics) is coupled to a coarsegrained (CG) system, which can move more quickly and can be computed more efficiently than the original system. We simulate this combined system (original + CG system) using Langevin dynamics where different heat baths are coupled to the two systems. When the coupling is strong enough, the original system is guided by the CG system, and is able to sample the configuration and path space with more efficiency. We need to correct the bias caused by the coupling, however, by employing the Hamiltonian replica exchange, where we prepare many path replicas with different coupling strengths. As a result, an unbiased path ensemble for the original system can be found in the weakest coupling path ensemble. This strategy is easily implemented because a weight for a path calculated by the OM action is formally the same as the Boltzmann weight if we properly define the path “Hamiltonian.” We apply this method to a model polymer with AsakuraOosawa interaction, and compare the results with the conventional transition path sampling method.

Dynamics of a twolevel system under the simultaneous influence of a spin bath and a boson bath
View Description Hide DescriptionWe study dynamics of a twolevel system coupled simultaneously to a pair of dissimilar reservoirs, namely, a spin bath and a boson bath, which are connected via finite interbath coupling. It is found that the steadystate population transfer in the twolevel system increases with its coupling to the spin bath, while optimal transfer occurs at intermediate coupling in the transient process. If the twolevel system is strongly coupled to the spin bath, the population transfer is unidirectional barring minor population oscillations of minute amplitudes. If the spin bath is viewed as an atomic ensemble, robust generation of macroscopic superposition states exists against parameter variations of the twolevel system and the boson bath.
 Advanced Experimental Techniques

Experimental observation and analysis of the 3ν_{1}(Σ_{g}) stretching vibrational state of acetylene using continuouswave infrared stimulated emission
View Description Hide DescriptionWe present a sensitive experimental method for molecular spectroscopy that can be used to determine rovibrational states using midinfrared stimulated emission. Our infrared stimulated emission probing (IRSEP) experiment is based on using a narrowline, continuouswave Ti:sapphire laser beam (pump) to excite the molecules to an upper vibrational state and a continuouswave, midinfrared beam from an optical parametric oscillator (probe) to detect the stimulated emission by the excited molecules. Spectroscopic data are gathered by tuning the wavelengths of the beams. The molecules are probed before their velocity distribution is disturbed by collisions, which leads to a subDoppler resolution. The full width at half maximum of the emission peaks is below 10 MHz. The stimulated emission lines are measured with an accuracy of at least 0.005 cm−1. We use the IRSEP experiment to observe and analyze the symmetric rovibrational state [21+] (3ν 1(Σg)) of acetylene (C2H2). This state is not accessible via one photon transitions from the ground vibrational state. We use the leastsquares method to determine that the band center is at 9991.9725 (12) cm−1 and the rotational parameters are B = 1.156145 (22) and D = 1.608 (87) × 10−6 cm−1, where the uncertainties in parentheses are onestandard errors in the least significant digit.