Lennard-Jones Potentials for the Interaction of ... - ACS Publications

Subscriber access provided by READING UNIV

Article

Lennard-Jones Potentials for the Interaction of CO With 5-Membered Aromatic Heterocycles 2

Ángel Vidal-Vidal, Carlos Silva López, and Olalla Nieto Faza J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.7b09382 • Publication Date (Web): 08 Nov 2017 Downloaded from http://pubs.acs.org on November 9, 2017

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 A 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 43

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

Lennard-Jones potentials for the interaction of CO2 with 5-membered aromatic heterocycles ´ Angel Vidal-Vidal,† Carlos Silva López,† and Olalla Nieto Faza∗,‡ †Departamento de Qu´ımica Org´ anica. Campus Lagoas-Marcosende, 36310, Vigo, Spain ‡Departamento de Qu´ımica Org´ anica, Facultade de Ciencias, Universidade de Vigo. Campus As Lagoas, 32004, Ourense, Spain E-mail: [email protected] Phone: +34 9883 6888 Abstract We have used M06-2X/Def2-TZVPP to calculate a broad set of rigid interaction profiles between CO2 and 30 different aromatic heterocycles, based on pyrrole, furan and thiophene with ring positions subsituted with up to four nitrogens. For each system, several orientations of the fragments have been explored in order to both, find the preferred interaction mode and have information about other interaction modes that can contribute to the binding energy when CO2 is captured by complex systems. From these data, Lennard-Jones potentials have been obtained, which can be used for the parametrization of force fields that correctly describe the multipolar and dispersion interactions at play between these kind of fragments. These results are expected to contribute to the development of new force fields for the study of chemical systems for the capture and sequestration of CO2 and also directly for the design of such systems.

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

1

Introduction

Anthropogenic global warming is at the root of several environmental issues that will challenge our way of life, 1–3 among which melting glaciers, severe droughts, disruption of habitats and fast rising of sea levels due to ice melting at the Earth’s pole (which will lead to coastal flooding) are counted. Although it does not have the largest greenhouse effect per mole or weight, carbon dioxide (CO2 ) is one of the most important greenhouse gases present in the atmosphere 4,5 that are responsible of global warming, because of the amount of it that is liberated to the atmosphere by human activity. Thus, reducing CO2 emissions and capturing CO2 , either from the atmosphere or directly from the emission sources, is one of the key strategies in preventing climate change. For this, old technologies such as basic aqueous solutions of NaOH and KOH or amine scrubbing, 6–8 have been superseded and new strategies for the capture and sequestration of polluting gases are being developed. 9–12 Ionic liquids, 13,14 inorganic-organic interface composites, 15 porous inorganic membranes (PIMs), 16 polymers with light organic functional groups, 17–20 boron nitride nanotubes, 21 zeolitic imidazolate frameworks (ZIFs), 22–24 metal organic frameworks (MOFs), 10,25–28 and covalent organic frameworks (COFs) 29–31 are some of the chemical structures that are commonly used for capture and storage of CO2 . Molecular simulation plays a fundamental role in the obtention of new compounds that have improved properties over the existing ones such as high capture capacity, high recyclability, good thermal stability and easy desorption protocols, as well as a high selectivity in the adsorption of certain pollutants. With the help of computer-aided simulation it is possible to study, prior to the synthesis of these systems, static and dynamical properties such as affinity for different gases, capture capacity, isotherms, dynamics of adsorption, length scale, etc. and it is also possible to identify the degree of occupation of the structure. 32–35 Bulky solutions of soluble compounds in a specific solvent or extended solids with periodic boundary conditions are often used for this. 2


Page 2 of 43

Page 3 of 43

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


Many of the capture systems that it would desirable to tackle through molecular modeling are, unfortunately too large to be studied by means of quantum mechanics. Even though some of the electrons of the system can be ignored (ECPs) or some integrals approximated (as with semiempirical methods), a large number of atoms must still be considered, often resulting in a prohibitive computational cost (time/number of processors). To deal with this problem, molecular mechanics (MM) is invariably used to perform calculations on systems containing a large number of atoms, including systems in which the smallest spatial dimension is higher than nanometers. 36–38 There are numerous examples in the bibliography highlighting the use of molecular mechanics with force fields (FF) in systems of biochemical interest such as proteins and DNA, as well as in small and simple organic molecules. 39–46 The force fields used in MM contain a set of parameters developed from a relatively small set of molecules; these parameters can then be applied to a much wider range of systems based on structural resemblance using standard arithmetic and geometric combination rules, thus transferability is a key attribute of a force field. As a result, the force fields in use for simple organic or biological molecules are well developed and provide accurate information in most instances of routine use. However the need may arise to (re)parameterize a force field for complex systems with less conventional structures or in cases where intermolecular interactions beyond a van der Waals potential or simple pairwise electrostatics are important in their description. This would be the case with complex systems such as MOFs, ZIFs, COFs among others. The energy of the system in a force field is usually computed as a sum of terms, classified in two sets of bonded and non-bonded interactions. The first describe the energy that is required for distorting a molecule from a equilibrium reference in a specific fashion, the most common are the stretching of bonds, the bending of angles, and the torsional energy for the rotations around a bond. Non-bonded forces are through-space interactions that are usually modeled as a function of some inverse power of the distance. Two main groups can be defined among non-bonded interactions: electrostatic and van der Waals. Many

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

of the more common force fields such as DREIDING 47 and UFF 48 use a Lennard-Joneslike potential to deal with van der Waals interactions. Since 1,2 and 1,3 interactions are included in the bond and angle energy terms, non-bonded interactions are only computed between 1,4-neighbors. Electrostatic interactions are almost always included by the use of a Coulombic term on the energy. The problem with this is how to obtain a realistic set of charges that properly describes the system. For small molecules it is possible both, to fit charges to the electrostatic potentials calculated from high-quality ab-initio methods or to compute them using different approaches such as the Mulliken, Hirschfeld, CHELPG or other common charges, but a general method that produces accurate charges for larger molecules is desperately needed. Force fields have advanced much in the last 25 years with the goal of achieving two main objectives: reducing the computational cost and increasing the accuracy. There are some MM force fields such as DREIDING 47 that use a special hydrogen bond term to describe the interactions involving a hydrogen atom on an electronegative atom (N, O, F) associated with hydrogen bonds. In this case, it is possible to include either a CHARMM-like hydrogen bonding potential 47,49 or to perform an explicit inclusion of all van der Waals and electrostatic interactions and charges including the hydrogen to have a better description of these particular systems. In any case, and despite all the advances, it must be kept in mind that molecular mechanics, as a rule, can not provide properties that depend upon the electronic distribution of the molecule, as parameters are associated to atom types at the beginning of the simulation, and can not change in time. Efforts to improve accuracy have lead to approximations such as the central multipole expansion, which has been included within the energy terms to take account of higher order electrostatic interactions. 50,51 There is also research dealing with the description and parameterization of many-body interaction terms, 52–54 but work on a correct description of the dispersion energy component must still be addressed. And this is a problem in molecular mechanics descriptions of the intermolecular interac-

4


Page 4 of 43

Page 5 of 43

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


tions involving CO2 . Carbon dioxide is a linear non polar molecule (D∞h ) in which a carbon atom with sp hybridization is bonded to two sp2 oxygens. Although symmetry makes its dipole moment zero, CO2 exhibits a high charge separation, with negatively charged oxygens (Mulliken charges of -0.259 each) and a central carbon with positive charge (+0.518) (charges computed at the MP2-FC/cc-pV(T+d)Z level of theory). This charge separation results in a significant quadrupole moment 55–58 so that CO2 can establish high order electrostatic interactions such as quadrupole-dipole, quadrupole-quadrupole, quadrupole-octapole with a wide variety of molecules. One of the problems when modeling the absorption and adsorption of CO2 by complex structures is based in the existence of these multipolar interactions, which sometimes are the dominant energy contribution. As a result, if we want to correctly describe the guesthost chemistry of CO2 for its capture or storage, specific interaction potentials are needed that take into account electron correlation, dispersion and high order electrostatic moments. These potentials can be derived from either the newer density functionals or with postHartree-Fock methods. The interest in obtaining interaction potentials between CO2 and different molecules of chemical interest is not new, as the number of studies for the CO2 dimer between 1996 and 2000, for example, can attest. Some of the existing potential energy surface studies are calculated using the levels MBTP2/6-311+G(2df), using Symmetry-Adapted Perturbation Theory (SAPT) and also with MP2 with different size basis sets such as: [8s6p4d1f], [5s4p2d], [5s4p3d2f1g]. 52,59–62 Among these studies, it is worth noting the work of Kaluniga et al. 63 who performed an explicit correlation treatment of the PES of the CO2 dimer using high-quality ab-initio calculations. The analytical representation of the 4D PES was tested by comparison with a set of experimental data such as the second virial coefficient and the first spectral moment. Boulmène et al. 64 have mapped in 2D the interaction potentials between CO2 and different triazole isomers considering as variables the intermonomeric separation distance and also the angle. Tang et al. 65 using DFT calculations, have systematically studied the cation

5



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

effects on the interactions between a group of highly conjugated N-containing heterocycles (including neutral, anionic, and carbenoid species) and CO2 . There are also some works that treat the optimal geometry of some CO2 -heterocycle complexes but the investigations do not focus on the interaction potentials. 66–68 In this context, a systematic study of the interaction potentials between CO2 and heterocycles using methods that take into account dispersion and electronic correlation is still needed. Therefore, in this work the interaction potentials between CO2 and thirty 5-membered heterocycles (10 pyrrole derivatives, 10 furan derivatives and 10 tiophene derivatives) are studied. Given that the presence of substituted aromatic or heteroaromatic systems can lead to enhanced CO2 binding energies, 17,69 the heterocycles under study are substituted with between 0 and 4 nitrogen atoms. Throughout the study, multiple orientations of the CO2 with respect to the heterocycle are considered. Among them, the orientations parallel and perpendicular to the plane defined by the heterocycle, all the possible hydrogen interactions in two limit orientations as well as diverse profiles of interaction between the atoms of O and C of the CO2 with the main heteroatom of the heterocycle have been taken into account.

2

Computational methods

We will use DFT calculations to describe the interaction between CO2 and the heterocycles in Figure 1, as they have been shown to provide good quality results with a reasonable computational cost, including by its design a certain amount of dynamic electron correlation. Previous ab-initio and DFT benchmarks have pointed to the necessity of taking into account long range correlation effects for the study of non-covalent interactions in complexes between CO2 and other aromatic and non-aromatic molecules. 70 Dispersion terms have proved to be essential in the optimization of complexes where π- π and Van der Waals interactions are relevant. 67,70 They can be included in two ways; either incorporating empirical correction for the dispersion forces, as with the DFT-D3 meth-

6


Page 6 of 43

Page 7 of 43

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


+ 1

+ 1

+ 1 1 1 1 + WHWUD]ROH

+ 1 1 1 1 1 + SHQWD]ROH

+ 1

+ 1

+ S\UUROH

1 + LPLGD]ROH

+ S\UD]ROH

1,

1,

1 ,,

1,

1 ,,

1 ,,,

1 ,9

1,

1 ,,

1,

2

2

2

2

2

2

2

2

2

2

IXUDQ 2,

1 R[D]ROH 2,

6

6

WKLRSKHQH 6,

1 WKLD]ROH 6,

1

1

+ 1

LVR[D]ROH 2 ,, 1

6

LVRWKLD]ROH 6 ,,

+ 1 1 1 +

1 WULD]ROH

1 1 R[DGLD]ROH 2, 6 1 1 WKLDGLD]ROH 6,

+ 1 1 WULD]ROH

+

1

1 R[DGLD]ROH 2 ,, 1

6

1 WKLDGLD]ROH 6 ,,

+

1

1 1 WULD]ROH

+

WULD]ROH

1

1 1 R[DGLD]ROH 2 ,,, 6

1

R[DGLD]ROH 2 ,9 1

1 1 WKLDGLD]ROH 6 ,,,

1

6

1

WKLDGLD]ROH 6 ,9

1 1 1 R[DWULD]ROH 2, 6 1 1 1 WKLDWULD]ROH 6,

+ 1

1 1 + WHWUD]ROH 1

1

1

1 1 R[DWULD]ROH 2 ,, 6

1 1 WKLDWULD]ROH 6 ,,

1 1 1 1 WHWUDQLWURJHQ ,,, PRQRR[LGH 2, 6 1 1 1 1 WHWUDQLWURJHQ ,,, PRQRVXOILGH 6,

Figure 1: Representation of the chemical structure of the 30 systems that will be studied in this work. In the first row there are heterocycles derived from pyrrole, in the second group there are structures derived from furan and in the third one, heterocycles derived from tiophene. Apart from the IUPAC nomenclature, the notation that will be used throughout the discussion of systems and results is included. In this naming scheme, the first figure indicates the number of heteroatoms contained in the system, the letter written next corresponds to the main heteroatom of the heterocycle from which the systems are derived and finally, the Roman numerals are only used to differentiate the different possible positional isomers. Numeration with Roman numerals is random, so, it does not follow any chemical criteria. A color code is used to highlight the orientation that maximizes the binding energy for each heterocycle: blue=parallel, magenta=perpendicular, orange=IPOH, red=PPCH and green = IPCH. ods, 71 or including the description of non-covalent interactions in the parameterization of the functional, as with PWB6K, 72 M05-2X, 73 M06, 74 or M06-2X. 74 Density Functional Theory in the Kohn-Sham formulation as implemented in Gaussian 09 75 was used to locate minimum structures on the potential energy surface of CO2 and the 30 heterocycles in Figure 1, and to obtain electronic energies for the complexes resulting from their interaction, built as explained below. All calculations were done in the gas phase with the dispersion corrected hybrid-meta GGA functional M06-2X in combination with Ahlrich’s triple-ξ quality basis set, including polarization functions def2-TZVPP. 76,77 This basis set uses 2p1 d polarization for H, and 2d1 f for atoms B-Ne and Al-Ar. M06-2X functional was selected because of the good performance in the study of hydrogen bonds, H-π, π- π and electrostatic interactions both in neutral and charged dimeric systems. 66,78 Aiming for high accuracy, tight self-consistent field (SCF) criteria and an ultrafine pruned grid with 99 radial shells and 590 angular points per shell, were used for the numerical integration. Recently Vidal et al. showed the good performance

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

of this functional for the study of the structure, topology and spectroscopic characterization of CO2 -heterocycle complexes. 79 With the aim of testing the general performance of force fields (FF) in the description of interactions with CO2 , the DREIDING, 47 and UFF 48 force fields were selected, to perform calculations both using a conventional parameterization and with charge equilibration on all atoms. 80 We could not use AMBER, 81 a very popular force field, because it lacks parameters for the cumulenic carbon of CO2 , at least in the version implemented in g09. The charge equilibration method (QEq) proposed by Rappé et al. 80 was used to assign charges to all atoms in both DREIDING and UFF. QEq method uses as input data are experimental atomic ionization potentials, electron affinities, and also atomic radii. Within the system, an atomic chemical potential is constructed by using these quantities plus shielded electrostatic interactions between all charges. Requiring equal chemical potentials leads to equilibrium charges that depend upon geometry.

2.1

Potential parameterization procedure

To obtain a rich description of the potential energy surface, several limit orientations have been selected. If the maximum and minimum interaction that CO2 can have with the heterocycle for a certain orientation are characterized, any of the intermediate possible arrangements must have an energy comprised between those limits. To ensure that potential energy curves correspond with a specific geometry, atomic positions have been fixed and the energy was obtained by performing single point calculations, only changing the distance between fragments. If the system were allowed to optimize it would adopt the minimum energy configuration, thus precluding the sampling of modes of interaction that can be preferred when CO2 is captured inside materials due to confinement. As a result, these interactions would be not well described. The profiles explored, span internuclear distances from 2 ˚ A to 10 ˚ A in 0.5 ˚ A increments. Once the region in which attractive and repulsive interactions are balanced is known, more 8


Page 8 of 43

Page 9 of 43

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


points are calculated around the minimum in 0.25 ˚ A increments (and also in 0.1 ˚ A increments in some systems) in order to better characterize the region around the minimum of the curve. Energy versus distance curves are plotted in terms of relative energies, with the zero energy of the intermolecular interaction potentials corresponding to infinite separation of the molecules. As the distance between the two fragments increases, the interaction energy will tend asymptotically to zero. Once the graphical representation is available, we have proceeded to parameterize the interaction potentials following a classical Lennard-Jones equation where exponents 6 and 12 are fixed and the only adjustable parameters are the potential well depth, ǫ, and σ, a parameter that indicates the finite distance at which the potential between the particles is zero. This can be done in three different ways: 1) Since we have a continuous and derivable function throughout the closed interval in which the potential is studied [2, 10], and since it takes opposite signs at the extremes, according to Bolzano’s Theorem, there must be at least one value c ∈ (2 , 10) such that f (c) = 0. We obtain the root of the function, and subsequently the σ value, by defining increasingly smaller intervals at whose ends the function shows different signs. When the size of the interval has reached the required precision, we use its midpoint as the root value. 2) If the equation of the interaction potential is known, performing its derivative, and finding its roots, allows us to locate the geometric coordinates of the minimum. 3) With the aim of validating the previous approaches, a numerical adjustment of the computed points to a LJ equation where only σ and ǫ are adjustable parameters was also performed. The mean of the difference for each intermolecular distance r between a calculated LJ curve for the system and the points computed by ab-initio methods must be minimal. The nonlinear Generalized Reduced Gradient (GRG) algorithm proposed by Lasdon et al. was used to perform the optimization of the parameters. 82 Comparison of the three methods described provides differences in the third decimal of the adjustment for σ ′ and ǫ parameters, place in which the experimental error is located, for that

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

reason, all the methods can be used interchangeably to obtain the mentioned parameters.

2.2

Description of the heterocycles considered

In Figure 1 are represented the 30 heterocycles under study. We organize them in three groups, presenting them as derivatives of pyrrole, furan and tiophene, to which nitrogen atoms are added at different positions. The heteroatom defining the parent compound in each series is considered to be the main heteroatom. In order to avoid using the lengthy IUPAC names shown in Figure 1, a simple naming scheme is presented, based on one number, one letter and one roman numeral. The beginning number, which varies from 1 to 5, indicates the number of heteroatoms in the system, the letter indicates the main heteroatom of the parent system (N for pyrrole derivatives, O for furan derivatives and S for tiophene derivatives), and the roman number is arbitrarily chosen to differentiate between positional isomers of a given system.

3 3.1

Results Testing the M06-2X functional

The inclusion of dispersion is essential for the study of complexes stabilized by non-covalent interactions. There are different strategies for treating dispersion in the framework of DFT, and different functionals use different approaches. In order to test the validity of our functional of choice, we have performed a preliminary test with dispersion-including functionals belonging to different families. For this, we have used as a reference the high level computational results that S. Dalbouha et al. 83 have reported for the interaction of CO2 with the π-electronic density of imidazole. The authors provide the complexation energy of this complex computed at the CCSD(T)-F12/VTZ-F12 level of theory, correcting the result for the BSSE, a method that can be considered one of the most precise for the calculation of energies. 10


Page 10 of 43

Page 11 of 43

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


First, we have optimized the structure of the complex found in the literature at the M06-2X/def2-TZVPP level of theory and calculated the complexation energy corrected for the BSSE. On this same geometry, we have performed single point calculations to calculate binding energies with the following functionals: APFD, 84 B2PLYPD3, 85 B97D3, 86,87 and WB97XD, 88 which include Grimme’s D2 or D3 dispersion corrections. Since the geometry optimization at the M06-2X/def2-TZVPP may be a source of error on its own, we have also used the def2-tzvpp basis set with the APFD, B2PLYPD3, B97D3 and WB97XD density functionals to fully optimize the structure of the complex on which the BSSE corrected complexation energies have been calculated at the same level. Complexation energies obtained following the two approaches (single point calculations of the binding energy or optimization, followed by calculation of the binding energy ) and the percentage of error with respect to the reference can be found o Table 1. Table 1: Complexation energies corrected for the BSSE (in kJ/mol) with different density functionals and the def2-TZVPP basis set. The first two columns correspond to single point calculations on the M06-2X geometry. The last two columns list the binding energies calculated with each functional on the geometry optimized with that same functional. Percentual errors (%) are computed with respect to the CCSD(T)-F12/VTZ-F12+ BSSE energy (-14.355 kJ/mol) from the work of S. Dalbouha et al. 83 Density functional M062X APFD B2PLYPD3 B97D3 WB97XD

BE (kJ/mol) -14.675 -15.523 -11.339 -10.418 -10.000

Error (%) 2.23 8.13 -21.01 -27.42 -30.34

BE (kJ/mol) -15.941 12.385 -12.385 -11.673

Error (%) 11.05 -13.73 -13.73 -18.68

As can be seen in Table 1, the M06-2X functional, with an error of 2.23% provides a slightly over-stabilized complex with respect to the reference method. Of the other functionals tested on the M06-2X geometries, APFD also overstabilizes the complex, with an error of 8.13% with respect to the reference method, while B2PLYPD3, B97D3, and WB97XD yield more unstable complexes and larger percentual errors of -21.01%, -27.42% and -30.34% respectively. If a previous optimization with the functional is performed before binding energies are calculated, the results are similar. APFD overstabilizes the complex to a greater 11



extent (the error rises to 11.05%) and the other functionals underestimate binding energies with errors that increase as we go down the table. B2PLYPD3, -13.73%, B97D3, -13.73% and WB97XD, -18.68%. When using the last three functionals the results are better with the fully optimized geometries, but they are still far from providing acceptable results with respect to the reference method. From these results we conclude that M06-2X is the best functional to describe the behavior of these systems among those tested, so M06-2X/def2TZVPP is going to be the method of choice in all the calculations in this work.

3.2

Testing molecular mechanics force fields 0ár -2ár

E (kJ/mol)

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 43

-4ár -6ár -8ár -10ár -12ár 2ár

3ár

4ár

5ár

6ár

7ár

8ár

d (Å)

Figure 2: Representation of the interaction energy (in kJ/mol) vs. separation distance between the center of the 1H-imidazole ring and the C atom of CO2 calculated with different methods. The blue diamond shows the profile that is going to be used as a reference, obtained with M06-2X/Def2-TZVPP. The red diamond and the black circle correspond to the force fields UFF and DREIDING, respectively using charge equillibration (QEq) in all the atoms. The green triangle represents the results of UFF without QEq, and the pink square those obtained with the DREIDING force field in the same conditions. To test the general performance of force fields in describing the interaction between heterocycles and CO2 , we have chosen as a model a complex of CO2 with imidazole where the CO2 fragment lies parallel to the heterocyclic ring at different values of d, being d the distance between the center of the ring and the carbon atom in CO2 (See Figure 3 A). In this orientation, the formation of an electron donor-acceptor (EDA 89–92 ) complex is the most 12


Page 13 of 43

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


favorable, and it is also a common way of interacting found in macromolecular systems like MOFs, ZIFs and COFs. In Figure 2 we have represented the interaction potentials for the UFF and DREIDING force fields with and without the inclusion of charge equilibration, together with the curve computed with M06-2X/Def2-TZVPP, which has been used as a reference. In the curves calculated with the force-fields we appreciate a generalized decrease of the potential well depth (ǫ) and also a displacement of the curve to higher values of σ. With DREIDING but also with UFF, the inclusion of charge equilibration 80 substantially improves the results, but despite this, the difference is still remarkable, and improvement is needed if these potentials are to be used in the design of new materials for the capture of CO2 . parameterization of these intermolecular interaction potentials allows a quantitative comparison between them. The M06-2X reference values are 2.922 ˚ A, for σ and -10.459 kJ/mol for the potential well depth. UFF with charge balance, leads to a value of 2.991 ˚ A for σ, which represents a +2.36 % deviation. This deviation rises to +2.52 % if DREIDING is used. The deviations in the depth of the potential well are significantly larger, with variations of -16.28 % and -14.91 % for UFF and DREIDING, respectively. As already mentioned, the force fields without charge balancing provide poorer results (-47.38 % and -45.51 % deviations in ǫ for UFF and DREIDING, respectively and of 3.27 % and 3.41 % for σ). Table 2: Binding energies (BEs) for the systems 5N-I (1H-pentazole), 2N-I (1H-imidazole) 3O-IV (1.3.4-oxadiazole) and 4S-I (1,2,3,4-thiatriazole) computed using the supermolecule approach (kJ/mol). M06-2X/Def2-TZVPP (reference value) includes the Bernardi and Boys BSSE correction. 93 The BEs for the force fields are computed through single point energy calculations on the optimized M06-2X geometry. System M062X Dreiding Dreiding (Qeq=all) UFF UFF (Qeq=all)

5N-I -10.325 1.151 -0.370 1.362 -0.159

2N-I 3O-IV -14.675 -12.038 3.371 2.837 18.561 5.636 4.211 2.670 19.401 5.469

4S-I -10.404 2.114 -0.694 0.968 -1.840

For a correct description of the absorption or adsorption of CO2 by capture and storage 13



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

materials, accurate binding energies (BEs) are a must. In order to check the general accuracy of molecular mechanics force fields, we have calculated the BE of diverse CO2 -heterocycle complexes. Three model systems have been chosen: 5N-I-CO2 , 3O-IV-CO2 and 4S-I-CO2 . We optimized their geometries without constraints at the M06-2X/Def2-TZVPP level of theory and used the resultant structures to calculate BEs with the supermolecule approach. 94 The results, summarized in Table 2, display very large differences between the DFT and force field values, even with the beneficial use of charge equilibration. With FFs, favorable complexation energies are obtained only with 5N-I and 4S-I with the charge equilibration method, 80 and even there, the BE values do not exceed -0.7 kJ/mol (a 93 % difference with respect to the DFT reference at best). If balanced charges are not used, the complexes are predicted to be unstable, with larger errors over 110 %. Due to the particular electronic properties of oxygenated heterocycles, slightly different data are obtained for 3O-IV, where all the BEs computed with FFs are unfavorable to complexation. It should be noted that this is not an unexpected result. Whenever EDA studies have been reported about this kind of complexes, 66,68 it is found that the energy component due to dispersion is very important, an interaction that is poorly treated in most FFs. Table 2 also includes the CO2 -imidazole (2N-I) system that has been used in the previous subsection to calibrate the use of density functionals and the need of dispersion corrections. The results obtained are similar to those found for the 3O-IV system and in neither case is a stable complex obtained. Using the charge equilibration model makes the results worse for both systems, instead of improving them.. Again, the poor parameterization of these force fields for the description of the interactions in CO2 -heterocycle systems is revealed. There is a large difference between the interaction values computed with FFs through rigid profiles (See Figure 2) and those obtained from complexes optimized with DFT (See Table 2). This difference can be attributed to the minimization process making the two fragments adopt a configuration in which stabilizing interactions are maximized, and the fact that in this stabilization, the higher order electrostatic interactions and dispersion,

14


Page 14 of 43

Page 15 of 43

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


which are not well treated in FFs, are dominant. In the following sections we will explore with DFT the interaction potentials between CO2 and rigid heterocycles in different set orientations, in order to • get accurate Lennard-Jones parameters that can be used to build better FFs, • reproduce the experimental behavior of these structures and understand their interactions, • achieve predictive capacity about the preferred interactions, thus facilitating a` la carte design of new CO2 capture systems that use these motifs.

3.3

Rigid intermolecular interaction potential analysis and parameterization

In this section we will analyze the interaction between CO2 and the fundamental building blocks of the capturing systems that we want to develop. This study will be divided in three blocks: the first dealing with the interactions between CO2 and the heterocyclic backbone (carbon skeleton, π electron density...), the second with the interaction of CO2 with the main heteroatom on the heterocycle, and the third with the interaction between the two fragments through hydrogen bonding. In all three cases, different representative orientations have been considered to obtain the potential energy curves. With this we cover a wide range of interaction possibilities that should include the description of almost any situation that could be found in the absorption/adsorption events. 3.3.1

Interactions with π electron density: parallel and perpendicular orientations

Two spatial arrangements have been studied for the interaction of CO2 with the π electronic density of the heterocycles: parallel and perpendicular (A and B, respectively, in Figure 3). In the first (A), the CO2 lies parallel to the plane defined by the heterocycle, the 15


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

1N-I 2.910 -12.655 4.268 -0.445 1O-I 2.949 -9.498 4.212 -0.934 1S-I 3.033 -6.680 4.017 -2.256

System σP arallel (˚ A) ǫP arallel (kJ/mol) σP erpendicular (˚ A) ǫP erpendicular (kJ/mol)




16 2S-I 3.066 -5.558 3.995 -2.885

2O-I 2.959 -7.769 4.134 -1.572

2N-I 2.922 -10.459 4.157 -1.108

2S-II 3.053 -5.197 3.992 -3.225

2O-II 2.964 -6.847 4.074 -1.697

2N-II 2.920 -11.094 4.207 -0.905

3S-I 3.062 -4.577 3.970 -4.259

3O-I 2.978 -5.433 3.992 -2.648

3N-I 2.929 -9.228 4.066 -1.775

3S-II 3.056 -4.372 3.964 -4.350

3O-II 2.982 -4.904 3.983 -2.78

3N-II 2.942 -7.750 3.994 -2.400

3S-III 3.009 -5.012 3.973 -4.279

3O-III 2.941 -6.948 3.996 -2.402

3N-III 2.919 -10.120 4.102 -1.502

3S-IV 3.106 -3.393 3.969 -4.690

3O-IV 2.993 -3.805 3.973 -3.21

3N-IV 2.934 -8.586 4.022 -1.862

4S-I 3.128 -2.702 3.945 -5.790

4O-I 3.047 -2.289 3.949 -4.21

4N-I 2.954 -6.369 3.974 -3.201

4S-II 2.903 -3.983 3.942 -5.635

4O-II 2.989 -3.987 3.957 -3.89

4N-2 2.934 -8.070 3.984 -2.749

5S-I 3.062 -2.414 3.902 -7.366

5O-I 3.126 -1.374 3.907 -5.438

5N-I 2.958 -5.424 3.948 -4.132

Table 3: Compilation of σ (in ˚ A ) and ǫ (in kJ / mol) parameters for the interaction potential between CO2 and the heterocycles in Figure 1 in parallel and perpendicular orientations (see Figure 3) obtained by adjustment from the rigid energy profiles.

The Journal of Physical Chemistry Page 16 of 43

Page 17 of 43


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


8,0 6,0 4,0 2,0 0,0 E (kJ/mol)

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 18 of 43

-2,0 -4,0 -6,0 -8,0 -10,0 -12,0 -14,0 2,5

3ár

3,5

4ár

4,5

5ár

5,5

6ár

d (Å)

1S-I

2S-II

3S-IV

1O-I

2O-II

3O-IV

1N-I

2N-I

3N-II

Figure 4: Rigid interaction profiles for 9 selected systems in a parallel CO2 -heterocycle complex (orientation A in Figure 3). structural isomers, individual structures are considered. The trends in σ and ǫ values upon increasing N substitution of the heterocycles based in pyrrole, differ in the perpendicular orientation, reflecting a large difference in the interaction modes: σ decreases with the number of nitrogen atoms, in a logarithmic trend, with mean values of 4.268, 4.182, 4.046 , 3.979 and 3.948 ˚ A for the 1N, 2N, 3N , 4N and 5N systems, respectively. The depth of the potential well (ǫ) ranges from -0.445 kJ/mol for the 1N system and -4.132 kJ/mol for the 5N. In this orientation, interactions of the oxygen atom from CO2 with poly-N-substituted rings are more favorable than with low-N-substitution. However, the highest value that is reached in perpendicular orientation is a 23, 8% lower than the weakest ǫ value of parallel interactions. Heterocycles with oxygen display a similar behavior. The σ parameter for the parallel orientation follows the sequence 2.949 ˚ A, 2.962 ˚ A , 2.974 ˚ A, 3.018 ˚ A and 3.126 ˚ A for the 1O, 2O, 3O, 4O and 5O systems, respectively. The evolution of this trend however is not linear as in the case of systems with nitrogen, but it follows a polynomial third degree function. As the substitution with N atoms increases, there is a greater variation between the σ values

18


Page 19 of 43

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


of isomers. The variation of ǫ, is similar, where the value decreases from -9.498 kJ/mol of the system 1O to -1.374 kJ/mol of the fully substituted 5N system. The opposite trends are exhibited in the perpendicular orientation. In this case, σ values vary between 4.212 ˚ A in 1O and 3.907 ˚ A for 5N. The remaining individual values can be found in Table 3. Finally, when sulfur is the main heteroatom, there is not a clear relationship between σ and the number of nitrogen atoms on the heterocycle in the parallel orientation. It is noteworthy the small difference between the extreme σ values in the S systems (only 0.029 ˚ A), which contrasts with the wider ranges of N and O systems (0.048 and 0.177 ˚ A, respectively). The variations between the highest and lowest values in the perpendicular orientation are 0.320 ˚ A for pyrrole derivatives, 0.305 ˚ A for furan derivatives and 0.115 ˚ A (approximately one third less) for systems containing S as the main heteroatom. The σ variation for the perpendicular orientation follows a very clear decreasing linear trend ranging from 4.017 ˚ A for the 1S system to 3.902 ˚ A going through 3.994 ˚ A (2S), 3.969 ˚ A (3S) and 3.943 ˚ A (4S) values for the intermediate systems. If average values of σ and ǫ are used to perform a group comparison between rings containing N, O and S as the main heteroatom, we find that in the parallel orientation σ increases when going from N to S from 2.932 ˚ A up to 3.048 ˚ A. Taking into account that σ represents the finite distance at which the potential between the fragments is zero, this trend seems reasonable, with larger distances as the size of the heteroatom increases. In this orientation ǫ varies from -8.976 kJ/mol in pyrrole derivative heterocycles to -4.389 kJ/mol in tiophene derivatives (with an average of -5.285 kJ/mol for oxygen-based heterocycles). Contrary to what could be expected, this tendency is not linear in one of the oxygen based systems. This is a result of higher order electrostatic interactions having a large influence on ǫ and the binding energy of these systems, especially in the case of furan-derived heterocycles. The averaged evolution of σ for the perpendicular orientation follows a decreasing linear trend. ( 4.072, 4.018 and 3.967 ˚ A for the pyrrole, furan and tiophene families respectively). Although is not so marked, the trend in ǫ variation is not exactly linear and presents a slight

19




Page 20 of 43

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

1O-I 2.500 -5.952 2.416 -12.210 1S-I 3.043 -4.120 2.850 -7.861

System

σIP CH (kJ/mol) ǫIP CH (˚ A) σP P CH (kJ/mol) ǫP P CH (˚ A)

System


21

σIP CH (kJ/mol) ǫIP CH (˚ A) σP P CH (kJ/mol) A) ǫP P CH (˚

3.049 -3.380 2.984 -5.445

2S-I

2.538 -5.440 2.452 -8.918

2O-I

3.096 -3.335 2.892 -6.795

2S-II

2.583 -4.737 2.432 -11.035

2O-II

2.960 -3.967 3.115 -2.265

3S-I

2.466 -7.275 2.534 -4.005

3O-I

2.936 -4.483 3.094 -2.604

3S-II

2.456 -7.895 2.528 -2.265

3O-II

3.048 -2.235 3.112 -2.268

3S-III

2.510 -4.325 2.527 -4.525

3O-III

3.160 -2.333 2.931 -5.605

3S-IV

2.665 -3.363 2.449 -9.820

3O-IV

3.220 -1.371 3.012 -2.881

4S-I

2.482 -6.100 2.561 -2.695

4O-I

3.154 -1.175 3.200 -1.305

4S-II

2.532 -3.194 2.553 -3.160

4O-II

3.485 -0.089 3.385 -0.415

5S-I

2.559 -1.964 2.760 -1.665

5O-I

Table 4: Values of σ (in ˚ A ) and ǫ (in kJ / mol) parameters for the interaction potential between CO2 and the heterocycles in Figure 1 in IPCH and PPCH orientations (see Figure 5, obtained by adjustment from the rigid energy profiles.

Page 21 of 43 The Journal of Physical Chemistry


6,0 4,0 2,0 0,0

E (kJ/mol)

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 22 of 43

-2,0 -4,0 -6,0 -8,0 -10,0 -12,0 -14,0 2,0

3,0

4,0

5,0

6,0 d (Å)

7,0

8,0

1S-I

2S-II

3S-IV

1O-I

2O-II

3O-IV

1N-I

2N-I

3N-II

9,0

10,0

Figure 6: Rigid profiles for the interaction between CO2 and 9 selected heterocycles in a IPCH orientation (A in Figure 5). The average value of σ for the IPCH orientation in furan derivatives is 2.529 ˚ A, while it is substantially larger in tiophene derivatives: 3.115 ˚ A. Conversely, the potential well depth is 89.6% higher in the IPCH orientation with oxygen (-5.024 kJ/mol) than with sulfur (-2.649 kJ/mol). As for the PPCH orientation, σ values are slightly lower for both the oxygen system (2.521 ˚ A) and for those containing sulfur (3.058 ˚ A), however, ǫ values are substantially higher. The depth of the potential well increases a 16.7 % in the PPCH orientation with respect to the IPCH reaching a value of -6.030 kJ/mol. In tiophene-based systems, the increase of ǫ in PPCH with respect to IPCH is similar to the one obtained on oxygen-based heterocycles (16.4 %) since it increases from an average of -2.649 kJ/mol to -3.168 kJ/mol. The remaining values are shown on Table 4. Further analysis of Figure 6, reveals that the interaction curves corresponding to the pyrrole derivatives in the IPCH orientation have no minimum, which implies that these two units are never going to be attached in this spatial disposition. The rigid orientation assumed forces the system to interact through the H atom attached to the pyrrole nitrogen

22


Page 23 of 43



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

1N-I 3.152 -1.321 2.903 -4.950 1O-I 2.851 -2.205 No min 1S-I 3.317 -1.595 3.885 -0.252

System σP P OH (˚ A) ǫP P OH (kJ/mol) σIP OH (˚ A) ǫIP OH (kJ/mol)




24 2S-I 3.398 -1.402 3.765 -0.445

2O-I 2.838 -2.206 No min

2N-I 3.172 -1.065 2.860 -5.750

2S-II 3.440 -1.204 3.750 -0.571

2O-II 3.180 -1.654 3.456 -0.195

2N-II 3.088 -1.643 2.832 -5.628

3S-I 3.484 -0.915 3.470 -1.041

3O-I 3.184 -1.587 3.355 -0.316

3N-I 3.099 -1.462 2.758 -6.993

3S-II 3.466 -1.035 3.494 -0.872

3O-II 3.173 -1.652 3.406 -0.210

3N-II 3.178 -0.845 2.865 -5.432

3S-III 3.416 -1.208 3.466 -0.922

3O-III 3.077 -2.233 No min

3N-III 2.986 -2.263 2.740 -7.350

3S-IV 3.515 -0.766 3.354 -1.469

3O-IV 3.018 -1.055 2.842 -2.736

3N-IV 3.101 -1.464 2.776 -6.625

4S-I 3.590 -0.628 3.380 -1.574

4O-I 2.987 -1.163 2.995 -0.928

4N-I 3.052 -1.676 2.726 -8.980

4S-II 3.495 -0.775 3.360 -1.474

4O-II 3.177 -1.602 3.215 -0.473

4N-2 2.982 -2.156 2.712 -9.200

5S-I 3.718 -0.321 3.185 -2.220

5O-I 3.018 -0.988 2.902 -1.212

5N-I 2.965 -2.314 2.669 -11.520

˚ ) and ǫ (in kJ / mol) parameters for PPOH and IPOH orientations of CO2 and different Table 5: Compilation of σ (in A heterocyles obtained by adjustment from the computed energy profiles. Note that there is no minimum for three of the systems. As a result of the absence of an intersection with the x axis, these parameters can not be computed. The Journal of Physical Chemistry Page 24 of 43

Page 25 of 43

8,0 6,0 4,0 2,0 0,0 E (kJ/mol)

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,0 -4,0 -6,0 -8,0 -10,0 -12,0 -14,0 2,5

3,0

3,5

4,0

4,5

5,0

5,5

6,0

d (Å) 1N-I

2N-I

2N-II

3N-I

3N-II

3N-III

3N-IV

4N-I

4N-2

5N-I

Figure 8: Representative rigid energy profiles for the interaction of CO2 with several pyrrole derivatives in the IPOH orientation. This orientation involves the interaction between the main heteroatom of heterocycle and an oxygen atom of the CO2 molecule. The carbon dioxide molecule is in the plane defined by the heterocycle and the C=O· · ·heteroatom moiety is collinear. systems under study. As shown in Figure 7, the interaction between the fragments in the PPOH orientation is established between one of the oxygens on CO2 and the H atom in pyrrole heterocycles and with O and S in the furan and tiophene heterocycles, respectively. In the case of pyrrole derivatives, a hydrogen bond is expected to be the most stabilizing interaction, but this orientation fixates an O–H–N angle very unfavorable for a hydrogen bond, which are highly directional. In these cases, the IPOH orientation, where the O–H–N angle is linear, favoring the hydrogen bond, should provide a stronger interaction between the fragments. A selection of the interaction potentials for this two new limit orientations are shown in Figure 8, and the results for the parameterization of the interaction potentials for all the PPOH and IPOH orientations are displayed in Table 5. From the data in Table 5 we find in the PPOH orientation a mean value of σ in nitrogenbased systems of 3.071 ˚ A, which slightly decreases in the furan derivatives to 3.015 ˚ A, to increase again up to 3.493 ˚ A in tiophene derivatives. In this situation, the electronic 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

Page 26 of 43

properties of the oxygen atom cause a change in the expected behavior of the system again. The depth of the potential well (a proxy for the strength of the interaction between the fragments) is only larger in the pyrrole systems, than in their furan counterparts when the nitrogen substitution is very high, evidencing that hydrogen bond is not the predominant interaction in this arrangement. The trends in the IPOH orientation, however, differ from this picture, as expected for systems with a favorable hydrogen bond. σ, for example, increases steadily with the atomic number of the main heteroatom: N, O, S (as happened with the IPCH and PPCH orientations). However, the larger divergence is found in the values of ǫ. In this case the depth of the potential well is significantly larger in all pyrrolederivatives, increasing with the number of nitrogen atoms in the heterocycle up to a very large value of -11.520 kJ/mol for 5N-I. 3.3.3

Hydrogen bond study

;

+

+

;

+

+

;

+

1

1 +

+

;

+

1 1

+

;

+

+

+ ;,

;,

; ,, +

; 1

1

1 1

;

; ,, +

1 1

+

+

+ ;,

; 1 1 1

+

; ,,,

+

; 1

1 1

+ ; ,9

;,

; ,,

Figure 9: Generic representation of the heterocycles under study highlighting the candidates for hydrogen bonding. In the naming convention used: 1) the letter X symbolizes the main heteroatom: N, O or S, 2) color is used to mark the hydrogens according to their environment (symmetry makes same colored hydrogens equivalent): red, blue and pink. Note that the 5X systems are completely substituted with nitrogen atoms and therefore incapable of hydrogen bonding. We omit here hydrogen bonding between the oxygen in CO2 and the hydrogen on the main heteroatom of the pyrrole derivatives, as it has already been studied in the previous section (see PPOH and IPOH configurations).

26


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

1N-I (red) 3.355 -1.821

3N-I (red) 3.204 -3.191

1O-I (red) 3.255 -2.350

3O-I (red) 3.094 -4.084

1S-I (red) 3.285 -2.135

3S-I (red) 3.136 -3.705

System σHB (˚ A) ǫHB (kJ/mol)





27



3S-I (blue) 3.193 -3.083

1S-I (blue) 3.430 -1.502

3O-I (blue) 3.164 -3.398

1O-I (blue) 3.380 -1.759

3N-I (blue) 3.252 -2.486

1N-I (blue) 3.506 -1.175

3S-II (red) 3.315 -1.947

2S-I (red) 3.216 -2.697

3O-II (red) 3.054 -4.472

2O-I (red) 3.108 -3.505

3N-II (red) 3.152 -3.431

2N-I (red) 3.228 -2.697

3S-III (red) 3.249 -2.537

2S-I (blue) 3.278 -2.210

3O-III (red) 3.126 -3.696

2O-I (blue) 3.276 -2.281

3N-III (red) 3.246 -2.518

2N-I (blue) 3.292 -2.261

3S-IV (red) 3.148 -3.390

2S-I (pink) 3.226 -2.781

3O-IV (red) 3.107 -3.600

2O-I (pink) 3.202 -2.893

3N-IV (red) 3.210 -3.431

2N-I (pink) 3.355 -1.674

3S-IV (blue) 3.178 -2.965

2S-II (red) 3.325 -2.022

3O-IV (blue) 3.042 -4.689

2O-II (red) 3.216 -2.781

3N-IV (blue) 3.152 -3.447

2N-II (red) 3.316 -1.852

4S-I (red) 3.054 -4.755

2S-II (blue) 3.124 -3.593

4O-I (red) 2.983 -5.967

2O-II (blue) 3.274 -2.433

4N-I (red) 3.069 -4.576

2N-II (blue) 3.390 -1.699

4S-II (red) 3.902 -3.919

2S-II (pink) 3.221 -2.840

4O-II (red) 3.024 -5.110

2O-II (pink) 3.170 -3.231

4N-II (red) 3.110 -3.567

2N-II (pink) 3.270 -2.436

Table 6: Parameterization of the hydrogen bond interaction potentials for the forty-eight possible positions in the thirty heterocycles under study. The σ parameter (distance at which the interaction potential is zero) is shown in ˚ A , while the parameter that expresses the potential well depth (ǫ) is expressed in kJ/mol. The type of hydrogen atom being defined is expressed as a colour as illustrated in Figure 9.

Page 27 of 43 The Journal of Physical Chemistry


3,0

2,0

1,0 E (kJ/mol)

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 28 of 43

0,0

-1,0

-2,0

-3,0 2,5

3,0

3,5

4,0

4,5

5,0

5,5

6,0

d (Å) 1N-I-2/5

1N-I-3/4

2N-I-2

2N-I-4

2N-I-5

2N-II-3

2N-II-4

2N-II-5

Figure 10: Rigid profiles for the interaction between CO2 and several representative heterocycles through hydrogen bonding between C-H and the oxygen in CO2 . The CO2 fragment is considered to be coplanar and collinear with the C-H group with which it is interacting (similar to IPOH in Figure 7). Depending on the structure of CO2 capture materials, there may be positions in which a hydrogen bond (HB) may be established between one of the two electronegative oxygen atoms of carbon dioxide and one of the CH groups of the heterocycle, so a good treatment of such hydrogen bonds is key to a correct description of complex CO2 binding structures such as MOFs, ZIFs and COFs. However, this kind of hydrogen bonds in conventional force fields are either poorly or not parameterized. Specific parameterization of C–H–O hydrogen bond in this kind of heterocycles, where the environment can be strongly affected by the electronic properties of the heterocycle, depending on its main atom and the pattern of substitution with nitrogens, is thus needed. As a result, in this part of our work we have analyzed the interaction between the heterocycles and a CO2 fragment on the ring plane, through hydrogen bonding with the different hydrogens available. It should be noted that the hydrogen bond is highly directional, as shown in the PPOH and IPOH profiles, so a linear orientation has been selected for this study, where the CO2 is collinear with the C–H bond (the most favorable interaction). Figure 9 represents the sixteen different positions in which hydrogen bond energy profiles 28


Page 29 of 43

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


between CO2 and a given family of heterocycles have been calculated. In this figure, X represents the major heteroatom of the heterocycle ring (defining these families): NH, if derived from pyrrole, O for furan derivatives and S for tiophene derivatives. In addition, a color notation is used to differentiate the non-equivalent (based in symmetry considerations) hydrogens on the ring. There are sixteen possible hydrogen bond profiles per heterocycle, that is, forty-eight different profiles to be computed. Due to this large number, only eight have been selected for the graphical representation shown in Figure 10. All the profiles have been parameterized following the same procedure as in the previous sections. The results of this parameterization are shown on Table 6. From the results in Table 6 it can be appreciated that the σ and ǫ parameters associated to hydrogen bond between CO2 and C–H have quite similar values. The values for σ oscillate between 2.983 ˚ A and 3.902 ˚ A. The mean σ is 3.257 ˚ A for pyrrole derivatives, 3.155 ˚ A for those derived from furan and 3.268 ˚ A for those derived from tiophene. The potential well depths follow a similar trend, although there are larger variations between the heterocycles depending on their main atom. ǫ values range from -1.175 kJ/mol to -5.967 kJ/mol with a mean of -2.641 kJ/mol for pyrrole based heterocycles, -3.516 kJ/mol for those containing oxygen and -2.880 kJ/mol for those containing sulfur as the main heteroatom. The most energetic hydrogen bonds are those that occur within the sulfur heterocycles, followed by those in furan derivatives and finally with those containing the NH group as the main atom. There is a trend in the σ parameters for all the heterocycles. When a HB is established with CH groups placed two bonds away from the main heteroatom, the intermolecular interaction potential has a lower σ value. This behavior is followed in all systems where there is a C–H in this position except for the 2X-I one. In this particular case, only the pyrrole derivative follows the trend and has the minimum σ value for the pink hydrogen in Figure 9. For furan and tiophene derivatives, the lowest σ corresponds to the hydrogen atom shown in blue. All values for the parameterized HB curves are shown in Table 6. From the well depth values in Table 6, we show that this kind of interaction of CO2 with

29



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 hydrogens on the heterocyclic systems is not going to be the preferred binding mode of CO2 to these aromatic structures. However, these interactions are non-negligible, and in complex materials where more than one ring can be present, they can significantly impact the performance of the capturing system. Because of the interest in exploring approximation trajectories in dynamic studies, these potential energy profiles are also relevant in terms of the parameterization of FFs, as these interactions can be essential in determining which paths are preferred.

4

Conclusions

We have used DFT with a modern functional including dispersion and a triple ζ basis set (M06-2X/def2-TZVPP) to calculate a comprehensive set of rigid potential energy curves (energy vs. intermonomeric distance) for the interaction of CO2 with 30 five-membered aromatic heterocycles, based on pyrrole, furan and tiophene substituted with between 0 and 4 nitrogen atoms. We have fitted these curves to Lennard-Jones potentials in order to provide σ and ǫ values that can be used in the parameterization of force fields that correctly describe the multipolar and dispersion interactions at play between these kind of fragments. We have shown that two popular force fields (UFF and DREIDING) predict weaker CO2 imidazole interactions in the parallel orientation, and larger interaction distances (lower values for ǫ and larger σ values), although these results improve if charge equilibration is included. When binding Energies have been calculated for three representative, fully optimized (at the MO6-2X/def2-tzvpp level) heterocycle-CO2 complexes, force fields provide completely unreliable results (in some cases they even produce positive binding energies). Through rigid scans, we have explored a wide range of orientations between CO2 and the heterocycles, providing detailed information about the optimal interaction distance and the binding energy associated to each interaction mode. In each case, the effects on σ and ǫ of the main heteroatom in the cycle and the number of

30


Page 30 of 43

Page 31 of 43

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


nitrogens and their location on the ring have been discussed and, where possible, rationalized. This study has allowed us to identify the preferred orientation of CO2 with respect to each of the heterocycles studied, as a function of the depth of the potential well (some variations from these predictions can occur in some cases where the values are close, due to the fact that we are not performing full optimizations), as depicted with a color code in Figure 1. The fact that this preferred orientation is not always the same, showcases the rich variety of interactions at play, which can be exploited in technological applications. These results are expected to contribute to the development of new force fields for the study of chemical systems for the capture and sequestration of CO2 and also directly for the design of such systems. In the case of cage-like systems, for example, the best choice for number, type and spatial disposition of heterocyclic fragments in the system can be inferred from the profiles analyzed in order to maximize the binding energy.

Acknowledgement The authors thank the Centro de Supercomputación de Galicia (CESGA) for the generous allocation of computer time, and the Ministerio de Econom´ıa y Competitividad (MINECO, PCTQ2016-75023-C2-2-P) for funding. A. Vidal-Vidal is grateful to the Universidade de Vigo for a predoctoral fellowship.

Supporting Information Available SCF energies, Cartesian coordinates, number of imaginary frequencies for all computed structures. All this material is available as supporting information. This material is available free of charge via the Internet at http://pubs.acs.org/.

31



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

References (1) Khatib, H. IEA World Energy Outlook 2011–A comment. Energy Pol. 2012, 48, 737– 743, Special Section: Frontiers of Sustainability. (2) Wang, H.; Cao, D. Diffusion and Separation of H2 , CH4 , CO2 , and N2 in Diamond-Like Frameworks. J. Phys. Chem. C 2015, 119, 6324–6330. (3) Venna, S. R.; Carreon, M. A. Highly Permeable Zeolite Imidazolate Framework-8 Membranes for CO2 /CH4 Separation. J. Am. Chem. Soc. 2010, 132, 76–78. (4) Solomon, S.; Plattner, G.-K.; Knutti, R.; Friedlingstein, P. Irreversible Climate Change Due to Carbon Dioxide Emissions. Proc. Natl. Acad. Sci. 2009, 106, 1704–1709. (5) Monastersky, R. Global Carbon Dioxide Levels Near Worrisome Milestone. Nature 2013, 497, 13–14. (6) Tontiwachwuthikul, P.; Meisen, A.; Lim, C. CO2 Absorption by NaOH, Monoethanolamine and 2-Amino-2-Methyl-1-Propanol Solutions in a Packed Column. Chem. Eng. Sci. 1992, 47, 381 – 390. (7) Aronu, U. E.; Svendsen, H. F.; Hoff, K. A.; Juliussen, O. Solvent Selection for Carbon Dioxide Absorption. Energy Procedia 2009, 1, 1051 – 1057. (8) Luis, P. Use of Monoethanolamine (MEA) for CO2 Capture in a Global Scenario: Consequences and Alternatives. Desalination 2016, 380, 93 – 99. (9) Espinal, L.; Poster, D. L.; Wong-Ng, W.; Allen, A. J.; Green, M. L. Measurement, Standards, and Data Needs for CO2 Capture Materials: A Critical Review. Environ. Sci. Technol. 2013, 47, 11960–11975. (10) Sumida, K.; Rogow, D. L.; Mason, J. A.; McDonald, T. M.; Bloch, E. D.; Herm, Z. R.; Bae, T.-H.; Long, J. R. Carbon Dioxide Capture in MetalOrganic Frameworks. Chem. Rev. 2012, 112, 724–781. 32


Page 32 of 43

Page 33 of 43

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


(11) Sanz-Pérez, E. S.; Murdock, C. R.; Didas, S. A.; Jones, C. W. Direct Capture of CO2 from Ambient Air. Chem. Rev. 2016, 116, 11840–11876. (12) Sau, S. C.; Bhattacharjee, R.; Vardhanapu, P. K.; Vijaykumar, G.; Datta, A.; Mandal, S. K. Metal-Free Reduction of CO2 to Methoxyborane Under Ambient Conditions Through Borondiformate Formation. Angew. Chem. Int. Ed. 2016, 55, 15147–15151. (13) Niedermaier, I.; Bahlmann, M.; Papp, C.; Kolbeck, C.; Wei, W.; Krick Calderón, S.; Grabau, M.; Schulz, P. S.; Wasserscheid, P.; Steinr¨ uck, H.-P.; Maier, F. Carbon Dioxide Capture by an Amine Functionalized Ionic Liquid: Fundamental Differences of Surface and Bulk Behavior. J. Am. Chem. Soc. 2014, 136, 436–441. (14) Luo, X.; Guo, Y.; Ding, F.; Zhao, H.; Cui, G.; Li, H.; Wang, C. Significant Improvements in CO2 Capture by Pyridine-Containing Anion-Functionalized Ionic Liquids through Multiple-Site Cooperative Interactions. Angew. Chem. Int. Ed. 2014, 53, 7053–7057. (15) Prakash, M.; Mathivon, K.; Benoit, D. M.; Chambaud, G.; Hochlaf, M. Carbon Dioxide Interaction With Isolated Imidazole or Attached on Gold Clusters and Surface: Competition Between [Sigma] H-Bond and [Small Pi] Stacking Interaction. Phys. Chem. Chem. Phys. 2014, 16, 12503–12509. (16) Pera-Titus, M. Porous Inorganic Membranes for CO2 Capture: Present and Prospects. Chem. Rev. 2014, 114, 1413–1492. (17) Saleh, M.; Lee, H. M.; Kemp, K. C.; Kim, K. S. Highly Stable CO2 /N2 and CO2 /CH4 Selectivity in Hyper-Cross-Linked Heterocyclic Porous Polymers. ACS Appl. Mater. Interfaces 2014, 6, 7325–7333. (18) Kuwahara, Y.; Kang, D.-Y.; Copeland, J. R.; Brunelli, N. A.; Didas, S. A.; Bollini, P.; Sievers, C.; Kamegawa, T.; Yamashita, H.; Jones, C. W. Dramatic Enhancement of

33



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

CO2 Uptake by Poly(ethyleneimine) Using Zirconosilicate Supports. J. Am. Chem. Soc. 2012, 134, 10757–10760. (19) Manoranjan, N.; Won, D. H.; Kim, J.; Woo, S. I. Amide Linked Conjugated Porous Polymers for Effective CO2 Capture and Separation. J. CO2 Util. 2016, 16, 486 – 491. (20) Zulfiqar, S.; Awan, S.; Karadas, F.; Atilhan, M.; Yavuz, C. T.; Sarwar, M. I. Amidoxime Porous Polymers for CO2 Capture. RSC Adv. 2013, 3, 17203–17213. (21) Choi, H.; Park, Y. C.; Kim, Y.-H.; Lee, Y. S. Ambient Carbon Dioxide Capture by Boron-Rich Boron Nitride Nanotube. J. Am. Chem. Soc. 2011, 133, 2084–2087. (22) Huang, X.-C.; Lin, Y.-Y.; Zhang, J.-P.; Chen, X.-M. Ligand-Directed Strategy for Zeolite-Type MetalOrganic Frameworks: Zinc(II) Imidazolates with Unusual Zeolitic Topologies. Angew. Chem. Int. Ed. 2006, 45, 1557–1559. (23) Phan, A.; Doonan, C. J.; Uribe-Romo, F. J.; Knobler, C. B.; O’Keeffe, M.; Yaghi, O. M. Synthesis, Structure, and Carbon Dioxide Capture Properties of Zeolitic Imidazolate Frameworks. Acc. Chem. Res. 2010, 43, 58–67. (24) Banerjee, R.; Phan, A.; Wang, B.; Knobler, C.; Furukawa, H.; O’Keeffe, M.; Yaghi, O. M. High-Throughput Synthesis of Zeolitic Imidazolate Frameworks and Application to CO2 Capture. Science 2008, 319, 939–943. (25) Simmons, J. M.; Wu, H.; Zhou, W.; Yildirim, T. Carbon Capture in Metal-Organic Frameworks-a Comparative Study. Energy Environ. Sci. 2011, 4, 2177–2185. (26) Poloni, R.; Smit, B.; Neaton, J. B. Ligand-Assisted Enhancement of CO2 Capture in MetalOrganic Frameworks. J. Am. Chem. Soc. 2012, 134, 6714–6719. (27) Xue, D.-X.; Cairns, A. J.; Belmabkhout, Y.; Wojtas, L.; Liu, Y.; Alkordi, M. H.; Eddaoudi, M. Tunable Rare-Earth fcu-MOFs: A Platform for Systematic Enhancement of CO2 Adsorption Energetics and Uptake. J. Am. Chem. Soc. 2013, 135, 7660–7667. 34


Page 34 of 43

Page 35 of 43

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


(28) Babarao, R.; Jiang, J. Molecular Screening of MetalOrganic Frameworks for CO2 Storage. Langmuir 2008, 24, 6270–6278. (29) Babarao, R.; Custelcean, R.; Hay, B. P.; Jiang, D.-e. Computer-Aided Design of Interpenetrated Tetrahydrofuran-Functionalized 3D Covalent Organic Frameworks for CO2 Capture. Crys. Growth Des. 2012, 12, 5349–5356. (30) Buyukcakir, O.; Je, S. H.; Talapaneni, S. N.; Kim, D.; Coskun, A. Charged Covalent Triazine Frameworks for CO2 Capture and Conversion. ACS Appl. Mater. Interfaces 2017, 9, 7209–7216. (31) Pyles, D. A.; Crowe, J. W.; Baldwin, L. A.; McGrier, P. L. Synthesis of BenzobisoxazoleLinked Two-Dimensional Covalent Organic Frameworks and Their Carbon Dioxide Capture Properties. ACS Macro Lett. 2016, 5, 1055–1058. (32) Cabrales-Navarro, F. A.; Gómez-Ballesteros, J. L.; Balbuena, P. B. Molecular Dynamics Simulations of Metal-Organic Frameworks as Membranes for Gas Mixtures Separation. J. Membr. Sci. 2013, 428, 241 – 250. (33) Skoulidas, A. I. Molecular Dynamics Simulations of Gas Diffusion in MetalOrganic Frameworks: Argon in CuBTC. J. Am. Chem. Soc. 2004, 126, 1356–1357. (34) Amirjalayer, S.; Tafipolsky, M.; Schmid, R. Molecular Dynamics Simulation of Benzene Diffusion in MOF-5: Importance of Lattice Dynamics. Angew. Chem. Int. Ed. 2007, 46, 463–466. (35) Han, S. S.; Choi, S.-H.; van Duin, A. C. T. Molecular Dynamics Simulations of Stability of Metal-Organic Frameworks Against H2 O Using the ReaxFF Reactive Force Field. Chem. Commun. 2010, 46, 5713–5715. (36) Wongsinlatam, W.; Remsungnen, T. Molecular Dynamics Simulations of CO2 Molecules in ZIF-11 Using Refined AMBER Force Field. J. Chem. 2013, 35



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

(37) Chokbunpiam, T.; Chanajaree, R.; Remsungnen, T.; Saengsawang, O.; Fritzsche, S.; Chmelik, C.; Caro, J.; Janke, W.; Hannongbua, S. N2 in ZIF-8: Sorbate Induced Structural Changes and Self-Diffusion. Micropor. Mesopor. Mat. 2014, 187, 1 – 6. (38) Oxford, G. A. E.; Snurr, R. Q.; Broadbelt, L. J. Hybrid Quantum Mechanics/Molecular Mechanics Investigation of (salen)Mn for use in MetalOrganic Frameworks. Ind. Eng. Chem. Res. 2010, 49, 10965–10973. (39) Dickson, C. J.; Madej, B. D.; Skjevik, . A.; Betz, R. M.; Teigen, K.; Gould, I. R.; Walker, R. C. Lipid14: The Amber Lipid Force Field. J. Chem. Theor. Comput. 2014, 10, 865–879. (40) Steinbrecher, T.; Latzer, J.; Case, D. A. Revised AMBER Parameters for Bioorganic Phosphates. J. Chem. Theor. Comput. 2012, 8, 4405–4412. (41) Galindo-Murillo, R.; Robertson, J. C.; Zgarbová, M.; Sponer, J.; Otyepka, M.; Jurecka, P.; Cheatham, T. E. Assessing the Current State of Amber Force Field Modifications for DNA. J. Chem. Theor. Comput. 2016, 12, 4114–4127. (42) Perez, F.; Jaime, C.; Sanchez-Ruiz, X. MM2 Calculations on Cyclodextrins:Multimodal Inclusion Complexes. J. Org. Chem. 1995, 60, 3840–3845. (43) Deeth, R. J. Comprehensive Molecular Mechanics Model for Oxidized Type I Copper Proteins: Active Site Structures, Strain Energies, and Entatic Bulging. Inorg. Chem. 2007, 46, 4492–4503. (44) Morozov, A. V.; Misura, K. M. S.; Tsemekhman, K.; Baker, D. Comparison of Quantum Mechanics and Molecular Mechanics Dimerization Energy Landscapes for Pairs of RingContaining Amino Acids in Proteins. J. Phys. Chem. B 2004, 108, 8489–8496. (45) Lii, J.-H. Molecular Mechanics (MM4) Studies of Carboxylic Acids, Esters, and Lactones. J. Phys. Chem. A 2002, 106, 8667–8679. 36


Page 36 of 43

Page 37 of 43

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


(46) Sakakibara, K.; Naka, K.; Yamaguchi, Y.; Asami, M.; Chen, K.-H.; Lii, J.-H.; Allinger, N. L. Molecular Mechanics (MM4) Calculations on [3.3]- and [4.4]Orthoparacyclophanes. J. Phys. Chem. A 2004, 108, 3048–3055. (47) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. DREIDING: A Generic Force Field for Molecular Simulations. J. Phys. Chem. 1990, 94, 8897–8909. (48) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. UFF, a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. Soc. 1992, 114, 10024–10035. (49) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983, 4, 187–217. (50) Cieplak, P.; Dupradeau, F.-Y.; Duan, Y.; Wang, J. Polarization Effects in Molecular Mechanical Force Fields. J. Phys. Condens. Matter 2009, 21, 333102. (51) Dykstra, C. E. Electrostatic Interaction Potentials in Molecular Force Fields. Chem. Rev. 1993, 93, 2339–2353. (52) Tsuzuki, S.; Uchimaru, T.; Tanabe, K.; Kuwajima, S.; Tajima, N.; Hirano, T. Refinement of Nonbonding Interaction Parameters for Carbon Dioxide on the Basis of the Pair Potentials Obtained by MP2/6-311+G(2df)-Level ab Initio Molecular Orbital Calculations. J. Phys. Chem. 1996, 100, 4400–4407. (53) Borodin, O.; Smith, G. D. Development of Quantum Chemistry-Based Force Fields for Poly(ethylene Oxide) With Many-Body Polarization Interactions. J. Phys. Chem. B 2003, 107, 6801–6812. (54) Beran, G. J. O. Approximating Quantum Many-Body Intermolecular Interactions in

37



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

Molecular Clusters Using Classical Polarizable Force Fields. J. Chem. Phys. 2009, 130, 164115. (55) Kauffman, J. F. Quadrupolar Solvent Effects on Solvation and Reactivity of Solutes Dissolved in Supercritical CO2 . J. Phys. Chem. A 2001, 105, 3433–3442. (56) Reynolds, L.; Gardecki, J. A.; Frankland, S. J. V.; Horng, M. L.; Maroncelli, M. Dipole Solvation in Nondipolar Solvents: Experimental Studies of Reorganization Energies and Solvation Dynamics. J. Phys. Chem. 1996, 100, 10337–10354. (57) Graham, C.; Pierrus, J.; Raab, R. Measurement of the Electric Quadrupole Moments of CO2 , CO and N2 . Mol. Phys. 1989, 67, 939–955. (58) Varanasi, P.; Tejwani, G. D. T. On the Quadrupole Moment of CO2 . J. Chem. Phys. 1970, 53, 4404–4405. (59) Bock, S.; Bich, E.; Vogel, E. A New Intermolecular Potential Energy Surface for Carbon Dioxide From Ab Initio Calculations. J. Chem. Phys. 2000, 257, 147 – 156. (60) Bukowski, R.; Sadlej, J.; Jeziorski, B.; Jankowski, P.; Szalewicz, K.; Kucharski, S. A.; Williams, H. L.; Rice, B. M. Intermolecular Potential of Carbon Dioxide Dimer From Symmetry-Adapted Perturbation Theory. J. Chem. Phys. 1999, 110, 3785–3803. (61) Welker, M.; Steinebrunner, G.; Solca, J.; Huber, H. Ab Initio Calculation of the Intermolecular Potential Energy Surface of (CO2 )2 and First Applications in Simulations of Fluid CO2 . J. Chem. Phys. 1996, 213, 253 – 261. (62) Tsuzuki, S.; Uchimaru, T.; Mikami, M.; Tanabe, K. Intermolecular Interaction Potential of the Carbon Dioxide Dimer. J. Chem. Phys. 1998, 109, 2169–2175. (63) Kalugina, Y. N.; Buryak, I. A.; Ajili, Y.; Vigasin, A. A.; Jaidane, N. E.; Hochlaf, M. Explicit Correlation Treatment of the Potential Energy Surface of CO2 Dimer. J. Chem. Phys. 2014, 140, 234310. 38


Page 38 of 43

Page 39 of 43

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


(64) Boulmene, R.; Prakash, M.; Hochlaf, M. Microscopic Investigations of Site and Functional Selectivity of Triazole for CO2 Capture and Catalytic Applications. Phys. Chem. Chem. Phys. 2016, 18, 29709–29720. (65) Tang, H.; Lu, D.; Wu, C. Cation-Assisted Interactions Between N-Heterocycles and CO2 . Phys. Chem. Chem. Phys. 2015, 17, 15725–15731. (66) Lee, H. M.; Youn, I. S.; Saleh, M.; Lee, J. W.; Kim, K. S. Interactions of CO2 With Various Functional Molecules. Phys. Chem. Chem. Phys. 2015, 17, 10925–10933. (67) Hernández-Mar´ın, E.; Lemus-Santana, A. A. Theoretical Study of the Formation of Complexes Between CO2 and Nitrogen Heterocycles. J. Mex. Chem. Soc. 2015, 59, 36 – 42. (68) Chen, L.; Cao, F.; Sun, H. Ab Initio Study of the Interactions Between CO2 and Benzene, Pyridine, and Pyrrole. Int. J. Quantum Chem. 2013, 113, 2261–2266. (69) Biswas, S.; Vanpoucke, D. E. P.; Verstraelen, T.; Vandichel, M.; Couck, S.; Leus, K.; Liu, Y.-Y.; Waroquier, M.; Van Speybroeck, V.; Denayer, J. F. M.; Van Der Voort, P. New Functionalized MetalOrganic Frameworks MIL-47-X (X = Cl, Br, CH3 , CF3, OH, OCH3 ): Synthesis, Characterization, and CO2 Adsorption Properties. J. Phys. Chem. C 2013, 117, 22784–22796. (70) Vogiatzis, K. D.; Mavrandonakis, A.; Klopper, W.; Froudakis, G. E. Ab initio Study of the Interactions between CO2 and N-Containing Organic Heterocycles. ChemPhysChem 2009, 10, 374–383. (71) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104.

39



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

(72) Zhao, Y.; Truhlar, D. G. Design of Density Functionals That Are Broadly Accurate for Thermochemistry, Thermochemical Kinetics, and Nonbonded Interactions. J. Phys. Chem. A 2005, 109, 5656–5667. (73) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. Design of Density Functionals by Combining the Method of Constraint Satisfaction With Parametrization for Thermochemistry, Thermochemical Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2006, 2, 364–382. (74) Zhao, Y.; Truhlar, D. G. 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. Theor. Chem. Acc. 2008, 120, 215–241. (75) Frisch, M. J. et al. Gaussian 09 Revision D.01. Gaussian Inc. Wallingford CT 2009. (76) Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305. (77) Weigend, F. Accurate Coulomb-Fitting Basis Sets for H to Rn. Phys. Chem. Chem. Phys. 2006, 8, 1057–1065. (78) Kolaski, M.; Kumar, A.; Singh, N. J.; Kim, K. S. Differences in Structure, Energy, and Spectrum Between Neutral, Protonated, and Deprotonated Phenol Dimers: Comparison of Various Density Functionals With Ab Initio Theory. Phys. Chem. Chem. Phys. 2011, 13, 991–1001. (79) Vidal-Vidal, A.; Faza, O. N.; Silva López, C. CO2 Complexes with 5-Membered Heterocycles: Structure, Topology and Spectroscopic Characterization. J. Phys. Chem. A 2017,

40


Page 40 of 43

Page 41 of 43

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


(80) Rappe, A. K.; Goddard, W. A. Charge Equilibration for Molecular Dynamics Simulations. J. Phys. Chem. 1991, 95, 3358–3363. (81) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179–5197. (82) Lasdon, L. S.; Waren, A. D.; Jain, A.; Ratner, M. Design and Testing of a Generalized Reduced Gradient Code for Nonlinear Programming. ACM Trans. Math. Softw. 1978, 4, 34–50. (83) Dalbouha, S.; Prakash, M.; Timón, V.; Komiha, N.; Hochlaf, M.; Senent, M. L. Explicitly Correlated Interaction Potential Energy Profile of Imidazole + CO2 Complex. Theoretical Chemistry Accounts 2015, 134, 63. (84) Austin, A.; Petersson, G. A.; Frisch, M. J.; Dobek, F. J.; Scalmani, G.; Throssell, K. A Density Functional With Spherical Atom Dispersion Terms. Journal of Chemical Theory and Computation 2012, 8, 4989–5007. (85) Goerigk, L.; Grimme, S. Efficient and Accurate Double-Hybrid-Meta-GGA Density FunctionalsEvaluation With the Extended GMTKN30 Database for General Main Group Thermochemistry, Kinetics, and Noncovalent Interactions. Journal of Chemical Theory and Computation 2011, 7, 291–309. (86) Grimme, S. Semiempirical GGA-Type Density Functional Constructed With a LongRange Dispersion Correction. Journal of Computational Chemistry 2006, 27, 1787– 1799. (87) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. Journal of Computational Chemistry 2011, 32, 1456–1465. 41



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

(88) Chai, J.-D.; Head-Gordon, M. Long-Range Corrected Hybrid Density Functionals With Damped Atom-Atom Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615–6620. (89) Consan, K. A.; Smith, R. D. Observations on the Solubility of Surfactants and Related Molecules in Carbon Dioxide at 50C. J. Supercrit. Fluids 1990, 3, 51 – 65. (90) Raveendran, P.; Wallen, S. L. Cooperative CHO Hydrogen Bonding in CO2 Lewis Base Complexes: Implications for Solvation in Supercritical CO2 . J. Am. Chem. Soc. 2002, 124, 12590–12599. (91) Meredith, J. C.; Johnston, K. P.; Seminario, J. M.; Kazarian, S. G.; Eckert, C. A. Quantitative Equilibrium Constants between CO2 and Lewis Bases from FTIR Spectroscopy. J. Phys. Chem. 1996, 100, 10837–10848. (92) Kazarian, S. G.; Vincent, M. F.; Bright, F. V.; Liotta, C. L.; Eckert, C. A. Specific Intermolecular Interaction of Carbon Dioxide with Polymers. J. Am. Chem. Soc. 1996, 118, 1729–1736. (93) Boys, S.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures With Reduced Errors. Mol. Phys. 1970, 19, 553–566. (94) Morokuma, K.; Kitaura, K. In Chemical Applications of Atomic and Molecular Electrostatic Potentials; Politzer, P., Truhlar, D. G., Eds.; Springer Berlin Heidelberg: Berlin, Heidelberg, 1981; pp 215–242.

42


Page 42 of 43

Page 43 of 43

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


Graphical TOC Entry H N

S

N

N N

N

S

O

S

O

N

S

O

N N H N N

N H N

O N

N

N N S

H N

N S

S

H N

S

N

S

O

S

H N

N O

O N

N

N N

N N N N

N

N

N N

N

O

N

N N H N N

O N

N

N

N N

N

N

N

N N

N

N

N

N

N

H N

O N

N N

N N H N

N N H N N

N N

N

43


N

Lennard-Jones Potentials for the Interaction of ... - ACS Publications

Recommend Documents