The exact simulation of quantum mechanical systems on classical computers generally scales exponentially with the size of the system N. Using quantum computers, the computational resources required to carry out the simulation are polynomial. Our group has been working in the development and characterization of quantum computational algorithms for the simulation of chemical systems. We will give a tutorial on our algorithms for the simulation of molecular electronic structure, molecular properties and quantum dynamics, and will discuss the opportunities, open questions and challenges in the field of simulation of physical systems using quantum computers or dedicated quantum devices.

The application of conventional GGA, and meta-GGA, density functionals to van der Waals complexes is fraught with difficulties. Conventional functionals do not contain the physics of the dispersion interaction. To make matters worse, the exchange part alone can yield anything from severe over “binding” to severe over repulsion depending on the choice of functional. We rectify these problems by

a) adding a dispersion term with nonempirical C6, C8, and C10 dispersion coefficients (the Becke-Johnson dispersion model),

b) selecting a GGA exchange functional (PW86, also nonempirical) that gives excellent agreement with exact Hartree-Fock exchange repulsion curves.

The result is a simple GGA+dispersion theory giving excellent noble-gas pair interaction energies for He through Kr with only two adjustable parameters in the dispersion cutoff.

Studies of molecular dynamics and molecular spectroscopy generally start from the Born-Oppenheimer approximation and require some form of analytical potential energy surface fitted to ab initio electronic structure calculations. We have used computational invariant theory and the MAGMA computer algebra system as an aid to develop representations for the potential energy and dipole moment surfaces that are fully invariant under permutations of like nuclei. We express the potential energy surface in terms of internuclear distances using basis functions that are manifestly invariant. A dipole moment is represented with use of effective charges at positions of the nuclei, which must transform as a covariant, rather than as an invariant, under permutations of like nuclei. Malonaldehyde (CHOHCHCHO) provides an illustrative application. The associated molecular permutational symmmetry group is of order 288 (4!3!2!) and the use of full permutational symmetry makes it possible to obtain a compact representation for the surface.

We explored the relative stability, structure, and conductance of crossed nanotube junctions with dispersion corrected density functional theory. We found that the most stable junction geometry, not studied before, displays the smallest conductance. While the conductance increases as force is applied, it levels off very rapidly. This behavior contrasts with a less stable junction geometry that show steady increase of the conductance as force is applied. Electromechanical sensing devices based on this effect should exploit the conductance changes close to equilibrium.

Recent work in my group has focussed on the semiclassical origins of density functional theory, and how much of modern DFT can be understood in these terms, including the limitations of present approximations. I will discuss this in detail for model systems, describing a method that avoids DFT altogether. This leads to a grand algorithmic challenge, whose solution could revolutionize electronic structure calculations, by allowing much larger numbers of electrons to be tackled.

Recent advances in treating electrons and nuclei (typically protons) quantum mechanically without the Born-Oppenheimer approximation using nuclear-electronic orbital (NEO) method and multicomponent density functional theory (MCDFT) is presented. Electron-proton dynamical correlation is highly significant because of the attractive electrostatic interaction between the electron and the proton. Inadequate treatment of electron-proton correlation produces nuclear densities that are too localized, resulting in abnormally high stretching frequencies, as well as inaccuracies in thermally averaged geometries and isotope effects. To address this problem, an explicitly correlated Hartree-Fock (NEO-XCHF) scheme has been formulated to include explicit electron-proton correlation directly into the nuclear-electronic orbital self-consistent-field framework. This approach is based on a general ansatz for the nuclear-electronic wavefunction in which explicit dependence on the electron-proton distance is incorporated into the total wavefunction using Gaussian-type geminal functions. A multicomponent density functional theory (MCDFT) has also been formulated by developing electron-proton functionals based on the explicitly correlated ansatz for the nuclear-electronic wavefunction. Benchmark calculations illustrate that these new methods significantly improve the description of the nuclear densities, thereby leading to more accurate hydrogen vibrational frequencies and vibrationally averaged geometries.

Density functional theory of electronic structure is widely and successfully applied in simulations throughout engineering and sciences. However, for many predicted properties there are spectacular failures that can be traced to the delocalization error and static correlation error of commonly used approximations. These errors can be characterized and understood through the perspective of fractional charges and fractional spins introduced recently. Reducing these errors will open new frontiers for applications of density functional theory.

Joint work with Reinhart Ahlrichs.

Density functional theory provides a systematic approach to the electronic structure of atoms, molecules and solids. It requires the repeated solution of single particle Schrodinger equations in a self consistent loop. Most techniques involve some sort of basis set, the most common ones being plane waves or Gaussians. In crystalline materials the most accurate solutions involve augmented basis sets. These combine numerical solutions of the Schrodinger equation in regions near the atomic nucleii with so called ‘tail functions’ in more distant regions. In the linear augmented plane wave (LAPW) method the tail functions are plane waves. This formulation has been incorporated into the WIEN2k code. With the current interest in nanoscale clusters, biomolecules, and other finite systems it is desirable to have a comparably accurate method for these. While it is always possible to build supercells, it is often convenient to have completely localized functions which eliminate interaction between periodic images. We recently proposed a finite cluster version of the linear augmented Slater-type orbital (LASTO) method [1]. STO’s have the correct behavior at large distances and possess an addition theorem – they can be re-expanded about other sites with analytic coefficients. We solve the Poisson equation by replacing the spherical part of the density near the nucleii with a smooth pseudo-density. The full potential, including the non-sphrical piece is then solved on a grid. Examples of small clusters and comparison with the Gaussian based program NWChem will be given.

[1] K. S. Kang, J. W. Davenport, J. Glimm, D. E. Keyes, and M. McGuigan, submitted to J. Computational Chemistry.

We have implemented relativisitic GVVPT2 using DKH integrals and ANO-RCC basis sets from Molcas package. It is done by developing an interface code accessing and transforming one- and two-electron integral array, making it available for any method implemented within the UNDMol package. The relativistic GVVPT2 is applied to calculations of ground and excited states potential energy curves of TiC and CrH.

Joint work with Alexander Gaenko and Jochen Autschbach.

We develop a time-dependence quasirelativistic density functional ther based on the ZORA approximation for computing frequency dependent linear response of molecules. Density fitting was used for the calculation of complex components of the frequency dependent dipole-dipole polarizability. CPKS equations based on 2-componenent ZORA response were derived. Using damping techniques excitation energy corresponding to the poles of the polarizability curves were calculated. We present the results of the calculations of complex dipole-dipole polarizability, of two and three dimensional gold clusters, and absorption spectra of heavy metal oxides.

Spectroscopic studies of the hydroxymethyl and 1-hydroxyethyl radicals have found an unusually large difference in their ionization energies (IE). The anticipated decrease in the IE of the latter radical due to its larger size does not fully account for the experimentally observed difference of 0.92 eV. Here we investigated the problem with the aid of electronic structure calculations. We found that the large drop in the IE of 1-hydroxyethyl radical is a result of the combined effects of the destabilization of its highest occupied molecular orbital and the stabilization of the corresponding cation due to hyperconjugation. This qualitative explanation agrees with a simple Huckel-like approach and is also consistent with Natural Bond Orbital calculation results.

We are interested in realizing the variational calculation with the 2-order Reduced Density Matrix (2-RDM) for fermionic systems. The known necessary N-representability conditions P, Q, G, T1, and T2' are imposed, resulting in an optimization problem called semidefinite programming (SDP) problem. The ground state energies for various small atoms and molecules were calculated. Additionally, we show some results of the one-dimensional Hubbard model for high correlation limit using multiple precision arithmetic version of the solver.

We present a new linear scaling method for electronic structure computations in the context of Kohn-Sham density functional theory (DFT). The method is based on a subspace iteration, and takes advantage of the non-orthogonal formulation of the Kohn-Sham functional, and the improved localization properties of non-orthogonal wave functions. We demonstrate the efficiency of the algorithm for practical applications by performing fully three-dimensional computations of the electronic density of alkane chains.

This is joint work with Jianfeng Lu, Yulin Xuan, and Weinan E, at Princeton University.

The Resolution of the Identity (RI) is widely used in many-body algorithms. It expresses the completeness of a set of functions that possess the familiar orthonormality property.

In the first part of my lecture, I will discuss functions that possess an analogous property called Coulomb-orthonormality and which permit us to resolve the two-particle Coulomb operator into a sum of products of one-particle functions. Connections to Cholesky decomposition and Kronecker-product approximation will be made.

In the second part of the lecture, I will present and discuss numerical applications of Coulomb resolution in the context of electronic structure theory.

The green fluorescent protein (GFP) is widely used in biochemical and medical studies as a biomarker in living cells. Modeling properties of GFP is an essential step in the efforts to enhance efficiency of this in vivo marker and to expand the area of its applications. We apply the methods of electronic structure calculations, including the quantum chemistry methods and the combined quantum mechanical – molecular mechanical (QM/MM) approaches, to describe the structure, spectra and transformations of the GFP chromophore, 4-hydroxybenzylidene-imidazolinone, in the gas phase, solutions and in the protein matrix. We compare the results of calculations for the cis-trans chromophore isomerization in the ground electronic state in the gas phase by using the DFT approach PBE0/6-31+G** and the CASSCF(12,11)/cc-pVDZ approximation. The energy profiles computed with both methods are markedly different in the vicinity of the saddle point. The isomerization paths computed in the QM/MM approach for the chromophore buried inside the water cluster show that the CASSCF results are better consistent with the experimental observations than the DFT findings. We also report the results of QM/MM calculations with the DFT approximations in the QM subsystem for the geometry configurations of the chromophore binding pocket inside the protein by assuming various protonation states of the chromophore unit, anionic, neutral, zwitterionic, cationic. These structures are employed for calculations of the photoexcitation pathways in GFP.

First-principles molecular dynamics (FPMD) is a simulation method that combines molecular dynamics with the accuracy of a quantum mechanical description of electronic structure. It is increasingly used to address problems of structure determination, statistical mechanics, and electronic structure of solids, liquids and nanoparticles. The high computational cost of this approach makes it a good candidate for use on large-scale computers. In order to achieve high performance on terascale and petascale computers, current FPMD algorithms have to be reexamined and redesigned. We present new, large-scale parallel algorithms developed for FPMD simulations on computers including O(10³) to O(10⁴) CPUs. Examples include the problem of simultaneous diagonalization of symmetric matrices used in the calculation of Maximally Localized Wannier Functions (MLWFs), and the Orthogonal Procrustes problem that arises in the context of Born-Oppenheimer molecular dynamics simulations.

Supported by NSF-OCI PetaApps through grant 0749217.

A new generation of density-functional methods is based on orbital-dependent funtionals. With orbital-dependent functionals long-standing problems like the occurence of unphysical Coulomb self-interactions or the qualitatively wrong description of charge-transfer excitation in time-dependent density-functional theory can be solved. Orbital-dependent functionals indeed may represent the future of density-functional theory. In oder to determine exchange-correlation potentials corresponding to orbital-dependent energy functionals, however, it is necessary to solve a numerically very demanding integral equation with the optimized effective potential (OEP) method. Methods to handle this integral equation thus are required for the further development of density-functional methods. The numerical stability of the OEP integral equation is investigated and method to solve it are presented.

We study convergence or divergence of the Moller–Plesset perturbation series for systems with two electrons and a single nucleus of charge Z > 0. This question is essentially to determine if the radius of convergence of a power series in the complex perturbation parameter lambda is greater than 1. We examine a simple one-dimensional model with delta functions in place of Coulomb potentials and the realistic three-dimensional model. For each model, we show rigorously that if the nuclear charge Z is sufficiently large, there are no singularities for real values of lambda between -1 and 1. Using a finite difference scheme, we present numerical results for the delta function model.

Wave function-based quantum chemistry has two traditional lines of development – one based on molecular orbitals (MO's), and the other on valence bond (VB) theory. Both offer advantages and disadvantages for the challenging problem of describing strong correlations, such as the breaking of chemical bonds, or the low-spin (antiferromagnetic) coupling of electrons on different centers.

Within MO methods, strong correlations can be viewed as those arising within a valence orbital active space. One reasonable definition of such a space is to supply one correlating orbital for each valence occupied orbital. Exact solution of the Schrodinger equation in this space is exponentially difficult with its size, and therefore approximations are imperative. The most common workaround is to truncate the number of orbitals defining the active space, and then solve the truncated problem, as is done in CASSCF. An important alternative is to systematically approximate the Schrödinger equation in the full valence space, for example by using coupled cluster theory ideas. I shall discuss progress in this direction.

Within spin-coupled VB theory, the target wave function consists of a set of non-orthogonal orbitals, one for each valence electron, that are spin-coupled together into a state of the desired overall spin-multiplicity. The number of active orbitals is identical with the valence space MO problem discussed above, though the problem is not identical. Exact solution of the VB problem is exponentially difficult with molecular size, and therefore approximations are imperative. Again, the most common approach is to seek the exact solution in a truncated valence orbital space, where other orbitals are simply treated in mean-field. It is possible, however, to also consider approximations that do not truncate the space, but rather reduce the complexity. A new way of doing this will be introduced and contrasted with the MO-based approaches.

We perform a rigorous mathematical analysis of the bending modes of a linear triatomic molecule that exhibits the Renner-Teller effect. Assuming the potentials are smooth, we prove that the wave functions and energy levels have asymptotic expansions in powers of epsilon, where the fourth power of epsilon is the ratio of an electron mass to the mass of a nucleus. To prove the validity of the expansion, we must prove various properties of the leading order equations and their solutions. The leading order eigenvalue problem is analyzed in terms of a parameter b, which is equivalent to the parameter originally used by Renner. Perturbation theory and finite difference calculations suggest that there is a crossing involving the ground bending vibrational state near b=0.925. The crossing involves two states with different degeneracy.

A recent algorithmic revision of second order Generalized van Vleck perturbation theory (GVVPT2) has proven to make the method efficacious for many challenging molecular systems. 1 An extension to third order (GVVPT3) has been demonstrated to be a close approximation to multireference configuration interaction including single and double excitations (MRCISD). 2 To improve the computing efficiency, new GVVPT codes have been developed to take advantage of recently implemented configuration-driven configuration interaction (CI) with unitary group approach (UGA).

Joint work with Jerzy Leszczynski, Computational Center for Molecular Structure and Interactions, Jackson State University, Jackson MS and Leonid Gorb, US Army ERDC, Vicksburg, MS.

With the June announcement that RoadRunner supercomputer is the first system to reach the petaflop level, the HPC community is entering a realm of unprecedented computing power. More petascale computing systems will soon be available to the scientific community. Recent studies in the productivity of HPC platforms point to better software as a key enabler to science on these systems.

The combination of computationally demanding electronic structure methods with molecular dynamics is highly dependent on high-performance computing resources. The availability of such applications constitutes a big opportunity to evaluate both capabilities and limits of any HPC system and software application within the framework of a real-life feasibility study. The performance of benchmarks from the AIMD and hybrid QM/MM simulations on two high performance computing platforms will be discussed. Looking toward maximizing the computational time/performance ratio, we analyzed performance data for the Cray XT3/XT4 architectures available at ERDC.

The difficultly of approximate density functionals in describing the energetics of Diels-Alder reactions and dimerization of aluminum complexes is analyzed. Both of these reaction classes involve formation of cyclic or bicyclic products, which are found to be under-bound by the majority of functionals considered. We present a consistent view of these results from the perspective of delocalization error. This error causes approximate functionals give too low energy for delocalized densities or too high energy for localized densities, as in the cyclic and bicyclic reaction products. This interpretation allows us to understand better a wide range of errors in main-group thermochemistry obtained with popular density functionals. In general, functionals with minimal delocalization error should be used for theoretical studies of reactions where there is a loss of extended conjugation or formation of highly branched, cyclic, and cage-like molecules.

The bimolecular interaction potentials for various configurations of the ethylene dimer computed with coupled-cluster and spin-component scaled MP2 are reported. Of particular interest is any bias for particular orientations of the sigma- and pi-bonds introduced by the scaling of correlation components.

Over the last decade the coupled-cluster (CC) methodology played a dominant role in highly accurate predictions of electronic structure. For this reason, the need for more efficient parallel implementations is obvious. Currently, the existing NWChem implementation can scale across thousand of CPUs and can be used in correlating 250 electrons. Additionally, the CC codes can be used as a quantum mechanical component of various multiscale approaches.

Joint work with Marat Valiev, Niri Govind, Peng-Dong Fan, W.A. de Jong (William R Wiley Environmental Molecular Sciences Laboratory and Chemical Sciences Division, Pacific Northwest National Laboratory P.O. Box 999, MS K1-96, Richland, WA 99352) and Jeff R. Hammond (The University of Chicago).

The coupled-cluster (CC) methodology has become a leading formalism not only in gas-phase calculations but also in modeling systems for which the inclusion of the surrounding environment is critical for a comprehensive understanding of complex photochemical reactions. At the same time it has been proven that high-level CC formalisms are capable of providing highly adequate characterization of excitation energies and excited-state potential energy surfaces. With the ever increasing power of computer platforms and highly scalable codes, very accurate QM/MM calculations for large molecules (defining the quantum region) can be routinely performed in the foreseeable future even with iterative methods accounting for the effect of triples ( CCSDT-n/EOMCCSDT-n). We will discuss several components of recently developed and implemented CC methodologies in NWChem. This includes: (1) Novel iterative/non-iterative methods accounting for the effect of triply excited configurations, (2) Massively parallel implementations of the CC theories based on the manifold of singly and doubly excited configurations. Several examples will illustrate how these approaches can be used in multiscale QM/MM framework.

Range-separated (screened) hybrid functionals provide a powerful strategy for incorporating nonlocal exact (Hartree-Fock-type) exchange into density functional theory. Existing implementations of range separation use a fixed, system-independent screening parameter. Here, we propose a novel method that uses a position-dependent screening. These locally range-separated (LRS) hybrids add substantial flexibility for describing diverse electronic structures and satisfy a high-density scaling constraint better than the fixed screening approximation does.

Joint work with P.U. Manohar.

A non-iterative N7 triples correction for the equation-of-motion coupled-cluster wave functions with single and double substitutions (EOM-CCSD) is presented. The correction is derived by second order perturbation treatment of the similarity-transformed CCSD Hamiltonian. The spin-conserving variant of the correction is identical to the triples correction of Piecuch and coworkers [Mol. Phys. 104, 2149 (2006)] derived within method-of-moments framework and is not size-intensive. The spin-flip variant of the correction is size-intensive. The performance of the correction is demonstrated by calculations of electronic excitation energies in methylene, nitrenium ion, cyclobutadiene, ortho-, meta-, and para- benzynes, 1,2,3-tridehydrobenzene, as well as C-C bond-breaking in ethane. In all cases except cyclobutadiene, the absolute values of the correction for energy differences were 0.1 eV or less. In cyclobutadiene, the absolute values of the correction were as large as 0.4 eV. In most cases, the corrections reduced the errors against the benchmark values by about a factor of 2 to 3, the absolute errors being less 0.04 eV.

The van der Waals density functional of Dion, Rydberg, Schroder, Langreth, and Lundqvist [Phys. Rev. Lett. 92, 246401 (2004)] will be reviewed, discussing implementations and applications by our group and others. New results relevalent for hydrogen storage in metal-organic framework (MOF) materials, as well for the intercalation of drug molecules in DNA will be presented.

I will overview some open mathematical questions related to the models and techniques of computational quantum chemistry. The talk is based upon a recent article coauthored with E. Cances and PL. Lions, and published in Nonlinearity, volume 21, T165-T176, 2008.

A number of exact relations are briefly discussed in terms of present developments and goals in DFT. In addition, conjectured relations are presented.

By means of rigorous thermodynamic limit arguments, we derive a new variational model providing exact embedding of local defects in insulating or semiconducting crystals. A natural way to obtain variational discretizations of this model is to expand the perturbation of the periodic density matrix generated by the defect in a basis of precomputed maximally localized Wannier functions of the host crystal. This approach can be used within any semi-empirical or Density Functional Theory framework. This is a joint work with Eric Cancès and Amélie Deleurence (Ecole Nationale des Ponts et Chaussées, France).

The quantal geometric phase [1-3] in a Born-Oppenheimer (adiabatic) electronic wavefunction is a net phase change for nuclear motion over a closed path. Effects of the quantal geometric phase have been observed in theoretical studies of the vibrational spectra of cyclic trinitrogen (N₃) molecule [4, 5]. Making a classical-quantum correspondence [6] relates the quantal geometric phase to a classical one for cyclic nuclear motion in N-body molecular dynamics. Each is described differential geometrically as the holonomy of a connection [7], physically in terms of the internal angular momentum, and with examples. The classical geometric phase [8] is a net angle of overall rotation in the center-of-mass frame. A net rotation of 20 degrees has been observed experimentally in a triatomic photodissociation and a net overall rotation of 42 degrees has been observed computationally in protein dynamics. The Hamiltonian operator in a generalized Born-Oppenheimer Schrodinger equation for the electronic wavefunction is related to a classical Hamiltonian for N-body molecular dynamics. Both the quantal and classical geometric phases arise due to non-zero internal angular momentum in N-body molecular dynamics.

References: [1] C. A. Mead and D. G. Truhlar, J. Chem. Phys. 70, 2284 (1979). [2] M. V. Berry, Proc. Royal Soc. London A 392, 45 (1984). [3] B. Simon, Phys. Rev. Lett. 51, 2167 (1983). [4] D. Babikov, B. K. Kendrick, P. Zhang, and K. Morokuma, J. Chem. Phys. 122, 044315 (2005). [5] D. Babikov, V. A. Mozhayskiy, and A. I. Krylov, J. Chem. Phys. 125, 084306 (2006). [6] F. J. Lin, Quantal and classical geometric phases, 2008. [7] J. E. Marsden, R. Montgomery, and T. Ratiu, Memoirs of the American Mathematical Society, Vol. 88, No. 436, American Mathematical Society, Providence, RI, 1990. [8] F. J. Lin, Discrete and Continuous Dynamical Systems, Supplement 2007, 655 (2007).

The quantal geometric phase [1-3] in a Born-Oppenheimer (adiabatic) electronic wavefunction is a net phase change for nuclear motion over a closed path. Effects of the quantal geometric phase have been observed in theoretical studies of the vibrational spectra of cyclic trinitrogen (N₃) molecule [4, 5]. Making a classical-quantum correspondence [6] relates the quantal geometric phase to a classical one for cyclic nuclear motion in N-body molecular dynamics. Each is described differential geometrically as the holonomy of a connection [7], physically in terms of the internal angular momentum, and with examples. The classical geometric phase [8] is a net angle of overall rotation in the center-of-mass frame. A net rotation of 20 degrees has been observed experimentally in a triatomic photodissociation and a net overall rotation of 42 degrees has been observed computationally in protein dynamics. The Hamiltonian operator in a generalized Born-Oppenheimer Schrodinger equation for the electronic wavefunction is related to a classical Hamiltonian for N-body molecular dynamics. Both the quantal and classical geometric phases arise due to non-zero internal angular momentum in N-body molecular dynamics.

References: [1] C. A. Mead and D. G. Truhlar, J. Chem. Phys. 70, 2284 (1979). [2] M. V. Berry, Proc. Royal Soc. London A 392, 45 (1984). [3] B. Simon, Phys. Rev. Lett. 51, 2167 (1983). [4] D. Babikov, B. K. Kendrick, P. Zhang, and K. Morokuma, J. Chem. Phys. 122, 044315 (2005). [5] D. Babikov, V. A. Mozhayskiy, and A. I. Krylov, J. Chem. Phys. 125, 084306 (2006). [6] F. J. Lin, Quantal and classical geometric phases, 2008. [7] J. E. Marsden, R. Montgomery, and T. Ratiu, Memoirs of the American Mathematical Society, Vol. 88, No. 436, American Mathematical Society, Providence, RI, 1990. [8] F. J. Lin, Discrete and Continuous Dynamical Systems, Supplement 2007, 655 (2007).

In this presentation I will give a review of the Cholesky Decomposition (CD) as it has been implemented in the MOLCAS program package. These examples will include conventional CD, as implemented for the HF, CASSCF, MP2, DFT, CASPT2 and CC methods, to the recent 1-center CD approximation. In addition, the aCD abd acCD techniques for the on-the-fly generation of RI auxiliary basis functions will be discussed. Analytic CD gradients will be introduced for CD-HF, CD-DFT(pure and hybrid), and CD-CASSCF. If time allows I will briefly discuss the use of CD technique for the fast evaluation of the exchange energy in CD-HF through the use of CD localized orbitals.

We present a novel multiscale modeling approach that can simulate multi-million atoms effectively via density functional theory. The method is based on the framework of the quasicontinuum (QC) approach with orbital-free density functional theory (OFDFT) as its sole energetics formulation. The local QC part is formulated by the Cauchy-Born hypothesis with OFDFT calculations for strain energy and stress. The nonlocal QC part is treated by an OFDFT-based embedding approach, which couples OFDFT nonlocal atoms to local region atoms. The method - QCDFT- is applied to a nanoindentation study of an Al thin film, and the results are compared to a conventional QC approach. The results suggest that QCDFT represents a new direction for the quantum simulation of materials at length scales that are relevant to experiments.

We study the general problem of mixing for ab-initio quantum-mechanical problems. Guided by general mathematical principles and the underlying physics, we propose a multisecant form of Broyden's second method for solving the self-consistent field equations of Kohn-Sham density functional theory. The algorithm is robust, requires relatively little fine-tuning and appears to outperform the current state of the art, converging for cases that defeat many other methods. We compare our technique to the conventional methods for problems ranging from simple to nearly pathological.

Optimization concepts will be reviewed with an eye on their proved or potential application in Electronic Structure Calculations and other Chemical Physics problems. We will discuss the role of trust-region schemes, line searches, linearly and nonlinearly constrained optimization, Inexact Restoration and SQP methods and the type of convergence theories that may be useful in order to explain the practical behavior of the methods. Emphasis will be given on general principles instead of algorithmic details.

In the quantum mechanical treatment of molecules we use the Born-Oppenheimer (adiabatic) approximation, in which the motion of nuclei and electrons is separated. In this approximation the coupling between different electronic states is neglected and nuclei move on a single electronic potential energy surface. Nevertheless, non-adiabatic processes where the coupling between different electronic states becomes large and important. These processes are facilitated by the close proximity of potential energy surfaces, and especially by the extreme case where the potential energy surfaces become degenerate forming conical intersections. Modeling non-adiabatic processes requires accurate calculation of electronic structure states and their coupling. Methods for calculating excited states, the non-adiabatic couplings and conical intersections will be discussed.

Standard approximations for the exchange-correlation functional have been found to give big errors for the linearity condition of fractional charges, leading to delocalization error, and the constancy condition of fractional spins, leading to static correlation error. These two conditions are now unified for states with both fractional charge and fractional spin: the exact energy functional is a plane, linear along the fractional charge coordinate and constant along the fractional spin coordinate with a line of discontinuity at the integer. This sheds light on the nature of the derivative discontinuity and illustrates the need for a discontinuous functional of the orbitals or density. This is key for the application of DFT to strongly correlated systems.

The quantum mechanical - molecular mechanical (QM/MM) potential energy functions are used in calculations of the potential of mean force (PMF) following the conventional molecular dynamics (MD) based procedure. Constant temperature MD simulations, in particular, allowing for rigid-body MD algorithms, are performed for the canonical (NVT) ensemble in conjunction with the Nose-Poincare thermostat. The umbrella sampling technique and the weighted histogram analysis method are applied for PMF calculations. Two versions of the QM/MM method, namely, the mechanically embedded cluster technique and the flexible effective fragment approach are considered for potential energy estimates. The QM subsystem can be described by various electronic structure approximations including the DFT and multiconfigurational CASSCF methods. Conventional force field parameters can be employed in the MM subsystem. The computer program utilizes the PC GAMESS (A. Granovsky) and TINKER (J. Ponder) molecular modeling packages. The first application considered the mechanisms of proton conduction in the gramicidin A ion channel. The chain of nine water molecules inside the channel constituted the QM part described by the B3LYP/6-31G* approximation. The peptide walls of the gramicidin channel and two clusters of 20 water molecules placed at both ends of the channel constituted the MM subsystem described by the AMBER force field parameters. For the energy consuming stage of the water file reorientation inside the channel we calculated the activation free energy barrier of 7.7 kcal/mol at 300K as compared to 6.5 kcal/mol in experimental studies.

The exponential localization of Wannier functions in two or three dimensions is proven for all insulators that display time-reversal symmetry, settling a long-standing conjecture. The proof make use of geometric techniques, which also imply that Chern insulators cannot display exponentially localized Wannier functions. Finally, a new algorithm to explicitly construct the exponentially localized Wannier functions is exhibited.

Principle Collaborator: Natarajan Sukumar (University of California, Davis)

Over the past few decades, the planewave (PW) pseudopotential method has established itself as the dominant method for large, accurate, density-functional calculations in condensed matter. However, due to its global Fourier basis, the PW method suffers from substantial inefficiencies in parallelization and applications involving highly localized states, such as those involving 1st-row or transition-metal atoms, or other atoms at extreme conditions. Modern real-space approaches, such as finite-difference (FD) and finite-element (FE) methods, can address these deficiencies without sacrificing rigorous, systematic improvability but have until now required much larger bases to attain the required accuracy. Here, we present a new real-space FE based method which employs modern partition-of-unity FE techniques to substantially reduce the number of basis functions required, by building known atomic physics into the Hilbert space basis, without sacrificing locality or systematic improvability. We discuss pseudopotential as well as all-electron applications. Initial results show order-of-magnitude improvements relative to current state-of-the-art PW and adaptive-mesh FE methods for systems involving localized states such as d- and f-electron metals and/or other atoms at extreme conditions.

This work performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.

Simulation of ballistic transport in nanodevices are usually computationally intensive. To improve the efficiency of this simulation, the subband decomposition method approximates the solution $psi(x_1,x_2)$ by $sum_{i=1}^{n_e} varphi_i(x_1) xi_i(x_2;x_1)$, where $xi_i(x_2;x_1)$, $1 leq i leq n_e$ are solutions to an eigenvalue problem in the $x_2$-direction and parameterized by $x_1$; $varphi_i(x_1)$ are solutions to a Schr"{o}dinger equation with open boundary conditions; and $(x_1,x_2)$ denotes a point in the simulation domain $Omega$. However, determination of $xi_i(cdot;x_1)$ based on, say, finite element method, at all grid points in the $x_1$ direction can still be expensive. Here we propose approximating $xi_i(cdot;x_1)$ using the reduced basis method. By exploiting {em a posteriori} error estimators and greedy sampling algorithm, we can construct very efficient reduced basis approximation space for $xi_i(cdot;x_1)$. The computational cost only grows marginally when mesh spacing in the $x_2$-direction decreases, compared to exponential increase for a finite element approximation.

The exchange-correlation energy of Kohn-Sham density functional theory can be written as the integral over all space of an exchange-correlation energy density, which is a function of various density-dependent ingredients. Adding ingredients can produce approximations that satisfy more exact constraints, or fit data better. In our vision, the ladder has five rungs: (1) the local spin density approximation (LSDA), which uses only the local electron spin densities as ingredients, (2) the generalized gradient approximation (GGA), which adds the density gradients, (3) the meta-GGA, which adds the orbital kinetic energy densities, (4) the hyper-GGA, which adds the fully-nonlocal exact exchange energy densities, and (5) the generalized random phase approximation, which adds the unoccupied Kohn-Sham orbitals. All rungs but the fourth have been constructed without empiricism. Some recent developments will be sketched: (a) The bifurcation of the second rung into standard GGA's and GGA's for solids. (b) Possible refinement of the meta-GGA by recovering the gradient expansion for exchange over a wide range of density gradients, as in the PBEsol GGA for solids. Since meta-GGA is not much more expensive than GGA, and is potentially much more accurate for systems near equilibrium, an improved meta-GGA could replace LSDA or GGA in applications. (c) A hyper-GGA that interpolates between meta-GGA exchange (in normal regions where the errors of meta-GGA exchange and correlation tend to cancel) and exact exchange (in abnormal regions where no such cancellation is possible). Extensive abnormal regions can occur in open subsystems connected by stretched bonds.

The first step in electronic-structure calculations is geometry optimization: finding an atomic configuration that minimizes the energy. A popular and successful model for the energy is the total-energy functional from density-functional theory. However, the evaluation of this functional at a given geometry is computationally expensive. Therefore, standard minimization techniques may be costly in practice. We propose a new approach for geometry optimization based on surrogate modeling. We describe our approach and present preliminary results.

Joint work with Filipp Furche.

We study the ability of TD-DFT to correctly predict the splitting of closely spaced singlet excited states in a system of two spatially separated chromophores. We find that functionals without at least a certain fraction of Hartree-Fock exchange kernel fare very poorly at this task because of a known problem these functionals have with the underestimation of charge-transfer excitation energies.

This presentation will address our current efforts to develop more accurate exchange-correlation forms for density functional theory. There are two leading themes in our current work: range separation and local weights. On the first theme, I will present a three-range hybrid functional and discuss the rationale for the success of screened functionals like HSE and LC-wPBE. On the second theme, the emphasis will be on new metrics for local hybridization and local range separation. Much of the focus will be on the seemingly dissimilar needs between solids and molecules, and on the computational challenge of including nonlocal (Hartree-Fock type) exchange efficiently in condensed systems.

Structures and binding in small water-benzene complexes (1-8 water molecules and 1-2 benzene molecules) are studied using the general effective fragment potential (EFP) method. The lowest energy conformers of the clusters were found using a Monte-Carlo technique. The EFP method accurately predicts structures and binding energies in the water-benzene complexes. Benzene is polarizable and consequently participates in hydrogen bond networking of water. Since the water-benzene interactions are only slightly weaker than water-water interactions, structures with different numbers of water-water, benzene-water, and benzene-benzene bonds often have very similar binding energies. This is a challenge for computational methods.

Conventional density functional approximations for the exchange-correlation energy fail to describe an important class of systems formed by the van der Waals (vdW) interaction, because they are unable to account for the long-range part of the vdW interaction, while they may describe the short-range part well. Here we first propose a density functional to simulate the coefficient C_6 of the leading term of the long-range part. Then we construct a nonempirical vdW-corrected meta-GGA functional by properly building the long-range part into a sophiscated meta-generalized gradient approximation (meta-GGA). Numerical tests on diverse atom pairs show that the proposed C_6 model is remarkably accurate. Applications of the vdW-corrected meta-GGA functional to rare-gas dimers show that the binding energy curves and bond lengths obtained with the vdW-corrected meta-GGA are well improved over those with the original meta-GGA, and agree fairly well with experiments.

Full configuration interaction (FCI) data are used to quantify the accuracy of approximate adiabatic connection (AC) forms in describing two challenging problems in density functional theory—the singlet ground state potential energy curve of H2 in a restricted formalism and the energies of the helium isoelectronic series, H− to Ne8+. For H2, an exponential-based form yields a potential energy curve that is virtually indistinguishable from the FCI curve, eliminating the unphysical barrier to dissociation observed previously with a [1,1]-Padé-based form and with the random phase approximation. For the helium isoelectronic series, the Padé-based form gives the best overall description, followed by the exponential form, with errors that are orders of magnitude smaller than those from a standard hybrid functional. Particular attention is paid to the limiting behavior of the AC forms with increasing bond distance in H2 and increasing atomic number in the isoelectronic series; several forms describe both limits correctly. The study illustrates the very high quality results that can be obtained using exchange-correlation functionals based on simple AC forms, when near-exact data are used to determine the parameters in the forms.

This lecture reports on work carried out in collaboration with Yan Zhao.

We have developed a suite of density functionals. All four functionals are accurate for noncovalent interactions and medium-range correlation energy. The functional with broadest capability, M06, is uniquely well suited for good performance on both transition-metal and main group-chemistry; it also gives good results for barrier heights. Another functional, M06-L has no Hartree-Fock exchange, which allows for very fast calculations on large systems, and it is especially good for transition-metal chemistry and NMR chemical shieldings. M08-2X and an earlier version, M06-2X, have the very best performance for main-group thermochemistry, barrier heights, and noncovalent interactions. M06-HF has no one-electron self-interaction error and is the best functional for charge transfer spectroscopy. A general characteristic of the whole suite is the optimized inclusion of kinetic energy density and higher separate accuracy of medium-range exchange and correlation contributions with less cancellation of errors than previous functionals [1-4]; for example, the functionals are compatible with a range of Hartree-Fock exchange and, although one or another of them may be more highly recommended for one or another property or application, all four are better on average than the very popular B3LYP functional. A few sample applications, including catalytic systems [5,6] and nanomaterials [7], will also be discussed. Recent work on lattice constants, band gaps, and an improved version of M06-2X will be summarized if time permits.

[1] "Design of Density Functionals by Combining the Method of Constraint Satisfaction with Parametrization for Thermochemistry, Thermochemical Kinetics, and Noncovalent Interactions," Zhao, Y. ; Schultz, N. E.; Truhlar, D. G.; J. Chem. Theory Comput. 2006, 2, 364-382. [2] "A New Local Density Functional for Main Group Thermochemistry, Transition Metal Bonding, Thermochemical Kinetics, and Noncovalent Interactions," Zhao, Y.; Truhlar, D. G. J. Chem. Phys. 2006, 125, 194101/1-18. [3] “The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals,” Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215-241. [4] "Density Functionals with Broad Applicability in Chemistry," Zhao, Y.; Truhlar, D. G. Acc. Chem. Res. 2008 41, 157-167. [5] “Attractive Noncovalent Interactions in Grubbs Second-Generation Ru Catalysts for Olefin Metathesis," Zhao, Y.; Truhlar, D. G. Org. Lett. 2007, 9, 1967-1970. [6] "Benchmark Data for Interactions in Zeolite Model Complexes and Their Use for Assessment and Validation of Electronic Structure Methods," Zhao, Y.; Truhlar, D. G. J. Phys. Chem. C 2008, 112, 6860-6868. [7] "Size-Selective Supramolecular Chemistry in a Hydrocarbon Nanoring," Zhao, Y.; Truhlar, D. G. J. Am. Chem. Soc.2007, 129, 8440-8442.

Constrained density functional theory methods (C-DFT) have been developed by several groups [1,2] to construct localized charge or spin states. At the same time, others have suggested ways of constructing variable-charge potential energy surfaces from certain valence-bond (VB) states [3,4]. In general, mainstream electronic structure codes do not generate the desired valence-bond states, and does not mix them in the manner needed to construct variable-charge potentials. C-DFT does provide estimates of the desired, individual valence-bond states and energies, but does not directly give the state mixing. While there are methods to overcome this limitation [5], here we explore an alternative. That alternative, notionally suggested in Ref. 3, is to deduce state-mixing from C-DFT with fractional charges. One assumes that the normal integer-charge VB states remain the same when describing a state composed of atoms with fractional charges. For instance, we use C-DFT to calculate the energies for the three charge distributions for a geometrically symmetric water molecule. Two of the constrained VB states are integer-charge distributions, one with all atoms being neutral and a second with the oxygen integer being anionic and one of the hydrogens cationic. The third VB state places half of a negative and positive charge on the oxygen and one hydrogen, respectively. We present our initial efforts to calculate the variational energy surface for a water molecule at several bond-lengths along the symmetric stretch vibrational mode. To capture the fractional charge behavior of the energies, the B3LYP functional is employed [5-7].

[1] Q Wu and T Van Voorhis, Phys. Rev. A 72, 024502 (2005). [2] J Behler, B Delley, K Reuter, and M Scheffler, Phys. Rev. B 75, 115409 (2007). [3] SM Valone and SR Atlas, J. Chem. Phys. 120, 7262 (2004). [4] SM Valone, J Li, and S Jindal, IJQC 108, 1452 (2008). [5] Q Wu and T Van Voorhis, J. Chem. Phys. 125, 164105 (2006). [6] AD Becke, J. Chem. Phys. 98, 5648 (1993). [7] C Lee, W Yang, and RG Parr, Phys. Rev. B 37, 785 (1988).

Nearly all common density functional approximations fail to adequately describe dispersion interactions responsible for binding in van der Waals complexes. One of the most promising new methods is the nonlocal van der Waals density functional (vdW-DF) of Ref. [1], which was derived from first principles, describes dispersion interactions in a seamless fashion, and yields the correct asymptotics. Recently we reported a self-consistent implementation of vdW-DF with Gaussian basis functions [2]. Our code includes analytic gradients of the energy with respect to nuclear displacements, enabling efficient geometry optimizations. vdW-DF tends to overbind molecular complexes, especially if used in combination with Hartree-Fock exchange. We propose a slightly simplified construction of the nonlocal vdW-DF correlation functional, for which we also derive a semilocal gradient correction. This correction reduces the overbinding tendency and improves the accuracy of vdW-DF. Adjusting an empirical parameter in the semilocal part, vdW-DF can be made compatible with different exchange approximations.

[1] M. Dion, H. Rydberg, E. Schroder, D.C. Langreth, and B.I. Lundqvist, Phys.Rev.Lett. 92, 246401 (2004).

[2] O.A. Vydrov, Q. Wu, and T. Van Voorhis, J.Chem.Phys. 129, 014106 (2008).

A new method for calculating effective atomic radii within the generalized Born (GB) model of implicit solvation is proposed, for use in computer simulations of bio-molecules. First, a new formulation for the GB radii is developed, in which smooth kernels are used to eliminate the divergence in volume integrals intrinsic in the model. Next, the Fast Fourier Transform (FFT) algorithm is applied to integrate smoothed functions, taking advantage of the rapid spectral decay provided by the smoothing. The total cost of the proposed algorithm scales as O(N^3logN+M) where M is the number of atoms comprised in a molecule, and N is the number of FFT grid points in one dimension, which depends only on the geometry of the molecule and the spectral decay of the smooth kernel but not on M. To validate our algorithm, numerical tests are performed for three solute models: one spherical object for which exact solutions exist and two protein molecules of differing size. The tests show that our algorithm is able to reach the accuracy of other existing GB implementations, while offering much lower computational cost.

I will present a direct constrained minimization (DCM) algorithm for solving the Kohn-Sham equations. The key ingredients of this algorithm involve projecting the Kohn-Sham total energy functional into a sequences of subspaces of small dimensions and seeking the minimizer of total energy functional within each subspace. The minimizer of a subspace energy functional not only provides a search direction along which the KS total energy functional decreases but also gives an optimal ``step-length" to move along this search direction. I will provide some numerical examples to demonstrate the efficiency and accuracy of this approach and compare it with the widely used method of self-consistent field (SCF) iteration. I will also discuss a few other numerical issues in algorithms designed to solve the Kohn-Sham equations.