MD approach is based on Newton’s laws currently implemented in the various nanostructures [28–29]. According to the computational view, the MD approach is an appropriate way for estimating the atom-based compound’s behavior (physical properties) as a function of passing time. In this way, atoms interacted for a defined time, giving the complete vision for the structural evolution. Previous atoms displacement are settled by solving Newton’s equation in the general version of MD method, at which inter-atomic expressions are estimated various interactions between defined sections of materials. Technically, we used the atomic arrangement to study the iron based matrix behavior for atoms separation from O2 environment. MD simulations in current work were done using the 2021 version of the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package released by Sandia National Laboratories [30–32].

MD outputs are dramatically based on the selected potential functions [33]. To achieve the actual results, one should choose this computational parameter based on the study’s atomic properties. In MD simulations of C, O, and Fe-based structures, atoms’ interatomic potential was based on the Embedded Atom Model (EAM) and DREIDING force fields [34–36]. In DREIDING force field, Lennard-Jones (LJ) interatomic function (potential) has been applied to describe the non-bond force between defined atoms in O2-CO2 mixture system (inside MD box) [37],

\(U({r_{ij}})=4{\varepsilon _{}}\left[ {{{\left( {\frac{\sigma }{{{r_{ij}}}}} \right)}^{12}} - {{\left( {\frac{{{\sigma _{}}}}{{{r_{ij}}}}} \right)}^6}} \right]\) \({r_{ij}} \leqslant {r_c}\) (1)

where epsilon, sigma, and rij parameters represents the depth of the LJ function, the finite distance, and the distance between particles (atoms), respectively. In this equation, the cut-off constant (rc), set to 12 Å for all atom types. These LJ potential parameters for various atoms in the MD box are reported in Table 1.

Table 1

The epsilon and sigma parameters for defined LJ potential function in O2-CO2 system in vicinity of porous iron matrix [36]. Cutoff radius in various defined atoms set to 12 Å.

Element | Sigma(Å) | Epsilon(kJ/mol) |

C | 4.18 | 0.305 |

O | 3.71 | 0.415 |

Fe | 4.54 | 0.055 |

By estimating simple and angular components of inter-atomic bonds, the bonded interactions can be described. Currently, the bonded interactions in O2-CO2 system defined with a simple oscillator as equations (2) and (3) [38]:

$$E{\text{ }}={\text{ }}\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} {\text{ }}{k_r}{\left( {r - {r_0}} \right)^2}$$

2

$$E{\text{ }}={\text{ }}\raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} {\text{ }}{k_\theta }{\left( {\theta - {\theta _0}} \right)^2}$$

3

where kr, r0, kθ, and θ0 represent the harmonic oscillator constant, the atomic bond length, the angular oscillator constant, and the angle’s equilibrium value, respectively. These interaction constants for various atoms in the MD simulation box reported in Table 2.

Table 2

The kr/kθ and r0/θ0 constants for O2-CO2 system simulations' simple/angular bonded interaction [36].

Physical Parameter (Bond/Angle) | kr/kθ(kJ/mol) | r0 (Å)/θ0(degree) |

O-O Bond | 700 | 1.31 |

O-C Bond | 700 | 1.42 |

O-C-O Angle | 100 | 109.471 |

Atoms interaction in porous iron matrix defined by EAM force field. In EAM function, particles interaction described as below [35]:

$${E}_{i}={F}_{\alpha }\left(\sum _{i\ne j}{\rho }_{\beta }\left({r}_{ij}\right)\right)+\frac{1}{2}\sum _{j\ne i}{\phi }_{\alpha \beta }\left({r}_{ij}\right)$$

4

In Eq. (1), F is the embedding energy. This computational parameter is a function of electron density for each atoms (ρ). Furthermore, φ is an atomic interaction parameter and α and β are the element types of atoms i and j, respectively.

After inter-atomic force fields setting for atomic compounds inside computational box, the MD process was fulfilled. Technically, Newton’s equation is solved as the gradient of the potential function to describe the particle evolution as a function of passed time steps [39],

$${F_i}=\sum\limits_{{i \ne j}} {{F_{ij}}={m_i}\frac{{{d^2}{r_i}}}{{d{t^2}}}={m_i}\frac{{d{v_i}}}{{dt}}}$$

5

$${F_{ij}}={\text{ }} - grad{\text{ }}{V_{ij}}$$

6

in these equations, mi is particle’s mass, ri is the coordination of particle i, dt is computational time step, and Vij is the force field function. Recently, the Velocity-Verlet approach solved motion equations’ association in various MD-based packages [40–42]. In stated MD formalism, various ensembles are implemented to creating initial conditions inside box. We used NPT ensemble for the thermal/pressure equilibration phase of atomic structures detection in current work. In this ensemble, the Noose-Hoover barostat was used by the LAMMPS software (2021 version) [43–44]. Numerically, this barostat was implemented to access the equilibrium phase for 10 ns; but MD simulations were carried out until 10 ns later to simulate the atomic filtering process by using micro-canonical ensemble (NVE) [45]. Finally, this ensemble is used for mechanical tests implementing to membrane structure. So, MD simulations in our computational research consists of 2 main steps:

**Step A) Atomic purification performance of iron based membrane**

Firstly, the pristine porous iron membrane was modeled inside simulation box with 150 Å lengths in X and Y directions, 300 Å in Z direction. Also, boundary conditions in X and Y directions are defined periodically, and in the direction of Z is considered fixed one. Our atomic arrangement for this step of MD simulations depicted in Fig. 1. In this step, Nose-Hoover algorithm was then implemented to equilibrate of defined compound’s temperature at 300 K and 1 bar as initial temperature and pressure, respectively. After equilibrium phase detection, CO2 molecules, separation from O2 molecules was done with the porous iron matrix. For this atomic process reporting, physical parameters such as temperature, total energy, and CO2 molecules absorption ratio were calculated.

**Step B) Mechanical behavior of iron based membrane**

Next, the mechanical deforming process was implemented to pristine and after filtration iron based matrix by using the NVE ensemble. After the tensile test, physical parameters such as stress-strain curve, Young’s modulus, and ultimate strength were reported to describe the atomic behavior of pristine and final metallic membranes. MD simulation details reported in Table 3.

Table 3

MD simulation settings in current computational research.

Computational Parameter | Value/Setting |

Computational Box Length | 150×150×300 Å3 |

Boundary Condition | P-P-F |

Initial Temperature | 300 K |

Initial Pressure | 1 bar |

Time Step | 1 fs |

Computational Algorithm | NPT |

Temperature Damping Ratio | 10 |

Pressure Damping Ratio | 100 |

Equilibrium Time | 10 ns |

Total Simulation Time | 20 ns |

Number of Pollution Molecules | 175 |