Molecular Dynamics Simulation Study of Solvent ... - ACS Publications

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Molecular Dynamics Simulation Study of Solvent and State of Charge Effects on Solid Phase Structure and Counterion Binding in a Nitroxide-Radical Containing Polymer Energy Storage Material Travis W. Kemper, Thomas Gennett, and Ross E. Larsen J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.6b07118 • Publication Date (Web): 19 Oct 2016 Downloaded from http://pubs.acs.org on October 20, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Molecular Dynamics Simulation Study of Solvent and State of Charge Effects on Solid Phase Structure and Counterion Binding in a Nitroxide-Radical Containing Polymer Energy Storage Material Travis W. Kemper,† Thomas Gennett,‡ and Ross E. Larsen∗,¶ †Computational Science Center ‡Chemical and Nanoscience Center ¶Computational Science Center National Renewable Energy Laboratory 15013 Denver West Parkway, Golden CO 80401 E-mail: [email protected] Phone: (303) 275-4422 .

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract We performed molecular dynamics simulations to understand the effects of solvent swelling and state of charge (SOC) on the redox active, organic radical cathode material poly (2,2,6,6-tetramethylpiperidinyloxy methacrylate) (PTMA). We show that the polar solvent acetonitrile primarily solvates the nitroxide radical without disrupting the packing of the (2,2,6,6-tetramethylpiperidin-1-yl)oxyl (TEMPO) pendant groups of PTMA. We also simulated bulk PTMA in different SOC: 25%, 50%, 75% and 100%, by converting the appropriate number of TEMPO groups to the cation charge state and adding BF4− counterions to the simulation. At each SOC the packing of PTMA, the solvent and the counterions was examined. The binding of the anion to the nitroxide cation site was examined using the potential of mean force, and found to be on the order of tens of meV , with a binding energy that decreased with increasing SOC. In addition, we found that the cation state is stabilized by the presence of a nearby anion by more than 1 eV and the implications of this stabilization on charge transport are discussed. Finally, we describe the implications of our results for how the SOC of an organic electrode affects electron and anion charge transport during the charging and discharging processes.

1

Introduction

Improving energy storage remains a crucial hurdle in achieving many technological advances. Overcoming those obstacles will require an assortment of energy storage options tailored to an application. 1 Cost efficient high energy density storage materials are needed for electric vehicle and plug in hybrid electric vehicle applications. 2 For mobile devices, in addition to having a high energy density, high charge and discharge rates are needed. 3 In the area of high energy density storage materials Li-ion based technologies have dominated the market but as storage needs diversify, cost, fabrication and form factor become increasingly critical so exploring new materials is vital. One avenue that is being explored is using organic materials

2


Page 2 of 25

Page 3 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for energy storage, such as stable organic radical polymers to make organic radical batteries (ORB). 4,5 Organic materials have the advantage of being more environmentally friendly than traditional ceramic/inorganic based electronics. First, organic materials from renewable sources, such as bio feed stock, can replace rare or expensive elements used in the cathodes of Li-ion batteries. Also, with many strides being made in “green-synthesis,” by reducing the temperatures needed for processing and the elimination of metal catalysts, the environmental impacts of production can be reduced. Furthermore, upon disposal the organic materials in ORBs can be incinerated without releasing toxic chemicals into the atmosphere. 6 The ORB system has been well explored with promising results in regards to voltages, power and cyclability. Specifically, these systems allow for high charge and discharge rates that are ideal for peak power and mobile device applications. 7,8 Most of the ORB research has focused on the cathode material poly (2,2,6,6-tetramethylpiperidinyloxy methacrylate) (PTMA), which is composed of a poly-methylmethacrylate (PMMA) backbone modified with neutral nitroxide radical containing (2,2,6,6-tetramethylpiperidin-1-yl)oxyl (TEMPO) groups, as shown in Figure 1. PTMA has a theoretical capacity of 111 mAh/g, which is comparable to the LiCoO2 ’s working capacity of ∼ 130 mAh/g. 1,4 Unlike traditional crystalline inorganic cathode materials, the short range order in the PTMA coupled with its conductivity of 1 × 10−5 Scm−1 , leads to a unique and still poorly understood charge transport mechanisms. 9 In previous work we have established how the molecules are arranged in the solid phase, and the nature of the radical-radical interactions both theoretically and empirically. 10–13 This work builds on our previous efforts with the goal of understanding the relationship between how charged the organic cathode is, or its state of charge (SOC), and the mass transport and charge transport that occurs within the material. It is known that the charge transport in LiCoO2 is governed by solid state band transport, whereas charge transport in PTMA is governed by the hopping of a charge between pendant TEMPO groups. 14,15 Mass transport in LiCoO2 occurs in an inorganic lattice, whereas mass transport in PTMA

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

occurs in a disordered polymer that has been swelled by a solvent. The diffusion of charges from nitroxide site to nitroxide site has been quantified by a cation (hole) diffusion coefficient (Dh ) of 10−9 cm2 /s and the diffusion coefficient of anions (Danion ) through a polymer matrix has been quantified similarly at about 10−9 cm2 /s. 16 Yet charge/discharge rates greater than 50 C are easily achieved with ORBs. 7,8 This seeming mismatch between low diffusivities and rapid charge/discharge rates raises the question: what are the fundamental mechanisms associated with charge transfer and mass transport in PTMA? In this work we explore the structure and charge transport properties of PTMA as functions of solvent swelling and the SOC of a composite PTMA/electrolyte system. Here the SOC relates to how charged a battery would be if the PTMA material were the cathode of a battery. For example, if 50% of the unpaired electrons in PTMA were removed under an applied field leaving 50% of the radical sites in the cation charge state, the battery would be 50% charged and we would say the cathode has a SOC of 50%. We have performed molecular dynamics (MD) simulations of bulk PTMA in acetonitrile, with BF4− counter anions to keep the system charge neutral, for SOC of 0%, 25%, 50%, 75% and 100%. At each simulated SOC the charge neutral system of PTMA with radical and cation sites and incorporated solvent and anions was equilibrated using heating and cooling cycles to give quasi equilibrium snap shots of the molecular configuration of the cathode. In the simulations reported here we focus on equilibrium properties of the PTMA/electrolyte system with the goal of understanding such questions as: How does the cathode material swell in response to the introduction of anions, how do the anions couple to cation sites, and what is the effect of the cations on local electric potential? We then discuss the implications these equilibrium results have for charge transport properties relevant to ORB operation and to charged organic films more generally.

4


Page 4 of 25

Page 5 of 25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2

Approach and Numerical Methods

We used classical molecular dynamics (MD) simulations to simulate the structure of PTMA in the presence of acetonitrile solvent and with different SOC compensated by BF4− anions to maintain overall charge neutrality. Acetonitrile was chosen because of the known stability of the cation/radical of PTMA in the solvent, and the ease of conducting half-cell solution cell chemistry to compare empirical and theoretical results. Syndiotactic oligomers consisting of 24 repeat units were used to represent the PTMA polymer. 10 For each SOC a fraction of TEMPO radicals were changed to cations, evenly dispersed over the oligomer to approximate an even distribution of radical and cation sites within the charged material. Initial random configurations were generated by placing PTMA oligomers with the desired SOC at random points and under random rotations in a 140.0 ˚ A simulation cell, with periodic boundaries in all three cartesian directions. Inter-atomic distances were checked after each placement, and if the inter-atomic separation of an atom already in the system and a newly placed atom was less then the summed Van der Waals radii of those two atoms a new position and rotation of the PTMA oligomer was chosen. The Van der Waals radii for each atom type were taken from Bondi. 17 Once the PTMA oligomers were randomly placed in the periodic simulation cell, the BF4− anions were added corresponding to the SOC of the PTMA to create a charge-neutral simulation cell, using the same Van der Waals radius check. Finally, acetonitrile molecules comprising about 10 mass percent (see discussion in Section 3.1 below) were added on a grid, also with the same atomic Van der Waals radius check. Then each system was equilibrated using MD simulations using classical force-fields. The MD simulations were performed using the LAMMPS molecular dynamics code and a time step of 0.5 f s. 18 Temperature and pressure were controlled with a Nose-Hoover thermostat and a Parrinello-Rahman barostat respectively, both with time constants of 1.0 ps. The Coulombic interactions were handled using a particle-particle particle-mesh (pppm) solver with a cutoff of 12.0 ˚ A, which was also applied to the Lennard-Jones interactions. The MD simulations used the following force fields. The AMBER force-field was used for 6


Page 6 of 25

Page 7 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


PTMA repeat unit, with additional parameters specific to TEMPO taken from Stendardo, et al. 19 Atomic site charges were taken from density functional theory (DFT) calculations that we reported previously for neutral PTMA. 10 The atomic site charges on a TEMPO cation were taken from atomic charges found from an electrostatic-potential (ESP) fit of a methyl-terminated positively charged PTMA repeat unit, and are given in the supporting information (SI). Gaussian09 was used to used to calculated the ESP charges using density functional theory with the three parameter Becke-style hybrid exchange correlation functional B3LYP, and the 6-311++G(d,p) basis set. 20 For the neutral radical TEMPO, the virtual sites representing the lone pair electrons had a charge of 0.12 e, whereas for the cation case the virtual sites were retained but given a charge of 0.0 e. In addition, when we optimized the geometry of a methyl terminated PTMA repeat unit in the cation charge state, with the same functional and basis set described above, we found that the major change in geometry was a tipping of the N-O bond out of the plane defined by the nitrogen atom and its neighboring carbon atoms. Because the N-O group is surrounded by methyl groups, we decided that the geometry change would have a negligible effect on overall bulk structures so we took the cation and neutral radical geometries to be represented by the same force field except for aforementioned changes in the partial charges on each atom. Force field parameters for acetonitrile were taken from Jorgensen and Briggs 21 and for BF4− and the parameters were taken from Andrade, et al. 22 To be consistent with the charges used for PTMA, for both acetonitrile and BF4− the atomic charges we used were again taken from an ESP fit using density functional theory with the B3LYP exchange correlation functional and the 6-311++G(d,p) basis set, and these charges are also reported in the SI. To equilibrate the simulated systems, each simulation was heated and cooled with two successive cycles under the constant number pressure and temperature (NPT) thermodynamic ensemble, with a typical temperature history profile shown in Figure S2 of the SI. The initial heat/cool cycle increased the temperature of the system to 700 K and we then

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 25

varied the temperature between 400 K and 600 K before cooling the system down to 100 K. This initial heat cool cycle allowed the low-density system to relax to density greater than 1 g/cm3 . Once the initial condensed solid system was achieved the system was re-heated to 700 K then kept at 600 K for 200 ps before being cooled to 300 K at a rate of 50 K/ns at a pressure of 1 atm. Finally, the system was equilibrated at 300 K for 10 ns at a pressure of 1 atm, and data was collected during the last 5 ns of the 10 ns run. The heating cooling cycles performed do not necessarily guarantee that the system will be in a global free energy minimum, however, we believe that our system is in a local configuration that is at least nominally equilibrated. To verify that the system was in a stable configuration we examined the density of each simulation and the radius of gyration of each PTMA chain to be sure there was no drift during the data collection period. For example, the system density and PTMA chain radius of gyration at a SOC of 50% simulation are shown in the SI. The density of the simulation for 50% SOC, Figures S3 and S4, was found to maintain a constant value with fluctuation around the average of less than 0.1%. The standard deviation for the fluctuations of the radius of gyration for each PTMA chain was also found to be less than one percent of the ensemble average for radius at a SOC of 50%. Finally, we note that the configurations generated here represent a kinetically trapped state, as indicated in Figure S4, which shows that the probability distribution of chain radii for a single snapshot matches the distribution averaged over 5 ns, and which also shows that the variation in radius of gyration for any single chain is less than 2 ˚ A2 . In sum, this suggests that we are examining a kinetically trapped state representing some approximation to the equilibrium distribution of configurations. As in our previous work 10,13 we assume that charge transfer from a radical TEMPO group to a TEMPO cation can be described using the Marcus formalism of electron transfer, which gives the rate of charge transfer, from a site i to a site j as, ij kET

1/2 Vij2 π (λ + ∆G◦ )2 = exp − , h ¯ kB T λ 4λkB T 8


(1)

Page 9 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where T is temperature, Vij is the electronic-coupling matrix element between the sites i and j, λ is the reorganization energy, ∆Go is the free energy associated with the chargetransfer reaction and h ¯ and kB are Planck’s constant over 2π and Boltzmann’s constant, respectively. 23,24 We do not explicitly compute Marcus rates in this paper, but we seek insight into charge transfer in our system by studying ∆G◦ . We approximate the change in Gibbs free energy (∆Gij ) by the change in enthalpy (∆Hij ), corresponding to the difference in vertical ionization potential (vIP LEP ), of sites i and j, including both conformational effects and the local electrostatic potential on each site, as described previously. 13 In this work, we focus on the local electric potential Φi because we found 13 that the electrostatic potential had the largest effect on site energies. Φ was calculated based on the atomic charges of the surrounding medium,

Φi =

1 X qn 4πǫǫo n6=i |rn−i ~ |

(2)

where ǫo is the permittivity of free space and qn is the atomic charge of atom n. The dielectric constant, ǫ, was taken to be 3.4 based on capacitive-voltage measurements on PTMA. 25 This differs from the ǫ of 1 used in the MD simulation, because, the many-body effects that lead to a larger dielectric constant are built into the MD simulation by the molecular motions. In estimating a local potential based on a static snapshot from the surroundings, however, we do not have the motions or fluctuating polarization that lead to a dielectric constant greater than 1. Hence, to take this physics into account for a single snapshot we chose to use the experimentally determined dielectric constant of PTMA. As explained in previous work, 13 the set of atomic charges for the TEMPO group at site i were not included in the sum. For the neutral radical sites the cumulative atomic charge of each TEMPO group was summed into the connecting oxygen atom of the ester of the PMMA backbone and for the cation sites the cumulative atomic charge of TEMPO group (with the single cation charge of 1e subtracted) was summed into the connecting oxygen atom of the ester of the PMMA

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

backbone. These steps effectively charge terminated the ester, allowing the each TEMPO group to experience a correct local potential for its given charge state.

3

Results

3.1

Incorporation of Solvent into PTMA

In order to establish the solubility of acetonitrile in PTMA a simple saturation experiment was conducted, wherein we measured the mass increase of a PTMA film in a saturated acetonitrile atmosphere. 10 milligram samples of PTMA were placed in a platinum thermogravimetric analysis (TGA) pan and set on a stainless steel rod. A small amount of acetonitrile was placed in the bottom of a glass 20 ml petri-dish, the rod was placed in the liquid and a large no-lip beaker was placed over to create the acetonitrile atmosphere. The TGA pan containing the PTMA was removed periodically and TGA of the saturated PTMA was performed to confirm the weight of solvent absorbed, until a saturation mass was achieved, typically within 24hrs. This experiment was repeated three times and the mass increase was consistently 10%, within experimental error. Based on the saturation experiment we introduced 1444 acetonitrile molecules along with 92 oligomers to represent the solvent swollen polymer at 0% SOC. The resulting simulation contained 99928 atoms. After equilibration the density was found to be 1.07 g/cm3 , which was an increase from the 1.06 g/cm3 simulated previously for neat PTMA. 10 The simulated PTMA swelled with the addition of the solvent such that volume per oligomer of the PTMA/acetonitrile system increased by 10% over that of neat PTMA, as shown in Table 1. In order to understand the incorporation of solvent into bulk PTMA, we have calculated radial distribution functions (RDF) among different atomic sites for each equilibrated simulation. Figure 2(a) shows that the inter nitroxide nitrogen RDF for the neat and solvent swelled PTMA systems underwent no major inter-TEMPO structural changes. The inter oligomer and intra oligomer components are shown in Figure S6 in the SI, and also showed no 10


Page 10 of 25

Page 11 of 25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 25

Table 1: Properties of PTMA with and without incorporated electrolyte shown for different mass percentages of solvent and different states of charge. Shown are: density; volume per oligomer; change in volume per oligomer compared to neat PTMA; and the barrier for BF4− to escape a cation, as defined in the text. SOC (%) 0 0 25 50 75 100

solvent (%) 0 10 10 10 10 10

density g/cm3 1.06 1.07 1.09 1.14 1.17 1.21

volume/oligomer (nm3 ) 9.1 10.0 10.6 11.0 11.4 11.8

dV (%) – 10 17 21 25 30

anion barrier (meV) – – 95 75 64 60

significant change in peak locations. This indicates that the TEMPO groups do not change their packing from the three primary motifs identified in our previous work when PTMA is swelled by acetonitrile. 10 The unchanged packing implies that the electronic coupling, Vij , between nitroxide sites that enters Equation 1 will be largely unchanged by the presence of the solvent. So where does the solvent go in swelled PTMA? The RDF of the solvent nitrogen sites to the components of the PTMA backbone shown in Figure 2(d) reveals that the solvent primarily solvates the methyl carbons of the TEMPO group. In addition, we found that the PTMA backbones were pushed apart slightly when PTMA was swelled, as can be seen in the PTMA backbone-backbone RDF in Figure 2(b). This, combined with the lack of change in the inter-TEMPO RDF with solvent incorporation, shows that the acetonitrile molecules preferentially occupy regions between TEMPO rings that are created by the steric repulsion of the TEMPO methyl groups on each other, without altering the inter-TEMPO packing by interposing themselves between these pendant groups. We also found that acetonitrile interacts strongly with itself within the PTMA matrix. This is seen in the inter-solvent N -N and inter-solvent N -C(sp3 ) RDF’s shown in Figure 2(c). The inter-solvent N -C(sp3 ) RDF has a large peak at 4 ˚ A, while the solvent nitrogens have two peaks between 4 ˚ A and 6 ˚ A. This indicates that the nearest neighbor acetonitrile molecules have dipoles pointed in opposite directions. Similar looking RDFs have been

12


Page 13 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


reported for liquid acetonitrile, demonstrating that some short range order not dissimilar to that in liquids is obtained within PTMA. 26 We also computed, and show as insets, the cumulative number of neighbors inside a given distance from each site (hNj (< r)i), defined as Nf 1 XX hNj (< r)i = N (rij ) Ni Nf f r

Molecular Dynamics Simulation Study of Solvent ... - ACS Publications

Recommend Documents