An Innovative Approach For Molecular Simulation Of Nano-Structured

Jan 16, 2018 - The Akaike information criterion is used as a part of the proposed method to determine the optimum number of required iterations. It wa...
0 downloads 10 Views 2MB Size
Subscriber access provided by READING UNIV

Article

An Innovative Approach For Molecular Simulation Of NanoStructured Adsorption Isotherms Via Ant Colony Method Leila LotfiKatooli, and Akbar Shahsavand J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b09476 • Publication Date (Web): 16 Jan 2018 Downloaded from http://pubs.acs.org on January 16, 2018

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

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

Page 1 of 23 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

An Innovative Approach for Molecular Simulation of Nano-Structured Adsorption Isotherms via Ant Colony Method L. LotfiKatooli and A. Shahsavand1 Chemical Engineering Department, Ferdowsi university of Mashhad, Mashhad, Iran

Abstract The original Grand Canonical Monte Carlo (GCMC) molecular simulation method heavily relies on pure random steps for its various phases of displacement, removal and insertion. A new algorithm is presented in this article which employs the global optimization Ant Colony technique in a state-of-the-art algorithm to predict the adsorption isotherm of any adsorbate inside nanostructured adsorbents. Several experimentally measured isotherms are used from literature to successfully validate the proposed method. A detailed flowchart is presented to describe the entire algorithm. The Akaike information criterion is used as a part of the proposed method to determine the optimum number of required iterations. It was clearly demonstrated that our in-house ant colony molecular simulation (ACMS) method outperforms conventional GCMC method for all case studies. The simulation results also indicate that the newly proposed algorithm performs order of magnitude faster than GCMC with the same or better accuracies.

Keywords: Molecular simulation, Ant colony Algorithm, Adsorption isotherm, Grand Canonical Monte Carlo.

1

Corresponding author;

Email: [email protected],

Tel. & Fax: +98 51 38816840

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

Introduction Computerized simulation methods such as Molecular simulation and molecular dynamics have been vastly employed as very powerful tools to predict the adsorption isotherms of heterogeneous adsorbents.1 The comparison of the simulation results with the available experimental data are quite helpful for correct interpretation of existing theoretical models. The well-known Grand Canonical Monte Carlo (GCMC) simulation method which is massively used for characterization of solid adsorption surfaces entirely relies on random procedures to sample a statistical mechanical ensemble and determine the average quantities such as the equilibrium uptake and enthalpy.2 The first Monte Carlo (MC) molecular simulation was introduced by Metropolis et al.3 in 1953 at Los Alamos National Laboratory. A Markov chain was utilized to produce a series of molecular configurations that randomly moved from one arrangement to the next one using the most proper probability rules.4 Theoretical and numerical studies on characterization of adsorption materials have been immensely accelerated in recent years due to the explosive growth of computational capacity of modern computers5,6 Snur et al.7 successfully applied GCMC technique in 1992 to predict various isotherms and corresponding heats of sorption for adsorption of methane on zeolite-silicalite over a wide range of temperatures. In the same year, Siepmann and Frenkel introduced a new version of GCMC for simulation of systems containing flexible chain molecules. The proposed method was based on self-avoiding of random walk algorithm of Rosenbluth and Rosenbluth.8 In 2003 Akten et al.9 recruited molecular simulation technique for prediction of Hydrogen, Nitrogen and Carbon dioxide adsorptions on zeolite Na-4A. The interatomic potentials were computed for both singlecomponent and binary-mixtures by comparing the GCMC predictions with the existing experimental data obtained at room temperature. Furthermore, the effects of both adsorption temperature and composition of gas phase on the selectivity of binary-systems were also investigated. In 2005, Maurin et al.10 coupled the Micro-calorimetry measurements with GCMC simulation to provide a better understanding of molecular interactions between Carbon dioxide and various Faujasitic Zeoliete structures. In 2006, Kumar et al.11 also coupled the modified Feynman-Hibbs variational technique with Monte Carlo and molecular dynamics simulations for adsorption of hydrogen and deuterium in RHO zeolite from 30K to 150 K. They observed that the deuterium molecules, usually diffuse much faster than hydrogen molecules and can lead to a reverse kinetic molecular sieving phenomenon. Their experimental measurements for adsorption of hydrogen on RHO Zeolite is used in the present article to validate the accuracy of the proposed method. In 2007 Shi and Maginn12 proposed a new open system Monte Carlo procedure to overcome difficulties with insertion and deletion of molecules. The method utilizes gradual insertions and deletions of molecules through the use of a continuous coupling parameter and an adaptive bias potential. The method expanded 2 ACS Paragon Plus Environment

Page 2 of 23

Page 3 of 23 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

ensemble Monte Carlo techniques and is applied to both the grand canonical and osmotic ensembles. It was shown that the method provides correct results for the volumetric properties of the Lennard-Jones fluid and water as well as the phase behavior of the CO2-ethanol binary system. Coasne et al.13 reported a molecular simulation to predict the mechanical properties of some microporous zeolites packed with certain guest molecules in 2011. They observed that the adsorption of molecules in the micropores of the porous material increases the corresponding bulk modulus. In 2013, Dresselhaus et al.14 presented a hybrid metaheuristics optimization (HMO) approach for nonlocal optimization of molecular systems. The introduced algorithm recruited particle swarm optimization (PSO) to optimize the internal parameters of ant colony optimization (ACO) technique. In the same year, Nickmand et al.15 used Grand Canonical Monte Carlo simulation to study the adsorption of argon at 87K in wedge shaped mesopores. Various structural parameters, including mean pore size, wedge angle and wall length were changed to investigate their effects on the different aspects of the hysteresis loop. It was claimed that the obtained results have far-reaching consequences for the characterization of pore size distribution. In constructing a kernel of local isotherms, the pore size is usually assumed to be Dirac delta function (impulse). Evidently, a small deviation from a constant pore width can shift the condensation and evaporation pressures significantly and thus lead to an incorrect estimation of pore size distribution.16 Li and Zhang17 in 2014 proposed a general multi-scale simulation procedure to accurately predict the uptakes of pollution gases such as CO2, SO2, H2S, and CO in one of the most investigated porous Zeolite with organic cages of CC3 by using the sophisticated force field of vdW3 fitted by double hybrid functional (B2PLYP) with a dispersion correction (D3) separately for gas–gas and CC3- gas interactions. The fitted vdW3 was used in grand canonical Monte Carlo simulations instead of conventional Lenard-Jones force field. In the same year, Pham et al.18 investigated the adsorption isotherms of carbon dioxide and nitrogen on various low-defect siliceous zeolites in a fluoride medium at different temperatures at sub-atmospheric pressures. Grand Canonical Monte Carlo (GCMC) technique was used with various force fields to simulate the adsorption of the above constituents on the corresponding zeolites. They reported excellent agreements between the simulation results and their experimental isotherms. These measurements will be used in the upcoming sections to validate our new in-house method using ant colony optimization technique for molecular simulation purpose. Wu et al.19 also studied adsorption and separation performances of ethane/ethylene on various forms of ZIF Metal-organic framework adsorbents with different topologies, in 2014. They combined GCMC simulations with the so called ideal adsorbed solution theory (IAST). They reported that the single component isotherm of ZIF-3 with the pore size of 8.02 Å has the most adsorption capacity for both pure ethane and ethylene at extremely low pressure of around 30 KPa. 3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

In 2015, Chang et al.20 predicted the adsorption isotherms of carbon dioxide, nitrogen and water molecules on MCM-41 meso-porous molecular sieves via GCMC. They concluded that according to the calculated adsorption energy distributions, the adsorption mechanisms of gas in MCM-41 can be divided into three types of: a) surface adsorption on the pore wall, b) multilayer adsorption of the adsorbed gas molecules and c) molecular self-aggregation near the pore center. One year later and in 2016, Prasetyo et al.21 compared the Monte Carlo simulation results with their adsorption isotherms collected at 77 K and 87 K for adsorption of krypton on the highly graphitized carbon black of Carbopack F. They verified their previous findings (Diao et al., 2015)22 for both horizontal and vertical hysteresis-loops. In our previous research, a new approach was proposed to predict the adsorption isotherms by coupling the genetic algorithm with molecular simulation and presenting the so called GAMS method.23 Several adsorption isotherms were borrowed from literature to clearly demonstrate that the proposed method successfully outperforms the conventional GCMC technique. While our new GAMS method predicted the corresponding adsorption isotherms very brilliantly, however, its complexity made its application too limited for non-mathematicians. In this article a similar approach is presented which uses Ant Colony optimization technique as a platform to construct and validate a much simpler technique named as "Ant Colony Molecular Simulation (ACMS)" method to successfully predict the adsorption isotherms of various crystalline (BEA, CHA, RHO) Zeolites adsorbents. AS the GAMS, this new procedure also uses a smart search algorithm instead of pure random search employed in conventional Grand Canonical Monte Carlo (GCMC) molecular simulation technique. The proposed algorithm is not only much simpler to understand and implement than our previous work (the so called GAMS) but it requires exceedingly less computational time to reach nearly the same accuracy. The first initial sections of the following passages will provide a brief overview of conventional GCMC and Ant Colony methods. Afterwards, a detailed description of our newly presented ACMS in-house procedure will be introduced. Finally, several superior performances of the ACMS method will be compared with the conventional GCMC technique by resorting to three experimental adsorption isotherm measurements borrowed from literature11,18 Methods A concise overview of conventional GCMC technique The purely random procedure of GCMC enables simulation in the grand canonical ensemble (at constant chemical potential μ, volume V, and temperature T). GCMC method employs three types of moves: a) the displacement of a randomly chosen atom, b) the removal or destruction of an atom selected at random, and c) the insertion of a new atom with a random location.24 In a GCMC step, the attempted displacement

4 ACS Paragon Plus Environment

Page 4 of 23

Page 5 of 23 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

of a single atom is followed by either an attempted insertion or removal of a single atom, chosen with equal probability. The following stepwise search algorithm is exploited to find the (sub) optimal positions of the adsorbed molecules inside the simulation box based on their corresponding adsorbate-adsorbate and adsorbent-adsorbate potential energy levels. 1. Specify the coordinates of both adsorbent and (initial) adsorbate atoms inside the simulation box. 2. Select one of the three original steps of GCMC using a random number (R≤⅓: displacement, ⅓⅔: insertion). 3. If the total potential energy of the system in decreased or increased below a certain value known as Boltzman factor (β), the change is accepted, otherwise, the change is discarded. 4. The above two steps are repeated until a certain number of predefined iterations is reached.4,5

Brief description of conventional Ant colony Algorithm The Ant Colony optimization algorithm is a metaheuristic class of optimization technique introduced by Dorigo and his co-workers in the early nineties.25 Natural ants use special kinds of pheromones to pinpoint the food sources around their living environment. This pheromone trail will aid other members to communicate and follow the path. Therefore, short paths with more visitation which leads to higher pheromone concentration are more attractive and ants choose such paths more frequently. Over the time and because of the evaporation process the pheromone level of less visited long paths will be reduced more profoundly and finally should be disappeared.26-28 The following ant colony optimization algorithm mimics the above natural behavior and searches inside a pre-specified domain to find the global optimum solution.

1. Select both grid type and grid size for the entire input domain. 2. Initialize pheromone level of all paths τ0 (including boundaries). 3. Position a pre-specified number of ants (na) inside the grid, at random. 4. For each ant (i=1,…,na ): a) Compute both individual and cumulative probabilities for all paths starting from ant position based on their corresponding pheromone levels. b) Move ant along a connection (path) inside the grid by comparing a random number with the calculated cumulative probability. c) Calculate the cost function value of the ant at the new position and update pheromone levels of both forward and reverse paths accordingly. 5. Apply the evaporation losses for all connections (paths) 6. Repeat from step 4 until the termination criterion (maximum number of iterations) is reached. 5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 6 of 23

In a later variant, Thomas Stützle et al.29 introduced the MAX-MIN Ant System (MMAS), in 2000 which limits the total pheromone level in every trip or sub-union to avoid local convergence.

Use of Ant Colony for Molecular Simulation The present article introduces a new technique which employs the ant colony (AC) algorithm as an optimization tool for molecular simulation (MS). The proposed algorithm of Ant Colony Molecular Simulation (ACMS) can be used successfully to predict any adsorption isotherm for structured adsorbent materials such as Zeolites. Figure 1 provides detailed description of our newly proposed flowchart for ACMS algorithm. As before (in GAMS), any conventional software (e.g. MERCURY) can be used to extract the exact locations of the adsorbent atoms inside the simulation box. A brief comparison of Figure 1 with our previously presented flowchart of GAMS method23, clearly reveals that the new ACMS method is not only much simpler to understand and develop, but also it will performs much faster as will be discussed later. As before, the following cost function should be optimized to find the proper number of molecules inside the simulation box at the prevailing temperature and pressure: 𝑌(𝑁) = −𝑁𝜇∗ + 𝐾𝐵 𝑇𝑙𝑛 (𝑁!) + 𝐸(𝑟 𝑁 )

(1)

𝜇∗ (𝑇, 𝑃, 𝑁) = 𝜇𝑒𝑥 (𝑇, 𝑃) + 𝐾𝐵 𝑇𝑙𝑛

(2)

Where 𝐾𝐵 is the Boltzmann constant and 𝜇𝑒𝑥 is the excess chemical potential of the adsorbate molecules which should be computed via any proper equation of state (EOS).4 Instead of using equation (2), the chemical potential value at any given temperature and pressure is computed by resorting to 𝜇 =

1 𝛽

ln(Λ3 𝛽φ𝑃), where β is the Boltzmann factor defined as 1/ (𝐾𝐵 𝑇), Λ is

the de Broglie thermal wavelength (ℎ/√2𝜋𝑚𝐾𝐵 𝑇), h is the Planck constant, m is the mass of a gas particle and φ is the real gas fugacity coefficient. In equation (1), the value of 𝐸(𝑟 𝑁 ) can be obtained by resorting to any appropriate potential function such as Lennard-Jones (LJ). It is one of the most popular potential energy functions which provides the potential energy of each individual adsorbate molecule with respect to other adsorbate and adsorbent atoms: 𝜎𝑖𝑗

𝐸𝑖𝑗 = 4𝜀𝑖𝑗 ((

𝑟𝑖𝑗

12

)

𝜎𝑖𝑗

−(

𝑟𝑖𝑗

6

) )

(3)

Where (𝜎𝑖𝑗 = (𝜎𝑖 + 𝜎𝑗 )⁄2) and (𝜀𝑖𝑗 = √𝜀𝑖 𝜀𝑗 ) are the LJ potential parameters (molecular distances at collision and corresponding potential energy release) and 𝑟𝑖𝑗 is the interatomic distance between the ith and jth atoms or molecules.30

6 ACS Paragon Plus Environment

Page 7 of 23 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 start

Data Entry Section: Input the following items a) Number of unit cells for each spatial direction, b) Coordinates of adsorbent atoms in the simulation box, c) All required physical and thermodynamical properties of adsorbate(s) and adsorbent(s) at prevailing T and minimum pressure increment (P) d)Internal parameters of ACO method (number of grids, promotion, punishment and evaporation rates

Initialize the pheromone level (T0) in all paths and zeroise all pheromone levels of simulation box boundaries Specify the smallest pertinent assumed number of ants (N)for the given T& P and position them at random inside the previously constructed spatial grid. Ensure that the primary positions of ants lead to negative total energies to enhance faster convergence. Choose initial positions of the adsorbate molecules as the best position. Compute potential energy values for all combinations of atoms and molecules using appropriate model.

Calculate the value of objective function (e.g. Y=-Nμ*+KTln(N!)+E(N)=f(T,P,N,E))

P = Pmin jflag= 0, j=0

Itr = 0 ,

Itr=Itr+1 ACO steps: For each Ant (i=1,...,N) a) Calculate individual and cumulative probabilities for all corresponding paths. b) Select new path using a random number and cumulative probability of all paths. c) Move ant to its new position along the chosen path. d) Update pheromone levels (punish or promote) of chosen path and its corresponding reverse connection, appropriately. a) Apply evaporation rate to update all pheromone levels of connections and their reverse paths. b)Recompute the value of objective function (Y) c) Select the current positions af all adsorbate molecules (Ants) as the best one (Ybest) if Y