EMSL: Eric Bylaska's Publications

Bylaska EJ, M Holst, and JH Weare. 2009. "Adaptive Finite Element Method for Solving the Exact Kohn-Sham Equation of Density Functional Theory." Journal of Chemical Theory and Computation 5(4):937-948. doi:10.1021/ct800350j Abstract Results of the application of an adaptive finite element (FE) based solution using the FETK library of M. Holst to Density Functional Theory (DFT) approximation to the electronic structure of atoms and molecules are reported. The severe problem associated with the rapid variation of the electronic wave functions in the near singular regions of the atomic centers is treated by implementing completely unstructured simplex meshes that resolve these features around atomic nuclei. This concentrates the computational work in the regions in which the shortest length scales are necessary and provides for low resolution in regions for which there is no electron density. The accuracy of the solutions significantly improved when adaptive mesh refinement was applied, and it was found that the essential difficulties of the Kohn-Sham eigenvalues equation were the result of the singular behavior of the atomic potentials. Even though the matrix representations of the discrete Hamiltonian operator in the adaptive finite element basis are always sparse with a linear complexity in the number of discretization points, the overall memory and computational requirements for the solver implemented were found to be quite high. The number of mesh vertices per atom as a function of the atomic number Z and the required accuracy e (in atomic units) was esitmated to be v (e;Z) = 122:37 * Z2:2346 /1:1173 , and the number of floating point operations per minimization step for a system of NA atoms was found to be 0(N3A*v(e,Z0) (e.g. Z=26, e=0.0015 au, and NA=100, the memory requirement and computational cost would be ~0.2 terabytes and ~25 petaflops). It was found that the high cost of the method could be reduced somewhat by using a geometric based refinement strategy to fix the error near the singularities.

Nichols PJ, N Govind, EJ Bylaska, and WA De Jong. 2009. "Gaussian Basis Set and Planewave Relativistic Spin-Orbit Methods in NWChem." Journal of Chemical Theory and Computation 5(3):491-499. doi:10.1021/ct8002892 Abstract Relativistic spin-orbit density functional theory (DFT) methods have been implemented in the molecular Gaussian DFT and pseudopotential plane-wave DFT modules of the NWChem electronic-structure program. The Gaussian basis set implementation is based upon the zeroth-order regular approximation (ZORA) while the planewave implementation uses spin-orbit pseudopotentials that are directly generated from the atomic Dirac-Kohn-Sham wavefunctions or atomic ZORA-Kohn-Sham wavefunctions. Compared to solving the full Dirac equation these methods are computationally efficient, but robust enough for a realistic description of relativistic effects such as spin-orbit splitting, molecular orbital hybridization, and core effects. Both methods have been applied to a variety of small molecules, including I$_{\text{2}}$, IF, HI, Br$_{\text{2}}$, Bi$_{\text{2}}$, AuH, and Au$_{\text{2}}$, using various exchange-correlation functionals. Our results are in good agreement with experiment and previously reported calculations.

Windus TL, EJ Bylaska, K Tsemekhman, J Andzelm, and N Govind. 2009. "Computational Nanoscience with NWChem." Journal of Computational and Theoretical Nanoscience 6(6 SP ISS):1297-1304. Abstract The NWChem software as been used to examine many nanoscale systems and their properties over the years. In this paper, an overiew of the general capabilities of NWChem is given as well as more specific details on the planewave and Gaussian based density functional codes usually used for nanoscale investigations. Examples are given of the scientific literature using NWChem, as well as two case studies: 1) Band gaps in oxides using exact-exchange based exchange-correlation functionals with the planewave DFT module, 2) Optical properties of chromophores using the Gaussian based DFT module.

Bylaska EJ, M Dupuis, and PG Tratnyek. 2008. "One-Electron-Transfer Reactions of Polychlorinated Ethylenes: Concerted and Stepwise Cleavages." Journal of Physical Chemistry A 112(16):3712-3721. Abstract Reaction barriers were calculated by using ab initio electronic structure methods for the reductive dechlorination of the polychlorinated ethylenes: C2CL4, C2Cl4, C2HCl3, trans-1,2-C2H2Cl2, cis-1,2-C2H2Cl2, 1,1-C2H2Cl2, and C2HCl3. Concerted and stepwise cleavages of R-Cl bonds were considered. Stepwise cleavages yielded lower activation barriers than concerted cleavages for the reduction of C2Cl4, C2HCl3, and trans-1,2-C2H2Cl2 via strong reducing agents. However, for typical ranges of reducing strength concerted cleavages were found to be favored. Both gas-phase and aqueous-phase calculations predicted C2Cl4 to have the lowest reaction barrier. Additionally, the reduction of C2HCl3 was predicted to have a significant amount of selectivity of cis-1,2-C2HCl2 over the corresponding reactions leading to the trans-1,2-C2HCl2, and 1,1-C2HCl2 radicals. These results illustrate how ab initio electronic structure methods, by providing experimentally inaccessible thermodynamics properties and activation energies, are able to sort out possible reactions mechanisms of reactions that have broad relevance in environmental chemistry.

Nichols P, EJ Bylaska, GK Schenter, and WA De Jong. 2008. "Equatorial and Apical Solvent Shells of the UO₂²⁺ Ion." Journal of Chemical Physics 128(12):124507. doi:10.1063/1.2884861 Abstract First principles molecular dynamics simulations of the hydration shells surrounding UO₂²⁺ ions are reported for temperatures near 300 K. Most of the simulations were done with 64 solvating water molecules (22 ps). Simulations with 122 water molecules (9 ps) were also carried out. The hydration structure predicted from the simulations was found to agree very well known results from X-ray data. The average U=O bond length was found to be 1.77Å . The first hydration shell contained five trigonally coordinated water molecules that were equatorially oriented about the O-U-O axis with the hydrogen atoms oriented away from the uranium atom. The five waters in the first shell were located at an average distance of 2.44Å (2.46Å - 122 water simulation). The second hydration shell was composed of distinct equatorial and apical regions resulting in a peak in the U-O radial distribution function at 4.59Å. The equatorial second shell contained 10 water molecules hydrogen-bonded to the five first shell molecules. Above and below the UO₂²⁺ ion, the water molecules were found to be significantly less structured. In these apical regions, water molecules were found to sporadically hydrogen bond to the oxygen atoms of the UO₂²⁺; oriented in such way as to have their protons pointed towards the cation. While the number of apical waters varied greatly, an average of 5-6 waters was found in this region. Many water transfers into and out of the equatorial and apical second solvation shells were observed to occur on a picosecond (ps) time scale via dissociative mechanisms. Beyond these shells, the bonding pattern substantially returned to the tetrahedral structure of bulk water.

Valiev M, EJ Bylaska, M Dupuis, and PG Tratnyek. 2008. "Combined Quantum Mechanical and Molecular Mechanics Studies of the Electron-Transfer Reactions Involving Carbon Tetrachloride in Solution." Journal of Physical Chemistry A 112(12):2713-2720. doi:10.1021/jp7104709 Abstract The reductive dechlorination of carbon tetrachloride, CC₄, was investigated using combined high level quantum mechanical and molecular mechanics (QM/MM) approach. The first electron transfer process was assumed to proceed by a concerted electron transfer-bond breaking mechanism, and reaction barriers for the first electron reduction were estimated by using the crossing point of the free energy profiles of CCl₃-Cl and CCl₃-Cl•- as a function of the CCl₃-Cl distance. The results of these calculations showed that the activation barriers for this reaction are reachable under a wide range of reduction potentials. In the gas-phase, the barrier to reduction varied from 0.8 kcal/mol for reducing agent with a -5 kcal/mol work function to 24.7 kcal/mol for a reducing agent with a 40 kcal/mol work function at the CCSD(T)/aug-cc-pVDZ level. In the aqueous phase, QM/MM calculations at the CCSD(T)/aug-cc-pVDZ level predicted that the barrier to reduction varied from 0.7 kcal/mol to 35.2 kcal/mol for -2.32 V and 0.93 V reduction potentials respectively. COSMO continuum solvation calculations were also performed for comparison. For strong reducing agents (EH < -1.5V) very little difference was seen between the QM/MM and COSMO activation barriers. For weak reducing agents (EH > 0V) the activation barriers differed by as much as 6 kcal/mol between the QM/MM and COSMO calculations. These results demonstrate that ab initio electronic structure methods coupled with explicit molecular mechanics representation of the aqueous environment offer an efficient and accurate way to calculate the free energy reaction barriers for dissociative electron transfer reactions of organochlorine compounds to identify the potentially important environmental degradation processes.

Bylaska EJ, M Valiev, JR Rustad, and JH Weare. 2007. "Structure and Dynamics of the Hydration Shells of the Al3+ Ion ." Journal of Chemical Physics 126(10):Art.no.104505. Abstract First principles simulations of the hydration shells surrounding Al3+ ions are reported for temperatures near 300oC. The predicted six waters in the octahedral first hydration shell were found to be trigonally coordinated via hydrogen-bonds to 12 second shell waters in agreement with the putative structure used to analyze the X-ray data, but in disagreement with results reported from conventional molecular dynamics using two- and three-body potentials. Bond lengths and angles of the water molecules in the first and second hydration shell and the average radii of these shells also agreed very well with the results of the X-ray analysis. Water transfers into and out of the 2nd solvation shell were observed to occur on a picosecond (ps) time scale via a dissociative mechanism. Beyond the second shell the bonding pattern substantially returned to the tetrahedral structure of bulk water. Most of the simulations were done with 64 solvating waters (20 ps). Limited simulations with 128 waters (5 ps) were also carried out. Results agreed as to the general structure of the solvation region and were essentially the same for the first and second shell. However, there were differences in hydrogen-bonding and Al-O radial distribution function in the region just beyond the second shell. At the end of the second shell a nearly zero minimum in the Al-O radial distribution was found for the 128 water system. This minimum is less pronounced minimum was found for the 64 water system, which may indicate that sizes larger than 64 may be required to reliably predict behavior in this region,

Du J, LR Corrales, K Tsemekhman, and EJ Bylaska. 2007. "Electron, hole and exciton self-trapping in germanium doped silica glass from DFT calculations with self-interactions correction." Nuclear Instruments and Methods in Physics Research. Section B, Beam Interactions with Materials and Atoms 255(1 (SP ISS)):188-194. Abstract We performed density functional theory (DFT) calculations of electron, hole and exciton self-trapping in germanium doped silica glass to understand the refractive index change in these glasses induced by UV irradiation. The local structure relaxation and excess electron density distribution upon trapping of the above species were calculated. The results show that both trapped exciton and electron are highly localized on germanium ion and, to some extent, on its oxygen neighbors. Exciton self-trapping is found to lead to the formation of Ge E’ center and non-bridging hole center. Electron trapping changes the GeO4 tetrahedron structure into trigonal bi-pyramid with the majority of the excess electron density located along the equatorial line. Self-trapped hole is localized on bridging oxygen ions that are not coordinated to germanium atoms and leads to elongation of the Si-O bonds and change of the Si-O-Si bond angles. We did comparative study of standard DFT vs. DFT with a hybrid PBE0 exchange and correlation functional. The results show that the two methods give qualitatively similar relaxed structure and charge distribution for the electron and exciton trapping in germanium doped silica glass; however, only using the PBE0 functional reproduces the hole self-trapping. This research is supported by the Divisions of Chemical Science, Office of Basic Energy Sciences, US Department of Energy. This research was performed in part using the Molecular Science Computing Facility in the William R. Wiley Environmental Molecular Sciences Laboratory (EMSL) at the Pacific Northwest National Laboratory (PNNL). The EMSL is funded by DOE’s Office of Biological and Environmental Research. The pacific Northwest National Laboratory is operated by Battelle for the US Department of Energy.

Gao F, J Du, EJ Bylaska, M Posselt, and WJ Weber. 2007. "Ab Initio Atomic Simulations of Antisite Pair Recovery in Cubic Silicon Carbide." Applied Physics Letters 90(22):Art. No. 221915. doi:10.1063/1.2743751 Abstract The thermal stability of an antisite pair in 3C-SiC is studied using ab initio molecular dynamics within the framework of density functional theory. The lifetime of the antisite pair configuration is calculated for temperatures between 1800 and 2250 K, and the effective activation energy for antisite pair recombination is determined to be 2.52 eV. The recombination energy path and static energy barrier are also calculated using the nudged elastic band method, along with the dimer method to accurately locate the transition states. The consistency of the results suggests that the antisite pair cannot be correlated with the DI photoluminescence center, as proposed by previously theoretical interpretations. An extended exchange mechanism is found for the antisite pair recombination, and this may be a dominant mechanism for antisite pair recombination and diffusion of impurities in compound semiconductors.

Bickmore BR, KM Rosso, CJ Tadanier, EJ Bylaska, and D Doud. 2006. "Bond-valence methods for pKa prediction. II. Bond-valence, electrostatic, molecular geometry, and solvation effects ." Geochimica et Cosmochimica Acta 70(16):4057-4071. doi:10.1016/j.gca.2006.06.006 Abstract In a previous contribution, we outlined a method for predicting (hydr)oxy-acid and oxide surface acidity constants based on three main factors: bond valence, Me–O bond ionicity, and molecular shape. Here electrostatics calculations and ab initio molecular dynamics simulations are used to qualitatively show that Me–O bond ionicity controls the extent to which the electrostatic work of proton removal departs from ideality, bond valence controls the extent of solvation of individual functional groups, and bond valence and molecular shape controls local dielectric response. These results are consistent with our model of acidity, but completely at odds with other methods of predicting acidity constants for use in multisite complexation models. In particular, our ab initio molecular dynamics simulations of solvated monomers clearly indicate that hydrogen bonding between (hydr)oxo-groups and water molecules adjusts to obey the valence sum rule, rather than maintaining a fixed valence based on the coordination of the oxygen atom as predicted by the standard MUSIC model.