Model interaction potentials for selenium from ab initio molecular simulations John C. Mauro NYS College of Ceramics at Alfred University, 2 Pine Street, Alfred, New York 14802, USA and Science and Technology Division, Corning Incorporated, SP-DV-02-08, Corning, New York 14831, USA Arun K. Varshneya NYS College of Ceramics at Alfred University, 2 Pine Street, Alfred, New York 14802, USA Received 10 November 2004; revised manuscript received 14 March 2005; published 21 June 2005 We develop model interaction potentials for elemental selenium based on ab initio molecular simulations and cluster expansion theory. Our potentials are used in classical Monte Carlo simulations to characterize the structure of selenium glass. DOI: 10.1103/PhysRevB.71.214105 PACS numbers: 61.43.Dq, 31.15.Ar, 31.50.Bc, 61.43.Fs I. INTRODUCTION Accurate descriptions of interatomic potentials are neces- sary to model and understand the structure and properties of materials. Ideally, these potentials are developed using a first-principles approach, which is most likely to accurately describe the interatomic forces. For larger atoms, however, this approach has proved difficult. As a result, researchers often use empirical or semiempircal models to describe the effective behavior of the atomic systems without necessarily capturing the true underlying physics. The most widely used potentials for selenium are based on a semiempirical model developed by Oligschleger and co-workers. 1 Following the approach of Stillinger et al. 2 in their model of sulfur, Oligschleger developed effective two- and three-body interaction potentials that reproduce known structures and energies of selenium clusters Se 2 -Se 8  from experimental data and density-functional calculations. The Oligschleger model of selenium has been used to determine selenium glass structure, 3 vibrational properties, 4 and quench behavior. 5–7 While the Oligschleger model has enabled much progress in the study of selenium, it is not based on funda- mental physics and hence is inherently limited in scope and applicability. This problem may be addressed using quantum- mechanical techniques. Shimojo et al. 8–10 investigated the structure and electronic properties of fluid selenium using density-functional theory, a quantum-mechanical technique by Hohenberg and Kohn 11 that relates the total electronic energy to the electron density. Zhang and Drabold 12,13 uti- lized additional localization techniques by Sankey and Niklewski 14 to study the impact of photon absorption on the structure of amorphous selenium. In addition, Shimizu et al. 15 and Nakamura and Ikawa 16 performed molecular-orbital calculations on amorphous selenium and found excellent agreement between their calculated structures and those de- termined by infrared/Raman spectrometry. However, none of this work at the quantum-mechanical level has attempted to bridge the gap to larger-scale molecular dynamics or Monte Carlo simulations. In this study, we use ab initio molecular simulations and a cluster expansion technique to derive effective interaction potentials for elemental selenium. We then use these poten- tials in classical Monte Carlo simulations to characterize the structure of selenium glass. II. SIMULATION DETAILS In our simulations, we employ second- and fourth-order Møller-Plesset perturbation theory 17 to calculate the energy of clusters of selenium atoms. We use the aug-cc-pVQZ ba- sis set of Dunning and co-workers, 18 where the acronym stands for “augmented correlation-consistent polarized va- lence quadruple .” For a single selenium atom, the aug-cc- pVQZ basis set consists of 93 basis functions composed of 343 primitive Gaussians. All ab initio simulations are per- formed using the GAUSSIAN 03 software. 19 In order to isolate the two- and three-body interaction potentials from the total system energy, we adopt the tech- nique of cluster expansion, 20 where the total potential of a system is the sum of the monomer energies and all combi- nations of higher-order interactions. Mathematically, we may write the total potential as V =  i=1 N V 1 +  i=1 N  ji N V 2 r ij  +  i=1 N  ji N  ki, j N V 3 r ij , r jk ,  ijk  + ¯ , 1 where V n refers to the n th -order interaction potential. In theory, the series of interactions terminates only with the V N term, where N is the total number of atoms in the system. However, since the magnitudes of the interactions typically decrease with increasing n and due to computational effi- ciency considerations, it is common to truncate the series after the second- or third-order terms. III. INTERACTION POTENTIALS In order to be used in higher-level classical simulations e.g., molecular dynamics and Monte Carlo, the ab initio potentials are fit to continuous functions that accurately re- produce the quantum data. While it is preferable to use fitting functions drawn from physical intuition, they do not neces- sarily provide the best fit for our particular systems. Here we PHYSICAL REVIEW B 71, 214105 2005 1098-0121/2005/7121/2141055/$23.00 ©2005 The American Physical Society 214105-1