key: cord-0058455-4rtvud2p
authors: de Aragão, Emília Valença Ferreira; Faginas-Lago, Noelia; Apriliyanto, Yusuf Bramastya; Lombardi, Andrea
title: Gas Adsorption on Graphtriyne Membrane: Impact of the Induction Interaction Term on the Computational Cost
date: 2020-08-26
journal: Computational Science and Its Applications - ICCSA 2020
DOI: 10.1007/978-3-030-58820-5_38
sha: 0324ba009850f2c96fdadce948b88ca26daeca9c
doc_id: 58455
cord_uid: 4rtvud2p

Graphynes are a family of porous carbon allotropes that are viewed as ideal 2D nanofilters. In this present work, the authors have modified the Improved Lennard-Jones (ILJ) semi-empirical potential used in the previous works by adding the induction term (iind) to define the full interaction. The evaluation of the computational cost was done comparing ILJ vs ILJ-iind and analyzing the adsorption of 1 gas (CO[Formula: see text]) and a small mixture of gases containing CO[Formula: see text], N[Formula: see text] and H[Formula: see text]O. The computational time of the different calculations is compared and possible improvements of the potential models are discussed.

Recent reports have shown that the concentration of CO 2 in the atmosphere has risen a lot the last few decades [1] . This trend is seen as a consequence of largescale human activity, whether it involves energy production or manufacturing materials (cement, iron, steel, etc) [2] . The excess of CO 2 in the atmosphere causes many problems, such as the more frequent apparition of toxic blue green algae in lakes during hot seasons and the rising of global temperatures [3] . This is why it is urgent to investigate strategies that can be implemented in order to, if not lower the CO 2 in the atmosphere, at least change the evolution trend and keep the CO 2 concentration at the current level. There are two main ways to do it: either capturing CO 2 (and later, the captured CO 2 can be converted into methane and other fuels [4] [5] [6] ). The second approach is generally seen as more efficient and energetically cheaper than the first, but there are some constraints: the flue gas is a mixture of water, carbon dioxide, oxygen and nitrogen molecules [7] . One way to selectively capture CO 2 in flue gas is through the adsorption using porous materials [8] [9] [10] [11] . The advantage of this method is that it is relatively cheaper and simpler to implement on existing power plants [2] . In the last few years, a range of porous materials have been evaluated in their ability to selectively capture CO 2 : a) nanoporous carbons [12] [13] [14] , b) zeolites and zeolitic imidazolate frameworks (ZIFs) [15, 16] , c) metal-organic frameworks (MOFs) [17] , d) porous polymer networks (PPNs) or covalent organic frameworks/polymers (COFs/COPs) [18, 19] , and e) a slurry made of solid adsorbents in a liquid absorbent [20] . In particular, carbon-based membranes have desirable physicochemical properties (e.g. hydrophobic, chemically inert and thermally stable) and are economically suitable and viable for carbon capture and sequestration (CCS) [21, 22] . In contrast, MOFs show permeability and good selectivity, but they are not resistant in the presence of water vapor nor in high temperature, which are usually the conditions of CO 2 combustion. What makes both of these classes effective for gas separation is their permeability and their selectivity. In addition, the thinner a membrane is, the more it is permeable, which makes single-layer membranes interesting objects of study in this context [23] .

In CCS the range of options of applicable materials is vast, so it is impossible to do the synthesis, the characterization and the evaluation of the ability to selectively capture CO 2 for all of the candidates [24] [25] [26] . This is why a preliminary investigation with computational modelling and simulation is crucial for narrowing the selection of molecules. Molecular Dynamics (MD) simulations are a theoretical chemistry method for analyzing the movement of atoms and molecules using potential functions. It can be employed to investigate the structural rearrangement of pure solvents, mixed solutions and combustion processes. [27] [28] [29] [30] [31] [32] . In MD, a set of potential functions is called a force field. Currently, a number of options for force fields are available in MD software, such as UFF [33] and AMBER [34] . However, both these force fields are limited in use: when studying a particular system, those force fields do not always have appropriate parameters to describe it since they are too generic. Therefore, the researcher interested in a particular system has to develop or modify parts of the force field, choosing better potential energy functions. The choice of the modifications must be based on experimental and theoretical data available from the literature or prior quantum chemical computations. While the parameterization of force fields is not a simple task, it is crucial for describing a system correctly. Recently, a number of force fields have been developed for evaluating the adsorption of gas molecules on different porous materials like zeolites [35] , MOFs [36, 37] , graphene and its derivatives [38] [39] [40] and other polymeric materials [21] . These force fields are used to describe molecular interactions between the gases and the porous layer in a quantitative manner, aiming to give predictions of the adsorption dynamics and transport properties of the gases.

The authors of this paper have been recently involved in the study of γgraphynes using MD tools, in particular the development of force fields related to gas adsorption in that class of carbon allotropes. [41] [42] [43] [44] [45] . γ-graphynes are atomic monolayers where the carbon atoms are arranged in a way that two adjacent hexagons are connected through C-C triple bonds. As nanoporous materials, γ-graphynes are interesting candidates for CCS. It has been reported in the literature that the pores are uniformly distributed and adjustable [46] and that they are not prone to form aggregates due to low dispersion forces. In this work, simulations were performed involving a form of γ-graphyne called graphtriyne and a mixture of CO 2 /N 2 /H 2 O in an attempt to reproduce the chemical environment of flue gases. A challenge the authors have been facing is the considerable increase of the computational time by many times because of a choice made in the potential function to represent the interaction between the gases. In this report, a discussion is hold on how this problem is being confronted.

In the next section of the present work, methods and construction of the present potential energy function are outlined. In Sect. 3 a discussion about the improvements of the code is hold and in Sect. 4 the paper brings up concluding remarks.

The MD simulations were performed in simulation boxes with dimensions 72.210 A x 62.523Å x 280.0Å. Inside each box, a graphtriyne membrane with dimensions 72.210Å x 62.523Å was placed. The membrane structure was taken from Ref. [47] , and had been previously optimized through periodic DFT calculations. Simulations were performed uniquely at the temperature of 333 K. Part of the simulations were performed with only CO 2 molecules, while the other part involved a CO 2 /N 2 /H 2 O gaseous mixture. This mixture was composed of an equal number of moles of all molecules. For the nitrogen and the carbon dioxide molecules, models taking into account the quadrupole moment were employed. Those models were a three-charge-site N 2 model [48] and a five-charge-site CO 2 model [49] . As for the water molecule, a model taken from Ref. [50] was used. In this model, the charge is distributed in a way that corresponds to the dipole moment of water in the gas phase (1.85 D) [51] . All the details of the geometries of those molecules are shown in Fig. 1 . For the MD simulation, the membranes were set as a frozen framework while the gas molecules were set as rigid bodies. The gas molecules were randomly distributed with equal amount into each region of the box. Figure 2 shows the relative sizes of the gas molecules, the graphtryine membrane and the simulation box in a qualitative way.

In a MD simulation, defining correctly the intermolecular forces at play is important for obtaining accurate results. In this system, the intermolecular forces of interest are those between gas molecules and between the gas molecules and the graphtriyne membrane. The intermolecular interaction energy is decomposed in terms of molecule-molecule pair contribution, which are electrostatic and nonelectrostatic contributions. The non-electrostatic contributions are measured by taking into consideration the strength of induced dipoles and the average molecular sizes. This can be done by assigning a value of polarizability to the both interacting centers, as shown in Fig. 1 . Here, the intermolecular forces were expressed using the Improved Lennard-Jones (ILJ) potential [52] [53] [54] [55] [56] [57] [58] [59] [60] .

In Eq. 1, ε, r 0 and m are parameters specific to the molecular pair involved, and r is the distance between the two interacting centers of the same molecular pair. In particular, m assumes the value of 6 for neutral-neutral pairs, 4 for ion-neutral pair and 1 for ion-ion pairs. The first term of the Eq. 1 represents the dependence of the repulsion in function of r, while the second term is the dependence of the long-range attraction in function of r. To modulate the decline of the repulsion and the strength of the attraction in Eq. 1, the n(r) term is employed (Eq. 2).

In Eq. 2, β is a factor that modulates the hardness of the interacting pair [61, 62] . This newly introduced parameter is what makes ILJ potential (Eq. 1) able to indirectly take into account some effects of atom clustering, induction and charge transfer and to improve the Lennard-Jones function in the asymptotic region. Nevertheless, while ILJ improves the description of the present system, other effects should be included in the intermolecular potential for enhancing the model of the system. Bearing in mind that charge transfer and induction effects may be important in the interaction between H 2 O and CO 2 and N 2 a careful separate characterization of each contribution was performed. Charge transfer effects in the perturbative limit were taken into account indirectly by lowering the value of β as discussed, for instance, in ref. [42] . In addition, induction due to the permanent water dipole was estimated and incorporated explicitly using the following semiempirical asymptotic expression (in meV)

which is applicable because of the small dimension of water with respect to the related intermolecular distances. In Eq. 3, the left coefficient -2140 (that incorporates the square of the water dipole moment value) is given in meV· A 3 , X i refers each to the C, N or to the O atoms of CO 2 and N 2 , α i is the polarizability (in A 3 ) associated with them, R OW −Xi is the distance vector between the oxygen atom from the water molecule and C, N or O atoms of CO 2 and N 2 and γ is the angle formed by the R OW −Xi vector and the dipole moment of H 2 O. In the present system, the ILJ potential and the electrostatic interactions cutoff distance was set to 15Å. The electrostatic interactions were calculated using the Ewald method, present in the DL POLY 2 code [63] . All molecular dynamics calculations were performed using the DL POLY 2. The system was studied in the canonical NVT ensemble using the Nose-Hoover thermostat and periodic boundary conditions in all directions. At first, two simulations with only CO 2 were run with using ILJ and ILJ coupled with induction. Both lasted 5.5 ns after 0.5 ns of equilibration period with a time step of 1 fs. It was observed that the simulation with ILJ and induction took about 10 times longer than the calculation where only the ILJ potential was involved. Then eight shorter calculations were run in order to analyze why the simulations with ILJ and induction were many times longer than the simulations employing only the ILJ potential. They all lasted 600 fs after 50 fs of equilibration period and a time step of 1 fs. Half of those calculations involved the box containing only CO 2 and the other part was done in the system with the gas mixture. The simulations were performed before and after a modification in one particular routine.

As stated in the introduction, the aim of this work is to compare the time for the different simulations. Table 1 reports the original simulation time before any modification was done to the code. At first the calculations were run with the gas mixture, and it was observed that the total CPU time of the simulation with the ILJ and the induction took about five times longer than the simulation with only the ILJ potential. Then calculations were run with only CO 2 and they were a little bit shorter compared to the ones with the gas mixture, but the simulations with ILJ+induction were still five times longer than with only ILJ. There are no other differences in both ILJ and ILJ+induction options in the code besides the inclusion of the induction part in the latter. This modification was done in a subroutine of the DL POLY2 code that calculates all the interatomic forces using the verlet neighbour list and is written in the forces modules.f file. The forces subroutine of DL POLY2 is a piece of code that calls other subroutines to calculate the intermolecular forces for each atom of the system. The subroutines subsequently called by forces change if one is using ILJ or ILJ+induction. If the calculation involves only the ILJ potential, the ewald1 subroutine is called. If the simulation involves ILJ+induction, then the ewaldm1 and induct subroutines are called. Both ewaldm1 and induct were written by members of the authors group in order to take the induction into account and are situated in ewald module.f and coulomb module.f respectively.

In Table 1 , the time is about the same for every simulation when it starts to go through forces. After the ewald module is called, the time of the ILJ+induction simulations is already 3.9 times higher than the time at the ILJ simulations. Between the ewald module and the coulomb module other subroutines are called, but there is not much difference in the CPU time. The induct subroutine is only used in ILJ+induction simulations and because of it the CPU time becomes five longer than the ILJ simulation. The forces subroutine is called for each of the 600 steps of the simulation. The last two lines in Table 1 report that there is not a significant time difference between leaving the forces module at the last step and the end of the simulation. Evidence points to something inside ewaldm1 and induct subroutines that makes the simulation take more time. At the moment this paper is being written, the ewaldm1 routine is still being analysed. So discussions of the modifications in this subroutine will be left for the future.

The induct subroutine was analysed for both the cases with CO 2 and the gas mixture. As said before, the forces module is called in every step of the simulation, so 600 times. However, the induct subroutine is not called only once in the Forces module. In fact, induct is run for every single atom (including pseudoatoms) in the system at every step of the simulation. In the calculation with only CO 2 , there are 909 atoms in total, so induct is called 545400 times in the simulation. As for the calculation with the gas mixture, the total number of atoms is 972 and the induct subroutine is called 583200 times. Calling that sub- routine one or ten times has reported no difference in the time of the simulation, since this piece of code contains only loops and arithmetic operations. However, when such subroutine is called hundreds of thousands of times in a simulation, it begins to heavily impact in the full CPU time. Moreover, a further analysis of the coulomb module has shown that calling the routine for most of the atoms in the system is useless: there are if statements inside induct that only applies to carbon and oxygen atoms of CO 2 , oxygen atoms of water and nitrogen atoms of N 2 . The proposed solution for turning this task less time consuming was to modify the forces subroutine in order to only call induct when treating those particular atoms. In the code, if statements were written to make the call of induct conditional: if the atom being treated at the moment was O from H 2 O or C from CO 2 , the induct subroutine was called. Else if the atom was O from CO 2 or N from N 2 , the induct subroutine was called. No else statement was made. Table 2 reports the result of this change in the code. Calculations were run again for ILJ and ILJ+induction, with CO 2 and the gas mixture. There were some fluctuations in the simulation times, but the most important is to look at the values of the CPU time after the last induct call cycle. The difference with it and the CPU time before the induct subroutine is less than 0.01 s in both CO 2 and the gas mixture simulations. Total CPU time for ILJ+induction is now around 4 times the Total CPU time of ILJ simulations. When running longer calculations with ILJ+induction, for instance a simulation with 6 million steps, it is expected that the modification made on induct will contribute to spare days of computer time.

In this work, the simulation running time of 600 fs is not long enough for estimating the improvement in the forces results. A simulation on the order of nanoseconds or millions steps is more adequate. In the current version of the code, a simulation of 6 million steps lasts around 20 days when using ILJ+induction, and some technical problems have happened to the cluster during the calculation. Further improvements are needed to be made to the code in order to have a faster running simulation using ILJ+induction, therefore the discussion on the improvement of the forces results will be made in future works. 

A small modification in part of the code has had a considerable impact on the total simulation time. Future modifications in the ewaldm1 subroutine are expected to diminish CPU time even further for calculations with ILJ and induction potentials. While the time will not be the same as simulations with only the ILJ potential, there is at least hope that it will not be 5 or 10 times longer. That way, future simulations with even larger quantities of gas in the box will be possible, and the system will be described better than in simulations with only ILJ.

Climate Change Indicators in the United States: Global Greenhouse Gas Emissions

Carbon capture and storage: introductory lecture

World Resources Institute (WRI): Climate Analysis Indicators Tool (CAIT) 2.0: WRI's climate data explorer

Methane production by CO2 hydrogenation reaction with and without solid phase catalysis

Nanosecond pulsed discharge for CO2 conversion: kinetic modeling to elucidate the chemistry and improve the performance

Fuel production from waste CO2 using renewable energies

Tri-reforming of methane over Ni catalysts for CO2 conversion to Syngas with desired H2/CO ratios using flue gas of power plants without CO2 separation. In: Carbon Dioxide Utilization for Global Sustainability

Evaluating different classes of porous materials for carbon capture

Carbon capture and storage (CCS): the way forward

Porous materials with pre-designed single-molecule traps for CO2 selective adsorption

Atomic and molecular data for spacecraft re-entry plasmas

Exceptional CO2 capture in a hierarchically porous carbon with simultaneous high surface area and pore volume

Activated graphene-derived porous carbon with exceptional gas adsorption properties

Defining a performance map of porous carbon sorbents for high-pressure carbon dioxide uptake and carbon dioxide-methane selectivity

Predicting large CO2 adsorption in aluminosilicate zeolites for postcombustion carbon dioxide capture

Molecular simulation studies of separation of CO2/N2, CO2/CH4, and CH4/N2 by ZIFs

Understanding CO2 dynamics in metal-organic frameworks with open metal sites

Carbon dioxide separation with a two-dimensional polymer membrane

Systematic tuning and multifunctionalization of covalent organic polymers for enhanced carbon capture

A hybrid absorption-adsorption method to efficiently capture carbon

Accurate force field development for modeling conjugated polymers

First principles investigation of hydrogen physical adsorption on graphynes' layers

Separation of hydrogen and nitrogen gases with porous graphene membrane

A bondbond portable approach to intermolecular interactions: simulations for nmethylacetamide and carbon dioxide dimers

Modeling of energy transfer from vibrationally excited CO2 molecules: cross sections and probabilities for kinetic modeling of atmospheres, flows, and plasmas

Modeling the intermolecular interactions and characterization of the dynamics of collisional autoionization processes

On the semiclassical initial value calculation of thermal rate coefficients for the N + N2 reaction

Thermal rate coefficients in collinear versus bent transition state reactions: the N + N2 case study

An extension of the grid empowered molecular simulator to quantum reactive scattering

A nonorthogonal coordinate approach to atom-diatom parallel reactive scattering calculations

Design and implementation of a Grid application for direct calculations of reactive rates

Grid calculation tools for massive applications of collision dynamics simulations: carbon dioxide energy transfer

UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations

AMBER. A package of computer-programs for applying molecular mechanics, normal-mode analysis, molecular-dynamics and free-energy calculations to simulate the structural and energetic properties of molecules

Transferability of CO2 force fields for prediction of adsorption properties in all-silica zeolites

Force-field prediction of materials properties in metal-organic frameworks

Force-field development from electronic structure calculations with periodic boundary conditions: applications to gaseous adsorption and transport in metal-organic frameworks

Potential models for the simulation of methane adsorption on graphene: development and CCSD(T) benchmarks

Energy transfer upon collision of selectively excited CO2 molecules: state-to-state cross sections and probabilities for modeling of atmospheres and gaseous flows

An innovative synergistic grid approach to the computational study of protein aggregation mechanisms

Full dimensional quantum versus semiclassical reactivity for the bent transition state reaction N + N2

Nanostructure selectivity for molecular adsorption and separation: the case of graphyne layers

Adsorption of hydrogen molecules on carbon nanotubes using quantum chemistry and molecular dynamics

Multi-scale theoretical investigation of molecular hydrogen adsorption over graphene: coronene as a case study

Molecular simulations of CO2/N2/H2O gaseous mixture separation in graphtriyne membrane

Graphynes: indispensable nanoporous architectures in carbon flatland

A novel nanoporous graphite based on graphynes: first-principles structure and carbon dioxide preferential physisorption

Energy transfer dynamics and kinetics of elementary processes (promoted) by gas-phase CO2-N2 collisions: selectivity control by the anisotropy of the interaction

A full dimensional grid empowered simulation of the CO2 + CO2 processes

On the development of an effective model potential to describe water interaction in neutral and ionic clusters

Carbon dioxide clathrate hydrates: selective role of intermolecular interactions and action of the SDS catalyst

Beyond the lennard-jones model: a simple and accurate potential function probed by high resolution scattering data useful for molecular dynamics simulations

Carbon oxides in gas flows and earth and planetary atmospheres: state-to-state simulations of energy transfer and dissociation reactions

The molecular stirrer catalytic effect in methane ice formation

A force field for acetone: the transition from small clusters to liquid phase investigated by molecular dynamics simulations

Ion-water cluster molecular dynamics using a semiempirical intermolecular potential

Water (H2O)m or Benzene (C6H6)n Aggregates to Solvate the K + ?

Accurate analytic intermolecular potential for the simulation of Na + and K + ion hydration in liquid water

Collisional energy exchange in CO2-N2 gaseous mixtures

Ion size influence on the Ar solvation shells of M + C6F6 clusters

Atombond pairwise additive representation for intermolecular potential energy surfaces

A highlevel ab initio study of the N2 + N2 reaction channel

DL POLY: application to molecularsimulation