Prediction of Structure and Properties of Boron-Based Covalent

Article pubs.acs.org/JPCC

Prediction of Structure and Properties of Boron-Based Covalent Organic Frameworks by a First-Principles Derived Force Field Saeed Amirjalayer,† Randall Q. Snurr,‡ and Rochus Schmid*,† †

Lehrstuhl für Anorganische Chemie 2, Organometallics and Materials Chemistry, Ruhr-Universität Bochum, Universitätsstr. 150, D-44780 Bochum, Germany ‡ Department of Chemical and Biological Engineering, Northwestern University, Evanston, Illinois 60208, United States S Supporting Information *

ABSTRACT: We present a method to theoretically predict structures in arbitrary network topologies for all currently known boron based covalent organic frameworks (COFs). This is particularly useful because these materials are accessible experimentally only as polycrystalline powders. The method is based on a new fully flexible molecular mechanics force field. The consistent parameter set is derived by a genetic algorithm optimization approach from first-principles reference computed data. To achieve high accuracy, the convergence with respect to the level of theory is carefully controlled for this reference. The force field is used to investigate the relative stability of the two high symmetry topologies ctn and bor. Interestingly, for all systems, the ctn topology is found to be energetically more stable. This preference is observed experimentally, too, with the single exception of COF-108, which forms the bor topology. This exception can thus be attributed to the different synthesis conditions, demonstrating that other topologies might be accessible in principle for all COFs. The force field is further used to compute first benchmark surface areas for ideal systems, thermal expansion coefficients, elastic constants, and CO2 adsorption isotherms for all systems in both topologies, which are experimentally unavailable. Our force field opens the way for theoretical structure prediction and prescreening of properties for these fascinating materials.

1. INTRODUCTION Covalent organic frameworks (COFs) have recently been introduced as a new class of functional porous materials.1−3 Similar to the so-called metal organic frameworks (MOFs),4−6 they are crystalline polymeric networks, assembled from molecular building blocks and maintain their crystallinity and porosity after removal of the guest molecules. In contrast to MOFs, COFs are not coordination polymers but are built from light main group elements connected by covalent bonds. Because of their high porosity, they show the lowest densities known for crystalline materials,2 which makes them particularly interesting for application in gas adsorption of, e.g., hydrogen, methane, or carbon dioxide.7−9 Unfortunately, COFs are obtained as polycrystalline powders, which makes structural characterization difficult, and one has to rely on structural assumptions as a starting point for the Rietveld refinement procedure. For the boron based 3D networked COFs with trigonal and tetrahedral building blocks, the formation of high symmetry (edge transitive) 3,4-networks can be assumed. Two such topologies are possible, in principle,10,11 namely bor and ctn, as shown in Figure 1 as augmented nets. The topology of all known 3D-COFs was deduced by a comparison of the measured and simulated PXRD pattern for these topologies. Surprisingly, all but one are found to form the ctn topology, and for unknown reasons, only COF-108 crystallizes in the bor form, however, under substantially different synthesis con© 2012 American Chemical Society

Figure 1. High symmetry network topologies (edge transitive) for 3,4heteroconnected nets with tetrahedral and trigonal building units.

ditions.1 Now the question arises, whether this is an intrinsic consequence of the building blocks of COF-108 or solely due to the solvothermal conditions as is usually the case for zeolites.12,13 Note that different phases have been crystallized also for MOFs, depending on the nature of the solvent.14,15 Thus, it would be desirable to assess the relative thermodynamic stability of the different topological isomers in order to Received: November 23, 2011 Revised: January 13, 2012 Published: January 23, 2012 4921

dx.doi.org/10.1021/jp211280m | J. Phys. Chem. C 2012, 116, 4921−4929

The Journal of Physical Chemistry C

Article

periodic systems at the necessary theoretical level is not straightforward.16 However, molecular mechanics force fields can be used for strain energy analysis since the topological isomers reside on the same molecular mechanics potential energy surface (PES). The condition is that an accurate parametrization must be available. For purely organic systems, accurate parametrizations like, for example, the MM3 and MM4 force fields30−32 or MMFF94,33 are available, whereas for the hybrid MOFs, only generic rule based force fields like UFF34 can be used. Here, the accuracy is neither known nor can it systematically be improved. To this end, we recently introduced a systematic parametrization strategy using first-principles reference data and a genetic algorithm (GA) optimization strategy with a novel objective function to derive parameters for different metal organic frameworks.35,36 The present contribution has two aims. By using a high-level parametrized force field, which is explicitly developed to treat this new class of porous materials, the structures of these threedimensional networks are studied on an atomistic level. Thus, an understanding of the intrinsic factors determining the formation of a specific topology can be gained, which is only possible with the computation of reliable steric energies also for large systems in a sequential multiscale fashion. Further, on the basis of these results, mechanical, thermal, and adsorption properties of the COFs are predicted. To our knowledge, this kind of bottom-up investigation has not been done for the class of covalent organic frameworks. Thus, for the first time, these materials can be analyzed independent from experimental data, so that an unrestricted prescreen is possible, which is crucial in order to develop rational materials with designated properties.

answer this question. In addition, in the case of COF-105 and -108, the solvent molecules could not be fully removed and currently neither experimental nor theoretical adsorption data for 3D-COFs in bor topology are available. Theoretical investigations are playing a key role in understanding the host−guest interactions of porous functional materials in general.16−18 Up to now, theoretical studies of COFs dealt mainly with adsorption of light gases like hydrogen or methane.8,19−21 Also, the possibility to improve hydrogen uptake by a doping of COFs with lithium atoms has been studied this way.22−25 However, in the majority of cases, the network structure was taken from experimental data, which limits the predictive capabilities of the theoretical methods. It is obvious that also the topology of the network influences the host−guest interactions in a drastic way, and, thus its tuning is an important but ambitious goal. The basic assumption is that the linking of secondary building units (SBUs) (Figure 2) with

Figure 2. SBUs considered in this work together with the augmented vertex figures. Note that the central atom in the tetrahedral SBU, shown on the left, is either carbon or silicon.

2. COMPUTATIONAL DETAILS All quantum mechanical calculations of the nonperiodic reference systems were performed using the TURBOMOLE suite of programs.37 Geometry optimization and normal-mode analysis were performed by density functional theory (DFT) using the hybrid functional B3LYP.38 A fine integration grid (m5) and a tight SCF criterion (10−10 a.u.) were used to accurately predict also the low vibrational modes. The basis sets cc-pVXZ (X = D, T, Q) were employed for all elements.39 Single point calculations on the double-hybrid level B2PLYP, as introduced by Grimme,40 and on the MP2 level were performed, in order to study the impact of dispersion effects. All MP2-type calculations (also in B2PLYP) were performed using the resolution of identity (RI-MP2) approximation.41 The efficient RICC2 module42 implemented in the TURBOMOLE package was used together with the standard auxiliary basis sets.43 Atomic charges were obtained by fitting the electrostatic potential (ESP) using the Merz−Kollman sampling schema44 as implemented in the TURBOMOLE code. Second derivatives were computed analytically with the AOFORCE module. The classical molecular mechanics calculations were performed using a locally modified version of the TINKER code.45 The force field energy expression used in this work is a slightly extended version of MM3,30 which was used previously in ref 35. In contrast to standard MM3, electrostatic interactions were modeled by atomic point charges. We used a cutoff of 12.0 Å for the van der Waals (vdW), and the shortrange electrostatic interactions and the long-range electrostatic interactions were computed by a smooth particle mesh Ewald summation (SPME). For details on the energy expression and a complete listing of the parameters, see the Supporting Information. Standard periodic boundary conditions were

a well-defined conformation and shape leads to a specific topology, and in case of shape-conserving variations, this results in a family of isoreticular materials.26 However, because of the conformational flexibility of the linkers, in principle, a wide variety of networks could be formed. In the case of zeolites, the rotational freedom of the T−O bond (T = Si, Al) allows for a plethora of structures, with only a small subset yet synthesized.27 In order to understand and control the formation of certain framework topologies, it remains to be clarified which factors govern the solvothermal growth process. Because of the slow solvothermal growth process of MOFs and COFs, it is plausible that thermodynamic factors dominate, and the topology with the lowest free energy (including the solvent) will be formed. Two major factors determining this free energy can be discriminated. First, there are extrinsic effects, governed by the synthesis conditions like, for example, the nature of the solvent, the temperature, or the concentrations of the building blocks. In contrast to that, intrinsic contributions depend only on the sterical and electronic preferences of the framework itself and are tightly connected to the conformational freedom of the building blocks. It is an important aspect of hybrid porous materials like MOFs that these intrinsic effects can dominate, which allows to alter the topology by substituting the organic linker.28 As we have shown recently for a MOF,29 the latter thermodynamic preference can be determined theoretically by accurate force field methods since the different topologies represent supramolecular isomers. The challenge is to model the respective systems with sufficient accuracy, in order to assess the strain energy and its origin with sufficient precision. In principle, quantum mechanical methods could be used, but treating such large 4922

dx.doi.org/10.1021/jp211280m | J. Phys. Chem. C 2012, 116, 4921−4929

The Journal of Physical Chemistry C

Article

Figure 3. Nonperiodic reference systems used for the force field parametrization of the trigonal units (hydrogen atoms, white; oyxgen atoms, red; carbon atoms, blue; boron atoms, yellow; silicon atoms, brown).

assumed for the bor and ctn networks (space group P1). The initial structures were generated by the WEAVER code,46 followed by a full energy minimization including the orthorhombic cell parameters. In order to avoid local minima, which are due to conformational flexibility combined with steric strain, we used a simple annealing scheme (minimizations of snapshots, generated by a molecular dynamics simulation in the NVT ensemble at an elevated temperature of 600 K). Molecular dynamics simulations in the NPT ensemble were performed to investigate the thermal properties, with an integration time step of 1.0 fs, using the Berendsen coupling method for the barostating47 and a Nose−Hoover thermostat.48 For the parametrization of the force field, our recently developed Genetic Algorithm (GA) was used.35 We performed 20 independent GA runs, with 100 generations and a population of 200. The best individuals of each run were collected for a final GA run. The crossover rate and mutation rate were 0.85 and 0.05, respectively. In certain cases, the relative energy of rotational transition states was included into the objective function. Solvent accessible surface (SAS) areas were computed using a code provided by Düren et al.49,50 Grand Canonical Monte Carlo (GCMC) calculations were performed to predict the carbon dioxide adsorption in various three-dimensional COFs. For every GCMC simulation, 50 000 equilibration cycles were done followed by 100 000 production cycles. The GCMC simulations include random insertion/ deletion, translation, and rotation moves of CO2 with equal probabilities. All isotherms were computed at T = 298 K up to 50 bar, and the discussed excess quantities were obtained using a Lennard-Jones potential (cutoff distance of 12.5 Å) and the Ewald sum technique for the electrostatic contributions. The force field parameters are given in the Supporting Information. During the GCMC runs, the framework structures were kept fixed at the optimized geometries. In order to calculate the fugacities, the Peng−Robinson equation was applied.

Figure 4. Nonperiodic reference systems used for the force field parametrization of the tetrahedral units D[X] (X = C, Si) (hydrogen atoms, white; silicon and carbon atoms, brown).

Both the minimum structures and the corresponding second derivatives are fundamental for the force field parameter fitting approach since they allow a proper extrapolation of the PES in the accessible energy range around the reference structures, and their convergence with respect to the theoretical level needs to be monitored carefully. A crucial point in the parametrization is the fact that in the two topologies (Figure 1), different local symmetries of the tetrahedral units are present. In the ctn net, the four coordinate vertex has a local S4 symmetry, whereas in bor, it has a D2d symmetry. In the real network, this is realized by two different conformations of D (Figure 4). Already in our previous work, we found it to be necessary to refit the standard MM3 parameter set in order to properly model the tetraphenyl molecule.51 In order to accurately describe the PES for the conformational freedom with respect to single bond rotations, we now decided to consider in addition also the rotational transition states of all the reference systems. In this context, we need to account for dispersion effects, which are not covered by standard DFT functionals. The importance of noncovalent intramolecular interactions was shown before.56,57 In the modeling of porous hybrid materials, dispersion contributions were just considered in the context of weak host−guest interactions like, for instance, in the case of hydrogen physisorption. Consequently, mostly wave function-based correlated methods are used in parametrizing these intermolecular potentials.7,16,19,21,58 However, to the best of our knowledge, these effects were not considered with respect to their influence on intramolecular steric interactions in porous materials. Using the structures obtained on the B3LYP level, improved relative energies were computed by single point calculations using the double hybrid functional B2PLYP. In addition, for the tetrahedral SBUs, the MP2 method was used for a full geometry optimization. In case of the models for the trigonal SBUs, a rotation around one of the carbon−boron single bonds was considered. The QM results of the trigonal units are described in the Supporting Information. For the tetrahedral models D[X], the concerted X−C bond rotation from S4 to D2d symmetry was

3. RESULTS AND DISCUSSION 3.1.1. Force Field Parametrization. Reference Systems. Motivated by our previous work35,36,51,52 and experimental concepts,10,53−55 we used a building block approach to generate appropriate reference systems for the force field parametrization. Thus, the necessary parameters for the currently known boron based 3D-COFs 102, 103, 105, 1082 and 2023 have been derived from six nonperiodic model systems shown in Figures 3 and 4. These fragments are derived from the SBUs of the real systems by just cutting carbon−carbon bonds. Thus, in the case of the trigonal SBUs, the boron atoms are always saturated with phenyl groups. Because of the size of the trigonal unit in COF105 and -108, we further divided it into two subunits (B1 and B2; see Figure 3). 4923

dx.doi.org/10.1021/jp211280m | J. Phys. Chem. C 2012, 116, 4921−4929

The Journal of Physical Chemistry C

Article

computed. Because of the relevance for the network strain energy, these results are discussed here in detail. For the tetraphenylsilane D[Si], the influence of basis set and theoretical level is not significant (Table 1). In contrast to Table 1. Energy Difference between the S4 and D2d Conformer of Tetraphenylmethane D[C] and Tetraphenylsilane D[Si] (in kcal/mol) D[C] cc-pVDZ cc-pVTZ cc-pVQZ CBS force field

D[Si]

B3LYP

B2PLYP

MP2

B3LYP

B2PLYP

MP2

0.69 0.82 0.84 0.85

1.13 1.11 1.09 1.06 1.10

1.76 1.56 1.49 1.45

1.84 1.83

1.91 1.87

1.85

Figure 5. Nonperiodic model system B used for the force field validation.

fragments. The results for the complete molecule B, using the combined force field fitted to B1 and B2, were compared with the corresponding reference data computed on the B3LYP/ccpVDZ level. The geometry is very well reproduced with deviations in the bond lengths below 1% (see Supporting Information). More importantly, also the normal modes below 400 cm−1, which are mainly collective motions of the whole system, like the out of plane deformation or twisting of the B2 fragments with respect to the B1 unit, compare very well (see Table S23, Supporting Information). For the normal mode describing the buckling of the system, we see excellent agreement between DFT reference (17 cm−1) and force field (18 cm−1) and an overlap of the normal mode vectors of 1.0. This confirms our decision to partition the periodic system into nonperiodic subunits for the parametrization. A similar accuracy of our force field was found for the butyl-substituted C as in COF-202 (see Figure S5 and Table S25, Supporting Information). 3.2.1. Prediction of Network Structures and Properties. Topologies and Strain Energy Analysis. It is striking that four of the five known boron containing 3D-COFs are experimentally obtained in the ctn topology. This holds also for the later published COF-202, which is unique in the sense that the phenylene units stand perpendicular to the trigonal vertex plane (see also model C in Figure 3), whereas in all other COFs, they are parallel to this plane. Only COF-108 shows the less dense bor topology, which differs from COF-105 only by substituting silicon with carbon in the center of the tetrahedral unit. However, from the first-principles reference calculations, we have seen that this substitution leads to noticeable differences in the conformational behavior and that the model D[C] was the hardest to parametrize accurately. Thus, it might indeed be possible that COF-108 is special in its conformational flexibility. However, the synthesis conditions for the preparation of COF-108 are substantially different from the silicon analogue COF-105 (different solvent mixtures, ratios of the monomers and crystallization times),2 and the observed structure could also be due to extrinsic factors. In order to shed light on this question, our well parametrized force field can be used to investigate structure and strain energy for the ideal frameworks in both topologies in the absence of any guest molecules (in other words, perfectly activated). Note, that especially for the larger systems COF-105 and -108, the experimentally obtained samples were synthesized under different conditions and could not be fully activated, resulting in a coloring of the materials due to the remaining oxidized monomer. Also a Brunauer−Emmett−Teller (BET) analysis could not be performed for them.2 At this point, it should be

1.90

this, the most difficult reference system was found to be the tetraphenylmethane D[C]. Replacing the central silicon by a carbon atom causes a decrease of the rotational barrier by about 1 kcal/mol to only 0.69 kcal/mol on the B3LYP/cc-pVDZ level, despite the shorter C−C versus Si−C bond lengths, showing the different nature of the D[C] and D[Si] units. In addition, an increase of the rotational energy barrier is observed when enlarging the basis set, converging at a value of 0.84 kcal/ mol for B3LYP. Furthermore, dispersion effects are very important for an accurate description of the conformational energies of D[C]. Since the MP2 values show a very slow convergence, we exponentially extrapolated the cc-pVDZ, ccpVTZ, and cc-pVQZ results to the complete basis set (CBS) limit. The MP2/CBS value of 1.45 kcal/mol is about 0.6 kcal/ mol higher than the B3LYP result. However, the use of the double hybrid functional B2PLYP, gives an intermediate value of about 1.1 kcal/mol for the barrier and converges much more quickly with the basis set size as mentioned before.40 Because of the known overestimation of MP2 for aromatic interactions59 and the good description of the B2PLYP functional for the isomerization energy of organic molecules,60 the double hybrid DFT method is considered as the most reasonable choice for the particular system. 3.1.2. Parameter Fit and Validation. In the case of A, B1, B2, and C, the minimum structure and the Hessian matrix computed on the B3LYP/cc-pVDZ level were used to fit the parameters by the GA optimization. In addition, to be able to reproduce the rotational profile around the carbon−boron bond, we found it sufficient to adjust the corresponding torsion parameter manually to reproduce the B2PLYP/cc-pVTZ energy barrier afterwards, since this term couples only very weakly. In the case of D[C] and D[Si], we decided to use the rotational energy barriers between conformers calculated on the B2PLYP/cc-pVTZ level (or B2PLYP/CBS where available) explicitly in the target function of the GA optimization. Overall, we obtained a very good correspondence between the final force field and the reference first-principles data for both structure (mean deviation: bond,