Tongan Yan,Minman Tong,Qingyuan Yang,Dahuan Liu,Yandong Guo,Chongli Zhong

1 State Key Laboratory of Organic-Inorganic Composites,Beijing University of Chemical Technology,Beijing 100029,China

2 School of Chemistry and Materials Science,Jiangsu Normal University,Xuzhou 221116,China

3 College of Mathematics Science,Bohai University,Jinzhou 121013,China

4 State Key Laboratory of Separation Membranes and Membrane Processes,Tiangong University,Tianjin 300387,China

Keywords:Metal-organic frameworks Open Cu sites Molecular simulations Carbon dioxide Diffusion

ABSTRACT Understanding CO2 diffusion behavior in functional nanoporous materials is beneficial for improving the CO2 adsorption,separation,and conversion performances.However,it is a great challenge for studying the diffusion process in experiments.Herein,CO2 diffusion in 962 metal-organic frameworks (MOFs)with open Cu sites was systematically investigated by theoretical methods in the combination of molecular dynamic simulations and density functional theory (DFT) calculations.A specific force field was derived from DFT-D2 method combined with Grimme’s dispersion-corrected (D2) density functional to well describe the interaction energies between Cu and CO2.It is observed that the suitable topology is conductive to CO2 diffusion,and 2D-MOFs are more flexible in tuning and balancing the CO2 adsorption and diffusion behaviors than 3D-MOFs.In addition,analysis of diffusive trajectories and the residence times on different positions indicate that CO2 diffusion is mainly along with the frameworks in these MOFs,jumping from one strong adsorption site to another.It is also influenced by the electrostatic interaction of the frameworks.Therefore,the obtained information may provide useful guidance for the rational design and synthesis of MOFs with enhanced CO2 diffusion performance for specific applications.

1.Introduction

With the continuous increase of anthropogenic CO2emissions,the atmospheric CO2content is at a high level,which has been widely considered as the primary cause of the increase in global temperature[1-5].Therefore,carbon capture,utilization,and storage (CCUS) technologies received extensive attention in order to provide a solution to both carbon emission-control and energysupply challenges [6].The key to CCUS is the advanced functional materials with excellent adsorption,separation and conversion performance of CO2[7].In recent years,metal-organic frameworks(MOFs)have exhibited attractive CCUS properties due to their ability to adjust the porosity,pore type,pore size and functionality by combining different metal nodes,organic ligands and synthesis strategies [8-18].In particular,MOFs with open Cu sites have shown great potential applications in adsorption [19-21],separation[22-24],and catalysis[25-27].Despite the amount of research on CO2adsorption and separations and catalysis using MOFs,the study on CO2diffusion is very insufficient up to now,which is the important factor affecting the overall performance of the above applications [28,29].For example,for the pressure swing adsorption process for CO2/CH4separation using the microporous Cu-MOF(Cu(hfipbb)(h2hfipbb)0.5),it is found that the difference in diffusion rates between guest molecules brings the possibility of kinetic separation of CO2/CH4mixture [30].In the threedimensional (3D) hierarchical electrode composed of d-Cu with abundant active sites prepared by Cu-MOF,the attainment of a higher current density of formate can be obtained by improving the diffusion of the CO2reactant [31].Liuet al.[32] prepared a Cu-MOF/NP catalyst that has high electroreduction of CO2to CO,and it is observed that the reduction of gas diffusion resistance can significantly enhance the catalytic activity.Therefore,the lack of fundamental understanding of CO2diffusion behavior tremendously limits the development of advanced materials for CCUS.One of the main reasons is that it is a great challenge to explore the dynamics of CO2in experiments,especially for examining a large number of potential materials.In this aspect,the largescale molecular simulation is a powerful tool that can gain insight into the dynamics of CO2in MOFs.Unfortunately,to the best of our knowledge,there are few researches on systematic investigations of CO2diffusion in MOFs at the moment,especially for MOFs with open Cu sites.

Herein,962 MOFs with open Cu sites were collected from the Cambridge Structural Database (CSD) and used for CO2diffusion at room temperature and pressure.Generally,the chemical environment imposed by an open metal site is very complex and could not be well described using the generic force fields,such as universal force field (UFF) [33].Thus,to accurately calculate the interaction of CO2molecules in the frameworks,a specific force field was derived from density functional theory (DFT) calculations.DFT-D2 method combined with Grimme’s dispersion-corrected density functional (D2) was utilized to optimize the Cu(II) paddle-wheel cluster to calculate the interaction energies between Cu and CO2at different distances,followed by the fit of discrete energy values of Morse potential function to obtain specific force field parameters.On the basis of the results obtained in large-scale simulations,the relationships between CO2diffusion behavior and geometric and chemical properties of materials are discussed,which may provide useful information for better understanding the diffusion of CO2in MOFs with open Cu sites.

2.Materials and Methods

2.1.MOF structures

In this work,the material database constructed in our previous work [34] was studied,which includes 962 MOF materials with open Cu sites.It is based on the CSD,collected by the algorithm developed by our research group.For several structural problems directly collected from CSD,including disordered structure,solvent occupation and the missing hydrogen atoms,the visual inspection and manual modification were performed referring to the original literatures to build a fully activated structure database,according to the structural processing principles proposed in Ref.[35] when building the computation-ready,experimental MOF database.Geometry optimizations were then carried out using the Forcite module in Materials Studio.The position of disordered atom was allowed to relax while the simulation lattice was fixed at the reported experimental constants.The process of removing solvent molecules simulates the experimental ‘‘a(chǎn)ctivation” process that is guaranteed to produce open metal sites(OMS)in MOFs.The influence of solvent molecules on the structural integrity of the framework is not considered to ensure the generation of the largest data set,which also helps to find promising MOFs[36].The open-source software Zeo++version 0.3[37]was used to calculate pore limiting diameter (PLD),largest cavity diameter (LCD),accessible surface area (Sacc),and free volume (Vfree).SaccandVfreewere calculated by using a probe with a diameter equal to N2kinetic diameter(0.364 nm) and 0.0 nm,respectively.Considering the molecular size,917 MOFs with open Cu sites and PLD greater than 0.330 nm (kinetic diameter of CO2) are included in the following calculations.

2.2.Force fields

Van der Waals force(vdW)and Coulomb potential were used to describing the interactions between materials and guest molecules.Van der Waals effect was expressed by Lennard Jones (LJ)equation.The required LJ parameters between different atoms are calculated by Lorentz Berthelot mixing rule.CO2molecules were adopted a three-point charged LJ model and had linear rigidity.The elementary physical model 2 (EPM2) force field [38] provided the energy parameters,and the C-O bond length was 0.1149 nm.

DFT calculations based on the quantum mechanics level were used to determine the precise interaction parameters between materials andgas molecules[39].Utilizing DMol3moduleofMaterialsStudio software (version 7.0),the GGA exchange functional Perdew-Burke-Ernzerhof (PBE) and the double numerical basis set with polarization(DNP)was used to optimize the model cluster obtained in Cu-BTC(BTC=benzene-1,3,5-tricarboxylic acid)and calculate the interaction energy between Cu and CO2at different distances based on the Grimme’s dispersion correction (D2) [40].In our previous work,the force field describing the special interaction between open Cu sites on MOFs and CO2had been developed,and the accuracy and reliability of the force field had been fully verified[41].

In this work,this force field was further used to explore the diffusion of CO2.Therefore,the vdW interactions between the atoms of CO2and MOFs (except Cu) were described using LJ potential with the parameters taken from UFF force field,while the specific Cu-O interactions were described using the Morse potential.The equations of the above potentials are defined as:

whereUis the potential energy of the interaction betweeniandj,andqi(qj) is the partial charge on the atom.For the LJ potential parameters,σ is the well depth,ε is the collision diameter,andrijis the atomic distance.For the Morse potential parameters,α is the potential well curvature,ris the distance between the Cu atom on the cluster and the O atom on one side of CO2,D0andR0represent the minimum energy and the distance corresponding to the minimum energy,respectively.

2.3.Simulation details

First,grand canonical Monte Carlo (GCMC) simulations [42]were employed to calculate the adsorption behaviors of CO2in MOFs at 298 K and 0.1 MPa.During the simulations,the molecules involve four types of trials including random replacement,random regeneration,new creation and existing deletion.All the atoms in the MOFs were frozen in their crystallographic positions to ensure that the frameworks are rigid.Part of the point charges on the frameworks were distributed by the extended charge balance method(EQeq)[43].In order to avoid the finite-size effect,several large simulation boxes were used to verify the simulation results,and the periodic boundary conditions in three-dimensional space were fully considered.The number of MOFs unit cells in the simulation boxes ranges from 1×1×1 to 5×3×3.For LJ interaction,the cut-off distance was set to 1.40 nm.The Ewald summation technique was used to deal with the long-range electrostatic interaction.Periodic boundary conditions are considered in threedimensional space.Fugacity in the simulations was transformed by Peng-Robinson equation of state.For each state point,GCMC simulation includes 1 × 107equilibrium steps and 1 × 107sampling steps to obtain the thermodynamic properties.

Then,based on the results of GCMC simulations,the canonical(NVT)ensemble was used to simulate the diffusion of CO2in MOFs with open Cu sites.Molecular dynamic(MD)simulation conditions were 298 K and 0.1 MPa.Nosé-Hoover chain (NHC) [44] thermostat was used to maintain constant temperature conditions,and velocity Verlet algorithm was used to integrate Newton’s equation of motion.The simulation steps and a time step of each MD simulation were 6 ns and 1 fs,respectively,and then preceded by 3 ns of equilibration.Finally,10 independent trajectories were output to calculate the self-diffusion coefficients on average.Our internal code HT-CADSS supported all simulations.

2.4.Properties calculations

The CO2molecules are tagged with Einstein equation [45] to obtain the mean square displacements (MSD),and the selfdiffusion coefficient is further determined by taking the slope calculation for a long time:

where 〈···〉 is the ensemble average,Nis the number of molecules,andri(t) describes the center-of-mass (COM) position vector of moleculeidetermined by timetin diffusion.dcorresponds to the size of the inspected system.For MOFs with three-dimensional pores,the average diffusion coefficient is taken in thex,yandzdirections,whered=3;when the material has one-dimensional pores,only thezdirection is considered,that is,d=1.

2.5.DFT calculations

B3LYP functional [46] was used to optimize the positions of H atoms used for frame clusters passivation.The SDD basis set [47]was used to describe the effective core potentials of Cu atoms,while the 6-311+G (d,p) basis set [48] was used to describe other atoms.This method has been proven effective in describing geometric structure,electronic distribution and weak interaction[49,50].The thermodynamic properties were obtained by frequency calculation,and optimizations were verified to reach the minimum energy value.The Gaussian 09 software [51] was used for DFT calculations,and the GaussView version 5.0.8 software[52] was used to obtain the three-dimensional structure model.

3.Results and Discussion

3.1.Validation of simulation method

Before the large-scale simulations,it is necessary to verify that the force field parameters used can accurately describe the interaction between adsorbent and adsorbate.Though many quantum calculations have been performed to explore the interactions between CO2and MOFs with open Cu sites,it is rarely reported for the related efforts paid on developing accurate force fields for further implementation of classical molecular simulations [53].In our previous work [41],a force field for the special interaction between open Cu sites on MOFs and CO2has been established,which has been verified to be universally applicable in the simulation for adsorption behaviors.However,to the best of our knowledge,CO2diffusion in MOFs with open Cu sites using specific force fields has not been exploited computationally.First,the MD simulation results of CO2were compared with the reported values of Cu(hfipbb)(H2hfipbb)0.5[54] and IISERP-MOF20 [55] in literature,and the results are shown in Fig.1.It can be seen that our calculated results using UFF force field are in good agreement with the values obtained by Watanabeet al.[54] and Nandiet al.[55],ensuring the reliability of the MD method used in this work.As we know,UFF force field cannot accurately describe the specific interaction between the guest molecule and the open metal sites contained in the MOFs [33].To further test the reliability of the developed force field for diffusion,the diffusion of CO2in Cu-BTC[28] at 5 × 10-4MPa and 298 K was simulated by using the UFF force field and the DFT-D2 derived FF.As shown in Table 1,the self-diffusion coefficient (Ds) of CO2is significantly reduced after considering the Morse potential parameters for the special interaction between open Cu sites and CO2in the newly derived force field,which is closer to the experimental data than that of UFF force field.In the diffusion process under low-pressure conditions,the main driving force is the interaction between CO2and the preferential adsorption sites (for example,open Cu sites).Therefore,the developed force field can be reliably verified to describe CO2diffusion behavior in these MOFs.

Table 1Comparison of the calculated and experimental coefficient of CO2 diffusion with different force fields

3.2.Relationships between diffusion and structure

It is important for understanding the influence of MOF structures on the CO2diffusion behaviors to guide the design and synthesis of promising candidates with specific performances.Fig.2 shows the relationship between the self-diffusion coefficient and adsorption capacities of CO2in 917 MOFs with open Cu sites.It can be seen that all of theDsof CO2as single-component gas are in the range of 10-7to 10-11m2·s-1.With the increase of CO2adsorption capacities,Dsis gradually decreased.In the case of high adsorption capacities,the diffusion space will be compressed for CO2molecules.As a result,the diffusive rate of CO2in the frameworks decreases since the interactions between adsorbate-adsorbate are increased.

In general,the MOFs studied in this work can be classified into 2D-MOFs and 3D-MOFs.From Fig.2,there is no significant difference in theDsof CO2in 2D-MOFs and 3D-MOFs.However,it is worth noting that the points of 3D-MOFs are more concentrated in the figure compared with those of 2D-MOFs,showing a more evident relationship.For 3D-MOFs,there is an obvious trade-off effect between CO2adsorption capacity andDs.The enhanced adsorption may generally decrease the diffusion of CO2,andvice versa.For example,3D-MOFs MOF-399 has the highestDs(1.39 × 10-7m2·s-1),while the CO2adsorption capacity is only 1.41 mmol·g-1.However,for Cu-MOF-1 (CSD refcode:BEXSII),Dsand CO2adsorption capacity are 7.51 × 10-10m2·s-1and 11.23 mmol·g-1,respectively.For 2D-MOFs,such a trade-off effect is not evident,and it is difficult to find a clear correlation trend between adsorption capacities andDsof CO2.For example,the CO2adsorption capacities of Cu_OBut-bdc (CSD refcode:GEWMON) and MCF-21 (CSD refcode:XOPLIY) are similar (5.08 and 4.96 mmol·g-1),while there is a large difference inDs(3.30×10-9and 9.74×10-11m2·s-1).For Cu-MOF-2 (CSD refcode:IGOJAR)and Cu-Tbip (CSD refcode:LEDCUU),Dsvalues are close(2.85 × 10-10and 2.36 × 10-10m2·s-1),while with different CO2adsorption capacities (1.90 and 6.20 mmol·g-1).This is maybe due to the unique and complex structural characteristics of 2DMOFs.For the 189 kinds of 2D Cu-MOFs studied in this work,180 of them have small apertures with the PLD and LCD less than 1.0 nm.In addition,the interlayer space also plays a significant role,which has been observed in experiments [56].Therefore,it is more flexible in tuning and balancing the CO2adsorption and diffusion behaviors in 2D-MOFs than in 3D-MOFs to meet the practical demands.

Fig.3 describes the correlations between the structural properties (PLD,LCD,andSacc) of MOFs and theDsof CO2.It can be seen that the pore size of MOFs plays a significant effect on CO2diffusions.The larger the pore size,the faster CO2diffuses,as shown in Fig.3(a).When the PLD is less than 0.66 nm,theDsvalue changes drastically with the change of pore size.While for the materials with PLD larger than 0.66 nm,the increment ofDsbecomes not evident.This may be caused by the fact that the collision of CO2molecules is great in the material with a small PLD,which is gradually weakened with the increase of PLD.As can be seen from Fig.3(b),theDsgradually increases with the increase of LCD,similar to that of PLD.

Fig.1.Comparison of CO2 diffusion coefficients calculated in two typical Cu-MOFs:(a) Cu(hfipbb)(H2hfipbb)0.5;(b) IISERP-MOF20.

Fig.2.Relationships between the CO2 uptake and Ds of CO2(Blue circles:2D-MOFs,red circles:3D-MOFs,green lines:trade-off relationship curve based on data of 3DMOFs).

As we all know,PLD is the maximum accessible guest molecule size,and LCD describes the largest pore size distribution [57].Because CO2molecules mainly diffuse in the pores of materials,LCD is more helpful to clarify the diffusion behavior.According to Fig.3(b),the diffusion behavior of CO2can be roughly divided into four stages with the change of LCD.(a)When LCD is less than～0.66 nm,only one CO2molecule is preferentially located in the channel in consideration of its kinetic diameter (0.33 nm).From the COM maps in Fig.4(a),the open Cu sites on both sides of the frameworks can synergistically interact with the same CO2molecule,making it difficult to diffuse (Dsis less than 10-9m2·s-1).(b) When the LCD is in the range of～0.66-0.99 nm,the channel can hold two CO2molecules.As shown in Fig.4(b),one CO2molecule can form a strong interaction with the open Cu sites on the one side of the channel,while the other CO2molecule interacts with the open Cu sites on the other side of the channel.As a result,CO2diffusion rate is enhanced.However,the steric hindrance effect is still significant,leading to the strong dependence ofDsvalue on PLD.(c) When the LCD is in the range of～0.99-1.32 nm,the third CO2molecule can be accommodated between two CO2molecules interacting with single open Cu site.Due to the enhanced interaction with these two CO2molecules,the diffusion will be limited.(d)When the LCD is larger than 1.32 nm,the steric hindrance effect is not evident because the pore size is large enough.CO2molecules are mainly affected by the open Cu site on the unilateral frameworks,inducing that theDsvalue changes not too much with further increasing the pore size.

Fig.3.Relationships between PLD,LCD, Sacc and Ds of CO2.

Fig.4.Two-dimensional color-contour plots of the center of mass(COM)probability densities of adsorbed CO2 in Cu-MOFs with different pore sizes by taking two MOFs as examples:(a) Cu-MOF-2 (CSD refcode:IGOJAR),LCD=0.447 nm;(b) Cu-MOF-3 (CSD refcode:MOYYEF),LCD=0.767 nm.

Table 2 lists top 14 MOFs with open Cu sites with largeDsvalues of CO2.It can be seen that 5 of them are tbo topology and 4 of them are pto topology.In essence,both tbo and pto are 4-3 connected,and these two topologies have high consistency in physical properties and can be obtained by the same ligand and Cu(II)paddle-wheel using different synthesis conditions in experiments[58].Therefore,it is believed that these kinds of topologies may be conducive to the CO2diffusion,providing reference information for screening materials from the CO2diffusion point of view.

3.3.Residence time of CO2 in MOFs

From Fig.3,it is observed that the increase ofDsof CO2may encounter a bottleneck when the LCD is larger than～1.0 nm.Among all 580 Cu-MOFs,only 14 of them(Table 2) exhibits a highDsvalue exceeding 3.0×10-8m2·s-1.Therefore,besides the structure of MOFs,it is worthy to further investigate the effect of chemical environments in pores on CO2diffusion.To this target,MOF-399,DUT-63 and PCN-6′(CSD refcode:NIBJAK) were selected(Fig.5),which have the same topology (tbo),similar CO2adsorption capacities (1.40,1.33 and 1.21 mmol·g-1) with the relatively close pore sizes (PLD:2.424,1.954,and 1.755 nm,LCD:4.279,3.753,and 3.199 nm).Despite the similarity in structure,there is a significant difference in theDsof CO2.For MOF-399,it is 1.39 × 10-7m2·s-1,3.65 times that of DUT-63 (3.82 × 10-8m2·s-1) and 18.25 times that of PCN-6′(7.61 × 10-9m2·s-1).This may be resulted from the different atoms in the frameworks of these three MOFs,as shown in Fig.5.

Table 2Properties and structural parameters of the top 14 Cu-MOFs with the large diffusion coefficients

To understand these results in-depth,the CO2diffusion movies obtained by MD simulations were decomposed frame by frame to help us study the diffusion behavior of CO2in these materials.At the first step,we counted the residence time of CO2at different positions in the frameworks,that is,open Cu sites,the center of organic ligand(COL),other positions of frameworks(OPF)and pore space(PS).It can be seen from Fig.6(a)that for the total statistical time of 3000 ps,the residence time of CO2at open Cu sites in MOF-399,DUT-63 and PCN-6′is about 1520,1520 and 1540 ps,respectively.The residence time near COL is 760,980 and 1300 ps,respectively.The sum of these two kinds of residence time accounts for more than 75% of the total time.It is interesting to found that in this kind of 4-3 connected 3D-MOFs,CO2molecules first adsorbed on the preferential adsorption sites,open Cu sites,and then jumped and diffused between these sites (Fig.6(b)).The organic unit on the frame can act as a ‘‘springboard”,which makes the movement of CO2in MOFs more directional to ensure that molecules can diffuse from one open Cu site to another open Cu site quickly.

However,if CO2molecules stay near the organic unit for a long time,it will play a negative effect on diffusion.As shown in Fig.6(a),the molecule only stayed for 760 ps near the ‘‘springboard” in MOF-399 with the fast CO2diffusion (Ds=1.39 × 10-7m2·s-1).In contrast,this residence time was extended to 1300 ps in PCN-6′,which is almost the same as that at open Cu sites positions(1540 ps).Thus,theDsvalue is sharply reduced to 7.61 × 10-9m2·s-1.According to the total simulation time and the time step,the diffusion movies obtained from MD simulations are decomposed.Furthermore,the distance between CO2and the typical atoms in the frameworks was counted one by one according to the radial distribution function (RDF) data of each simulation frame.Taking Cu and N as indicators,they are denoted as CO2_OCu and CO2_O-N.Fig.7 shows the change of atomic distance as a function of diffusion time.As can be seen from Fig.7(a) and (c),both of the shortest distances of CO2_O-Cu in DUT-63 and PCN-6′are about 0.30 nm,demonstrating no significant difference in the role of open Cu sites during the diffusion of CO2in these two MOFs.While for the distance of CO2_O-N,the range is 0.30-1.70 nm in DUT-63,which is slightly reduced to 0.30-1.0 nm in PCN-6′(Fig.7(b) and (d)).This is maybe attributed to the large number of N atoms in PCN-6′comparing to that in DUT-63,exhibiting a relatively strong interaction with CO2.

3.4.Influence of electrostatic interactions on diffusion behavior

From the previous works [59],it is observed that electrostatic interaction plays an important role in the adsorption of CO2in MOFs.Therefore,we further investigated the influence of such interaction on diffusion behavior,by switching off the charges of the atoms in frameworks for additional simulations.As shown in Fig.8(a),theDsvalue of CO2in DUT-63 is increased from 3.82 × 10-8to 7.83 × 10-8m2·s-1.For PCN-6′,the increment is more evident,from 7.61 × 10-9to 7.60 × 10-8m2·s-1.The residence time of CO2at each position in MOFs was also counted,and the results are listed in Fig.8(b).Obviously,it is greatlyreduced for open Cu sites and becomes almost equal at each position of the frameworks.In contrast,it is still very short in pore space.That is,CO2molecules still diffuse along with the frameworks,while no jumping behavior occurs.This directly causes theDsvalue to decrease from 1.39 × 10-7to 7.45 × 10-8m2·s-1in MOF-399.On the other hand,without the interactions of Ncontaining ligand,CO2molecules can diffuse more easily along with the frameworks,leading to the increase ofDsin DUT-63 and PCN-6′.

Fig.6.CO2 diffusion movies analysis:(a) statistical chart of residence time;(b) schematic diagram of diffusion path of CO2 in MOFs.

In order to better understand this phenomenon,the free energy profile for CO2is further calculated at 77 K and infinite dilution.The interaction coordinateqis the center of mass position of the molecule along with the frameworks,depicted as the red line on the right of Fig.9.Along with the frameworks,the fluctuation of free energy in PCN-6′with the lowest CO2diffusion coefficient is much more intense than that in MOF-399 and DUT-63.It is the smallest in MOF-399 among the three MOFs,inducing the best CO2diffusion ability.The results indicate that different organic units will lead to energy heterogeneity on the frameworks,which in turn leads to the difference in CO2diffusion capacity.

In addition,molecular electrostatic potential (MESP) is used to study the interaction between CO2and frameworks,in order to further explore the effect of different organic units on CO2diffusion.As shown in Fig.10,the O atoms at both ends of CO2generate a negative potential region,while the C atom in the middle generates a positive potential region.Previous studies have shown that the O atom at one end of the CO2molecule will form a strong interaction with Cu on the MOFs,while the O atom at the other end is exposed[53].It should be noted that a strong positive potential region is formed near the organic ligands on the MOF-399 frameworks,and a strong negative potential region is formed in PCN-6′.Accordingly,the diffusion coefficient in the former structure is much larger than that in the latter.Since the organic unit of PCN-6′contains more N atoms than that of DUT-63,the electronegativity is strong,leading to the strong interaction of the electrostatic-sensitive CO2guest molecules.During the diffusion of CO2molecule in MOF-399,the O atoms at both ends will be attracted by the open Cu sites and organic ligands respectively,resulting in ‘‘jumping” diffusion behavior.However,in PCN-6′,the O atoms at both ends of the CO2molecule are attracted by the open Cu sites and repelled by the organic ligands respectively,thus exhibiting a ‘‘holding” diffusion behavior.

Fig.7.Relationships between atomic distances range and diffusion time:(a)CO2_O-Cu;(b)CO2_O-N.Mean radial distribution function of diffusion:(c)CO2_O-Cu;(d)CO2_ON.

Fig.8.(a)Simulations of mean square displacement diagram with and without electrostatic forces between frameworks and CO2 molecules.(b)Statistical chart of residence time without frameworks charge.

The calculated binding energy of the CO2molecule with the clusters of PCN-6′(34.05 kJ·mol-1) is larger than that of DUT-63(24.21 kJ·mol-1) and MOF-399 (18.59 kJ·mol-1).Thus,a strong potential energy well can be formed to ‘‘hold” CO2molecules within a relatively short distance for a long time.The diffusion is hindered to result in a lowDsvalue.In short,CO2has better diffusion ability in MOFs with weakly electronegative organic ligands,which should have an internal correlation with the energetic environment in the frameworks.

Fig.9.(a)Free-energy profile for CO2 in 3 MOFs along with the frameworks at 77 K and infinite dilution,the yellow area on the left shows the free energy profile along the interaction path.(b) The red line in the right figure represents the interaction path.

Fig.10.Molecular electrostatic potential diagram of 3 MOFs clusters and CO2 molecule.

4.Conclusions

In this work,CO2diffusion behavior in 962 MOFs with open Cu sites was systematically studied using large-scale MD simulations.With the aid of the specific force field derived from DFT calculations,it is suggested that tbo and pto topologies are conductive to CO2diffusion and 2D-MOFs are more flexible in tuning and balancing the CO2adsorption and diffusion behaviors than 3D-MOFs.Through careful analysis of diffusive trajectories and the residence times on different positions,it is observed that CO2diffusions in these MOFs are mainly along with the frameworks,jumping from one strong adsorption site to another.The electrostatic interaction of the frameworks also plays an important role in the diffusion process.Therefore,this work helps to understand the CO2diffusion behavior in MOFs with open Cu sites by presenting a computational evaluation.Obtained information may provide useful guidance for the rational design of MOFs with enhanced CO2diffusion performance.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work is supported by National Key Research and Development Program of China(2016YFB0600901)and the Natural Science Foundation of China (22038010,21878229,22078024 and 21978005).