A method for calculating the contact angle of carbon dioxide based on two-scale molecular simulation

Abstract

The present invention discloses a method for calculating the contact angle of carbon dioxide based on two-scale molecular simulation, belonging to the fields of molecular simulation and interface science, and comprising the following steps: Step 1, performing all-atom molecular dynamics calculation; Step 2, constructing a coarse-grained molecular dynamics potential energy model, performing coarse-grained simulation, and predicting the wetting behavior of the CO2 interface; Step 3, constructing a molecular dynamics-Monte Carlo algorithm by combining molecular dynamics and Monte Carlo methods for two-scale molecular simulation. The present invention can efficiently and accurately calculate the contact angle of CO2 on the solid surface, and further effectively simulate the wetting behavior of CO2 on the solid surface.

Description

Technical Field

[0001] The present invention belongs to the field of molecular simulation and interfacial science, and particularly relates to a method for calculating the contact angle of carbon dioxide based on two-scale molecular simulation. Background Art

[0002] The wetting behavior of carbon dioxide (CO2) on solid surfaces plays a crucial role in applications such as enhanced oil recovery (EOR), carbon dioxide geological storage (CCUS), and membrane separation. The contact angle is a core parameter for measuring the wetting performance of fluids on solid surfaces, but there are still certain limitations in current calculation methods: all-atom molecular dynamics (AA-MD) simulations can accurately describe interfacial properties, but the computational cost is high and it is difficult to scale up to large-scale systems; the continuum model (CM) cannot accurately capture interfacial phenomena at the microscale, affecting simulation accuracy; traditional Monte Carlo methods (MC) do not fully combine molecular trajectory information, resulting in limited computational efficiency and accuracy. Summary of the Invention

[0003] To solve the above problems, the present invention proposes a method for calculating the contact angle of carbon dioxide based on two-scale molecular simulation. This method combines molecular dynamics (MD) and Monte Carlo method (MC) to construct a molecular dynamics-Monte Carlo (MDMC) algorithm for two-scale molecular simulation. By constructing a state transition matrix, this algorithm can efficiently and accurately calculate the contact angle of CO2 on solid surfaces, and thus effectively simulate the wetting behavior of CO2 on solid surfaces.

[0004] The technical solution of the present invention is as follows:

[0005] A method for calculating the contact angle of carbon dioxide based on two-scale molecular simulation, comprising the following steps:

[0006] Step 1: Perform all-atom molecular dynamics calculations;

[0007] Step 2: Construct a coarse-grained molecular dynamics potential energy model and perform coarse-grained simulations to predict the wetting behavior of the CO2 interface;

[0008] Step 3: Combine molecular dynamics and Monte Carlo methods to construct a molecular dynamics-Monte Carlo algorithm for two-scale molecular simulation.

[0009] Further, the specific process of Step 1 is as follows:

[0010] Step 1.1: Select the CO2 molecular force field of the phase equilibrium transferable potential model and construct a CO2 / solid interface system;

[0011] Step 1.2: Perform all-atom scale molecular dynamics simulations to obtain interfacial properties; the interfacial properties include the adsorption energy of CO2 molecules and the contact angle of CO2 molecules on the solid surface.

[0012] Furthermore, the specific process of step 1.1 is as follows:

[0013] Step 1.1.1: The selected CO2 molecular force field includes Lennard-Jones potential energy parameters and related charge distribution information; the Lennard-Jones potential energy parameters include the well depth and the distance between particles.

[0014] Step 1.1.2: Use the Lennard-Jones potential energy parameters to construct the CO2 intermolecular interaction potential energy function and calculate the intermolecular interaction potential energy:

[0015] (1);

[0016] where is the intermolecular interaction potential energy; is the distance between two particles; is the well depth; is the distance between particles;

[0017] Step 1.1.3: Determine the crystal structure of the solid surface and establish an atomic model.

[0018] Step 1.1.4: Use a large-scale atomic / molecular massively parallel simulator to perform all-atom scale structural modeling and energy minimization operations. After this operation is completed, the CO2 / solid interface system is constructed; the large-scale atomic / molecular massively parallel simulator synchronously generates a trajectory file.

[0019] Furthermore, in step 1.2, the calculation formula for the adsorption energy of CO2 molecules is:

[0020] (2);

[0021] where is the adsorption energy of CO2 molecules; is the total energy of the entire simulation system; is the independent energy of CO2 molecules; is the energy of the solid surface alone;

[0022] The calculation formula for the contact angle of CO2 molecules on the solid surface is:

[0023] (3);

[0024] where is the contact angle of CO2 molecules on the solid surface; is the solid-gas interfacial tension; is the solid-liquid interfacial tension; is the liquid-gas interfacial tension.

[0025] Furthermore, the specific process of step 2 is as follows:

[0026] Step 2.1: Use the Boltzmann inversion method for potential energy fitting, construct a coarse-grained molecular dynamics potential energy model and optimize it; the formula of the coarse-grained molecular dynamics potential energy model is:

[0027] (4);

[0028] where is the effective interaction potential energy function between coarse particles; is the Boltzmann constant; is the system temperature; is the radial distribution function obtained from all-atom simulation;

[0029] Optimize the coarse-grained molecular dynamics potential energy model through the force matching method; during the optimization process, construct the following objective function to measure the difference in particle forces between the all-atom model and the coarse model:

[0030] (5);

[0031] where is the objective function; is the real force calculated for the th particle in the all-atom model simulation; is the predicted force calculated for the th particle under the current coarse model parameters ; is the total number of simulated particles;

[0032] Step 2.2: Conduct a coarse-grained simulation to predict the wetting behavior of CO2 molecules at the interface; the wetting behavior of CO2 molecules at the interface is specifically the contact angle of CO2 molecules on the solid surface; calculate the contact angle formed by CO2 molecules on the solid surface according to formula (3);

[0033] Pre-optimize the initial spatial distribution of CO2 molecules using the Monte Carlo method. The specific process is as follows:

[0034] Under the condition of fixing the effective interaction potential energy function between coarse particles, randomly perturb the positions of CO2 molecules through an acceptance-rejection strategy. Combining with the Metropolis criterion, the probability of accepting the new state is:

[0035] (6);

[0036] where is the probability of accepting the new state; is the exponential function with as the base; , are the current position of the particle and the candidate position after perturbation, respectively.

[0037] Furthermore, the specific process of step 3 is as follows:

[0038] Step 3.1: Molecular dynamics simulation data collection and preprocessing; The data collection process is as follows: By parsing the trajectory file generated in step 1.1.4, the three-dimensional coordinate information of all particles and the simulation boundary information at each time step in the molecular dynamics simulation are extracted; The preprocessing process is as follows: First, eliminate the non-numeric tags in the three-dimensional coordinate information, then perform format conversion, and finally store it in a three-dimensional data structure;

[0039] Step 3.2: Discretized grid construction and state division, perform molecular dynamics simulation trajectory state statistics, and construct the original state statistical matrix;

[0040] Step 3.3: Construct the coordinate expansion matrix;

[0041] Step 3.4: Construct the state transition probability matrix to record the frequency of particles transitioning from one state to another;

[0042] Step 3.5: Use the Monte Carlo method to reconstruct the microscopic trajectory of CO2 molecules to simulate its multi-step transition process in the discrete state space.

[0043] Furthermore, the specific process of step 3.2 is as follows:

[0044] Step 3.2.1: In order to discretize the motion trajectory of molecules in three-dimensional space, first, a unified grid system needs to be constructed and all particle states need to be numbered; First, traverse the z-direction coordinate data of all particles in the trajectory file, and use the minimum value function and the maximum value function to determine the lowest position and the highest position that the particles have appeared in throughout the simulation; Then perform a unified translation process on the z-direction coordinate data of all particles, that is, subtract from each z-direction coordinate, so as to translate the entire coordinate interval to a positive number interval starting from zero ;

[0045] Step 3.2.2: Discretely divide the entire three-dimensional space according to the simulation box size and the set grid size. After division, several states are obtained, and the state matrix is initialized accordingly; Specifically, use the grid division method to convert the CO2 molecular trajectory into discrete states; At the same time, calculate the probability of a particle transitioning from one state to another, which is used to guide the statistics of state transition frequencies. The formula is:

[0046] (7);

[0047] Among them, is the probability that the particle jumps from state to state ; is the number of times the particle jumps from state to state ; is the total number of states;

[0048] Step 3.2.3: Read the molecular dynamics simulation trajectory data, traverse the spatial coordinates of each particle at each time step, and sequentially extract its position data in the x, y, and z directions; perform periodic boundary condition processing and normalized coordinate transformation on the atomic coordinates. Then, calculate the spatial state number where the particle is located according to the grid size in each direction, and store this number in the molecular trajectory array for marking. At the same time, accumulate the number of occurrences in the statistical matrix under the corresponding state index. Finally, the original state statistical matrix is obtained. The specific definition of the original state statistical matrix is , indicating the cumulative number of occurrences of the discrete state numbered in the molecular dynamics simulation;

[0049] The formula for performing periodic boundary condition processing on atomic coordinates is:

[0050] (8);

[0051] Among them, , , respectively represent the new atomic coordinates in the x, y, and z directions after processing; , , are the initial atomic coordinates in the x, y, and z directions respectively; , , are the lengths of the simulation box in the x, y, and z directions respectively; represents the rounding operation on floating-point numbers.

[0052] Furthermore, the specific process of step 3.3 is as follows: Add three columns to the original state statistical matrix to store the spatial center coordinates of each state, and the expanded coordinate-expanded matrix is denoted as ; Among them, , , are the x, y, and z coordinates of the grid cell corresponding to the discrete state numbered respectively;

[0053] Perform the coordinate reduction operation on each state number in sequence. By using the set boundary size of the simulation box and the grid cell size, call the state-coordinate mapping function to map the state number to the geometric center coordinates in the three-dimensional grid, and store the corresponding coordinate values in the x, y, and z directions into the coordinate extension matrix;

[0054] The calculation formula of the state-coordinate mapping function is:

[0055] (9);

[0056] Among them, is the state number of the particle; , , are the discrete indices of the current grid cell where the particle is located in the x, y, and z directions respectively; , are the total number of cells in the x and y directions of the grid respectively.

[0057] Furthermore, the specific process of step 3.4 is as follows:

[0058] Step 3.4.1: Traverse the state changes of each atom between consecutive time steps, count the number of transitions from state to state , and initialize the sparse transition matrix according to the total number of states divided in step 3.2.2;

[0059] (10);

[0060] Among them, represents the total number of transitions from state to state ; is the total number of time steps; represents the state number of the particle at time step ; is the state transition time interval; is the Kronecker function;

[0061] Step 3.4.2: Traverse the state evolution trajectories of all particles at each time step, extract adjacent two state numbers at a fixed time interval, and use them as the starting state and target state of the transition respectively. Calculate the probability of each particle transitioning from one state to another according to formula (5); every time it is recognized that a particle transitions from state to state , accumulate the value at the corresponding position in the matrix, indicating that this transition occurs once; this process is repeated within the range of all particles and the entire time series, and finally a complete state transition probability matrix is formed;

[0062] Step 3.4.3, Normalize the state transition probability matrix to describe state transitions in the form of probabilities. The specific process is as follows: Normalize the number of transitions for each row to obtain the normalized state transition probability matrix ;

[0063] Step 3.4.4, Calculate the steady-state probability distribution of the normalized state transition probability matrix:

[0064] (11);

[0065] where, is the steady-state probability distribution; represents the steady-state probability distribution in each state.

[0066] Furthermore, the specific process of Step 3.5 is as follows:

[0067] Step 3.5.1, Initialize the spatial state of each particle at the initial time step and set a random number within the interval [0, 1] as the transition criterion;

[0068] Step 3.5.2, Progress each particle step by step in the time series. Using the normalized state transition probability matrix in Step 3.4.3, calculate the cumulative distribution function corresponding to the current state row; by comparing the random number with the cumulative distribution function, determine whether the particle makes a transition and its target state number; obtain the transition probability corresponding to the current state from the normalized state transition probability matrix. If the transition probability is zero, retain the previous state to avoid trajectory interruption;

[0069] Step 3.5.3, Use the Metropolis-Hastings algorithm for probability sampling to simulate the transition process of CO2 molecules; the Metropolis-Hastings algorithm determines the state transition of molecules by sampling probabilities, and the formula is:

[0070] (12);

[0071] where, is the probability of accepting the new state; , are the transition probabilities of the new and old states respectively;

[0072] Calculate the residence time of CO2 molecules on the solid surface by statistically analyzing the state transitions of particles in the time series:

[0073] (13);

[0074] where, is the total number of time steps during which the particle persists in a certain specific state; represents the time step and is the state number of the particle at that time step; is the initial state number of the particle; is the Dirac function;

[0075] Calculate the transition rate:

[0076] (14);

[0077] wherein, is the transition rate;

[0078] Finally, based on the transition rate, the dynamic behavior of CO2 molecules and their interaction with the solid surface are estimated.

[0079] The beneficial technical effects brought by the present invention: The present invention constructs a state transition matrix based on the MDMC algorithm, and uses MD trajectory data to construct a Markov transition matrix to ensure the accuracy of the transition probability. A state merging and coarse-graining strategy is proposed, which merges similar states in the transition matrix, reduces the number of states, improves the calculation efficiency, and at the same time preserves the key information of the system. The periodic boundary condition (PBC) is optimized to ensure that the movement of CO2 molecules in the finite simulation region conforms to the real physical environment. Calculate the CO2 contact angle based on the state transition probability and compare and optimize it with the direct MD simulation results. The MDMC algorithm of the present invention has good migration performance and is applicable to oil and gas exploitation, carbon dioxide geological storage (CCUS), ion transport in protein channels, reservoir stimulation, and research on microscopic wetting behavior, etc. Description of the Drawings

[0080] Figure 1 is the flow chart of the method for calculating the carbon dioxide contact angle based on two-scale molecular simulation of the present invention.

[0081] Figure 2 is a schematic diagram of simulating the most basic state transition process through a simplified "toy model" in the embodiment of the present invention.

[0082] Figure 3 is a schematic diagram of an original Markov transition matrix containing 100 microscopic states in the embodiment of the present invention.

[0083] Figure 4 is a schematic diagram of the evolution of the original state trajectory of CO2 molecules over time steps in the microscopic state space in the embodiment of the present invention.

[0084] Figure 5 is the 10×10 state transition probability matrix after coarse-graining treatment in the embodiment of the present invention.

[0085] Figure 6It is the time trajectory diagram of the coarse-grained state in the embodiment of the present invention.

[0086] Figure 7 It is the schematic diagram of the contact angle configuration in the wetting behavior of CO2 molecules on the quartz surface in the embodiment of the present invention. Detailed implementation manners

[0087] The present invention will be further described in detail below in conjunction with the accompanying drawings and specific implementation manners:

[0088] As Figure 1 shown, a microscopic calculation method for the contact angle of carbon dioxide based on two-scale molecular simulation includes the following steps:

[0089] Step 1: Perform all-atom molecular dynamics (AA-MD) calculation; the specific process of AA-MD calculation is:

[0090] Step 1.1: Select the CO2 molecular force field of the phase equilibrium transferable potential model (TraPPE, Transferable Potentials for Phase Equilibria) to construct the CO2 / solid interface system. The specific process is:

[0091] Step 1.1.1: The selected CO2 molecular force field mainly includes the Lennard-Jones potential energy parameters (such as the potential well depth and the inter-particle distance ) in the TraPPE potential function, as well as the relevant charge distribution information, which is used to accurately simulate the intermolecular interactions and thermodynamic behaviors of CO2 molecules. Among them, the potential well depth and the inter-particle distance need to be obtained through experimental fitting or quantum chemical calculation and are preset in the force field library of the phase equilibrium transferable potential model for subsequent simulation calls;

[0092] The present invention specifically includes the Lennard-Jones potential energy parameters of the carbon atoms and oxygen atoms constituting the CO2 molecule and the relevant charge distribution information; among them, the inter-particle distance of the carbon atoms is 2.80 Å, and the potential well depth / k6 = 27.0 K for the interaction between carbon atoms, and the charge of the carbon atoms is +0.70e; the inter-particle distance of the oxygen atoms is 3.05 Å, and the potential well depth / k6 = 79.0 K for the oxygen atoms, and the charge of the oxygen atoms is –0.35e. The above parameters are used for subsequent potential energy function definition and inter-atomic interaction modeling in the interface system; k6 is a constant characterizing the interaction strength.

[0093] Step 1.1.2: Construct the potential energy function of the intermolecular interaction of CO2 molecules using the Lennard-Jones potential energy parameters, and calculate the potential energy of the interaction between particles:

[0094] (1);

[0095] Wherein, is the potential energy of the interaction between particles; is the distance between two particles; is the depth of the potential well, indicating the strength of the intermolecular interaction; is the distance between particles, specifically referring to the distance when the potential energy between two particles is zero, usually representing the molecular diameter;

[0096] Step 1.1.3: Determine the crystal structure of the solid surface and establish an atomic model; specifically, α-quartz (SiO2) is used for the solid crystal surface of the present invention. Its crystal structure is a hexagonal lattice, the space group is P3121, the first lattice constant a = 4.913 Å, and the second lattice constant c = 5.405 Å. Use this structure to construct a slab structure on the (001) plane, with a thickness of about 10 Å, and place CO2 molecules on its surface for wetting behavior simulation. The wetting behavior simulation specifically includes interface contact and adsorption simulation. The slab structure is the established atomic model, which contains Si and O atoms. Some atoms are fixed to simulate a rigid surface, and the rest participate in energy minimization and kinetic evolution;

[0097] Step 1.1.4: Use the Large-scale Atomic / Molecular Massively Parallel Simulator LAMMPS to perform all-atom scale structural modeling and energy minimization operations. After this operation, the CO2 / solid interface system is constructed; LAMMPS synchronously generates a trajectory file.

[0098] Step 1.2: Perform all-atom scale molecular dynamics (NVT / NPT) simulations to obtain interface properties. The interface properties include the adsorption energy of CO2 molecules and the contact angle of CO2 molecules on the solid surface; set the simulation temperature, pressure, and boundary conditions to ensure the stability of the simulation system composed of CO2 molecules, solid surface atoms, and the boundary conditions surrounding them. Run all-atom scale molecular dynamics simulations, balance the simulation system, and collect molecular dynamics simulation trajectory data.

[0099] The calculation formula for the adsorption energy of CO2 molecules is:

[0100] (2);

[0101] Wherein, is the adsorption energy of CO2 molecules; is the total energy of the entire simulation system; is the independent energy of the CO2 molecule; is the energy of the solid surface alone;

[0102] The calculation formula for the contact angle of the CO2 molecule on the solid surface is:

[0103] (3);

[0104] where, is the contact angle of the CO2 molecule on the solid surface; is the solid–gas interfacial tension; is the solid–liquid interfacial tension; is the liquid–gas interfacial tension;

[0105] Step 2: Construct a coarse-grained molecular dynamics potential model (CG-MD potential model) and perform coarse-grained simulation to predict the CO2 interfacial wetting behavior; the specific process is as follows:

[0106] Step 2.1: Use the Boltzmann Inversion (BI) method for potential energy fitting, construct a coarse-grained molecular dynamics potential model and optimize it. The specific process is as follows:

[0107] The coarse-grained potential energy function can be deduced from the radial distribution function obtained by all-atom simulation; the formula for the coarse-grained molecular dynamics potential model is:

[0108] (4);

[0109] where, is the effective interaction potential energy function between coarse particles; is the Boltzmann constant; is the system temperature; is the radial distribution function obtained by all-atom simulation, specifically representing the relative occurrence probability density when the distance between coarse particles is ;

[0110] Further optimize the above initially obtained CG-MD potential model by the force matching method to improve the physical consistency of the coarse-grained simulation in terms of force. Specifically: construct the following objective function to measure the difference in particle forces between the all-atom (AA) model and the coarse-grained (CG) model:

[0111] (5);

[0112] where, is the objective function; is the true force calculated for the -th particle in the all-atom model simulation; is the predicted force calculated for the -th particle under the current coarse-grained model parameters ; are the parameters to be optimized in the coarse-grained model, which may be the , and other parameters or potential function point values in formula (1) depending on the specific form of the potential energy function; is the total number of simulated particles;

[0113] Step 2.2. Conduct coarse-grained simulation to predict the wetting behavior of CO2 molecules at the interface; the wetting behavior of CO2 molecules at the interface specifically refers to the contact angle of CO2 molecules on the solid surface.

[0114] Set the system parameters of the CG-MD potential model, including the time step, temperature control method, etc., to ensure that the simulation system reaches a thermodynamically stable state. After the system energy converges, calculate the contact angle formed by CO2 molecules on the solid surface according to formula (4).

[0115] In Step 2.2, to further improve the balance and stability of the CG-MD potential model simulation system, the Monte Carlo method is used to pre-optimize the initial spatial distribution of CO2 molecules. The pre-optimization process is as follows:

[0116] Under the condition of fixing the effective interaction potential function between coarse particles, randomly perturb the positions of CO2 molecules through an acceptance-rejection strategy. Combining with the Metropolis criterion, the probability of accepting the new state is given by the following formula:

[0117] (6);

[0118] where is the probability of accepting the new state; is the exponential function with as the base; , are the current position and the candidate position after perturbation of the particle, respectively;

[0119] Through multiple iterative samplings, the rapid relaxation of particles in the CG potential field is realized, making them tend to the equilibrium configuration with the lowest free energy, thereby ensuring the stability and convergence of the subsequent coarse-grained kinetic simulation process. This strategy effectively avoids the non-physical particle accumulation or density fluctuation caused by the initialization of the CG-MD potential model, and is particularly suitable for simulating the wetting behavior of CO2 molecules in complex interface systems such as in porous structures or on solid surfaces.

[0120] Step 3. Use the MDMC algorithm to discretize and statistically sample the motion trajectory of CO2 molecules on the solid surface. The main process includes trajectory reading and state encoding, state transition relationship construction, transition matrix normalization and steady-state analysis, and Monte Carlo method for trajectory reconstruction. The above process not only improves the simulation efficiency, but also retains the microscopic state transfer characteristics of the system, providing a data basis for subsequent contact angle calculation and interface property evaluation. Combining molecular dynamics (MD) and Monte Carlo method (MC) to construct a molecular dynamics-Monte Carlo (MDMC) algorithm for two-scale molecular simulation, which is used to calculate the contact angle of CO2 on the solid surface; the molecular dynamics-Monte Carlo (MDMC) algorithm mainly includes MD trajectory reading, MD trajectory state statistics, state transition matrix construction, Monte Carlo simulation, and periodic boundary condition (PBC) processing. The specific process is:

[0121] Step 3.1, molecular dynamics (MD) simulation data acquisition and preprocessing;

[0122] In this step, the three-dimensional coordinate information and simulation boundary information of all particles in each time step of the molecular dynamics simulation are extracted by parsing the trajectory file (.lammpstrj format) generated by LAMMPS. In the specific implementation process, the program reads the content of the trajectory file in a line-by-line scanning manner. After identifying the iconic keywords (such as "ITEM:TIMESTEP", "ITEM: NUMBER OF ATOMS", "ITEM: BOX BOUNDS" and "ITEM: ATOMS"), the corresponding time step number, total number of particles, boundary data of the simulation box in the x, y, and z directions, and the spatial coordinate value of each particle in the current time step are extracted respectively. The preprocessing of the particle coordinates is performed by eliminating non-numeric labels and converting the format, and storing it in a three-dimensional data structure, that is, the molecular trajectory matrix is constructed in the order of "particle number, coordinate axis number, time step number". The preprocessing process ends after the entire trajectory file is traversed, completing the numerical reconstruction of the time-evolving behavior of CO2 molecules on the solid surface, laying a data foundation for subsequent state division, state transition matrix construction and Monte Carlo simulation.

[0123] Step 3.2: Discrete grid construction and state division, perform molecular dynamics simulation trajectory state statistics, and construct the original state statistical matrix; the specific process is:

[0124] Step 3.2.1: To discretize the molecular motion trajectories in three-dimensional space, a unified grid system needs to be constructed first and all particle states need to be numbered. For the z-direction coordinates (i.e., the positions of particles in the vertical direction), it is necessary to ensure that all particle coordinates are within a unified positive value interval, which is convenient for subsequent division into discrete states. For this purpose, first traverse the z-direction coordinate data of all particles in all trajectory files, and use the minimum value function and the maximum value function to determine the lowest position and the highest position that the particles have appeared in throughout the simulation process. Then perform a unified translation process on the z-direction coordinate data of all particles, that is, subtract from each z-direction coordinate, so as to translate the entire coordinate interval to a positive number interval starting from zero .

[0125] Step 3.2.2: Discretely divide the entire three-dimensional space according to the simulation box size and the set grid size. After division, several states are obtained, and the state matrix is initialized accordingly; specifically, the grid division method is used to convert the CO2 molecular trajectories into discrete states, and the set grid sizes are , , , , , , which are the lengths of the simulation box in the x, y, and z directions respectively. At the same time, calculate the probability of a particle jumping from one state to another state, which is used to guide the statistics of the state transition frequency. The formula is:

[0126] (7);

[0127] where, is the probability that a particle jumps from state to state ; is the number of times a particle jumps from state to state ; is the total number of states.

[0128] Step 3.2.3: Using the auxiliary function Coordinate_to_State in molecular dynamics, map the continuous coordinates of each atom at each time step to discrete state numbers, and accumulate the number of atoms in each state, and store them in the MD_state file. The specific process is as follows: Read the molecular dynamics simulation trajectory data, traverse the spatial coordinates of each particle at each time step, and sequentially extract its position data in the x, y, and z directions; perform periodic boundary condition (PBC) processing and normalized coordinate transformation on the atomic coordinates. Then, calculate the spatial state number where the particle is located according to the grid size in each direction, and store this number in the molecular trajectory array for marking. At the same time, accumulate the number of occurrences in the statistical matrix under the corresponding state index. Finally, the original state statistical matrix is obtained. This matrix is the statistical result of the access frequency of all discrete states during the entire molecular dynamics simulation, and is specifically used to reflect the spatial distribution density of each discrete state. The specific definition is , indicating the cumulative number of occurrences of the discrete state in the molecular dynamics simulation.

[0129] The formula for processing atomic coordinates using periodic boundary conditions (PBC) is:

[0130] (8);

[0131] where , , respectively represent the new atomic coordinates in the x, y, and z directions after processing; , , are the initial atomic coordinates in the x, y, and z directions respectively; , , are the lengths of the simulation box in the x, y, and z directions respectively; represents rounding the floating-point number;

[0132] Through this method, the spatial distribution state encoding and frequency statistics of particles in the entire trajectory can be efficiently completed, providing an accurate state basis for subsequent construction of the state transition matrix and Monte Carlo path sampling.

[0133] Step 3.3: Construct a coordinate expansion matrix; the specific process is as follows:

[0134] Add three new columns to the original state statistical matrix to store the spatial center coordinates of each state. After expansion, the coordinate expansion matrix is obtained, usually denoted as ; where , , They are the coordinates in the x, y, and z directions of the grid cells corresponding to the discrete states numbered respectively. The coordinate expansion matrix adds three columns of spatial position parameters on the basis of the original statistical frequency, which is used for subsequent density distribution analysis, state visualization, and coarse-graining processing.

[0135] For the convenience of visual analysis of the spatial distribution of particles and subsequent density calculation, it is necessary to reverse-convert the three-dimensional discrete state numbers obtained in step 3.2 into the corresponding spatial coordinate values. The specific process of obtaining the spatial center coordinates of each state is as follows:

[0136] Perform the coordinate reduction operation on each state number in turn. Through the set boundary size of the simulation box and the size of the grid cells, call the state-coordinate mapping function to map the state number to its geometric center coordinates in the three-dimensional grid, and store the corresponding coordinate values in the x, y, and z directions into the coordinate expansion matrix.

[0137] The calculation formula of the state-coordinate mapping function is:

[0138] (9);

[0139] where is the state number of the particle; , , are the discrete indices of the current grid cell where the particle is located in the x, y, and z directions respectively; , are the total number of grid cells in the x and y directions respectively.

[0140] Before calculating the state, call the function PBC to handle the coordinate jump of atoms caused by the periodic boundary condition to ensure mapping to the correct grid. This processing process realizes the correspondence between the state number and the real space position, providing the necessary data support for subsequent analysis steps such as drawing distribution diagrams, density contour diagrams, and state trajectory diagrams.

[0141] Step 3.4: To characterize the transition behavior between different microscopic states of CO2 molecules during the molecular dynamics simulation, construct a state transition probability matrix to record the frequency of particles transitioning from one state to another. The construction process of the state transition probability matrix is as follows:

[0142] Step 3.4.1: First, traverse the state changes of each atom between consecutive time steps, and count the number of transitions from state to state The number of transitions is used to initialize the sparse transition matrix Trans_Matrix according to the total number of states divided in Step 3.2.2. Calculate the number of probability transitions according to the following formula;

[0143] (10);

[0144] where, represents the total number of transitions from state to state ; is the total number of time steps; represents the state number of the particle at time step ; is the state transition time interval; is the Kronecker function, which is 1 when the condition is satisfied and 0 otherwise.

[0145] For each particle, record its state number state1 at time step and its number state2 at time step . Consider this as an act of transitioning from state1 to state2, and increment the value of the corresponding element in the state transition matrix by one.

[0146] Step 3.4.2. Subsequently, traverse the state evolution trajectories of all particles at each time step, extract the adjacent two state numbers at a fixed time interval (i.e., the sampling step size) as the starting state and the target state of the transition respectively, and calculate the probability of each particle transitioning from one state to another state according to formula (7). Every time it is recognized that a particle transitions from state to state , increment the value at the corresponding position in the matrix, indicating that this transition occurs once. This process is repeated within the range of all particles and the entire time series, and finally a complete state transition probability matrix is formed, providing a quantitative basis for subsequent normalization processing, steady-state probability calculation, and Monte Carlo trajectory reconstruction.

[0147] Step 3.4.3. Normalize the state transition probability matrix to describe state transitions in terms of probabilities. The specific process is as follows: Normalize the number of transitions in each row (divide by the sum of that row) to obtain the normalized state transition probability matrix . The normalized matrix is used to reflect the probability distribution of a particle transitioning from any state to other states, laying the foundation for the subsequent Monte Carlo (MC) method.

[0148] Step 3.4.4. Calculate the steady-state probability distribution of the normalized state transition probability matrix, that is, the state distribution after long-term evolution:

[0149] (11);

[0150] wherein, is the steady-state probability distribution; is the normalized state transition probability matrix; represents the steady-state probability distribution in each state. Equation (11) means that: in the long-term simulation process, starting from any initial state, after multiple-step transitions, it will tend to the steady-state probability distribution . This equation is the mathematical expression of the steady-state condition, reflecting the convergence of the state probability distribution.

[0151] Step 3.5. Reconstruct the microscopic trajectory of CO2 molecules using the Monte Carlo method to simulate its multi-step transition process in the discrete state space; the specific process is as follows:

[0152] Step 3.5.1. First, initialize the spatial state of each particle at the initial time step and set a random number in the interval [0,1] as the transition criterion;

[0153] Step 3.5.2. Subsequently, step by step for each particle in the time series, use the normalized state transition probability matrix in the previous step 3.4.3 to calculate the cumulative distribution function (CDF) of the current state corresponding row. By comparing the random number with the cumulative transition probability, determine whether the particle makes a transition and its target state number. Obtain the transition probability corresponding to the current state from the normalized state transition probability matrix. If the transition probability is zero (i.e., the particle falls into a non-transition state), then retain the previous state to avoid trajectory interruption and ensure the continuity and rationality of the simulation. This Monte Carlo sampling method is based on the Markov process, allowing particles to freely jump in the state space according to statistical laws, thereby quickly generating molecular motion trajectories close to the actual distribution and can be used to further calculate interface properties such as residence time, diffusion path, and contact angle.

[0154] Specifically in implementation, starting from the initial state obtained from MD, using the normalized transition matrix, through the cumulative distribution function (cumsum) and random number judgment, simulate the random jump of atoms between states, thereby generating the MC trajectory (stored in MC_state_eachatom). Count the occupancy times of each discrete state under the MC simulation and store them in the MC_state_acum array. Determine the target state of the next jump by comparing the random number with the cumulative transition probability function (i.e., the cumsum result).

[0155] Using the auxiliary function State_to_Coordinate, the discrete state numbers are mapped to specific coordinates in the physical space. The density distributions in two dimensions (e.g., the y-z plane) and one dimension (e.g., along the z-axis) are respectively counted to obtain the spatial distribution information under the MD and MC simulations. This provides a data basis for the subsequent analysis of the interfacial contact angle.

[0156] Step 3.5.3. Use the Metropolis-Hastings algorithm to perform probability sampling to simulate the transition process of CO2 molecules. This algorithm determines the state transition of molecules through the sampling probability. Calculate using formula (12):

[0157] (12);

[0158] where, is the probability of accepting the new state; , are the transition probabilities of the new and old states respectively;

[0159] Calculate the residence time of CO2 molecules on the solid surface, and statistically analyze the state transition of particles through the time series:

[0160] (13);

[0161] where, is the total number of time steps that the particle continuously appears in a certain specific state; represents the time step when the state number of the particle; is the initial state number of the particle; is the Dirac function; is the total number of time steps.

[0162] Statistically analyze the total number of time steps that the particle continuously resides in a certain specific state as a measure of the stability of the particle in this state.

[0163] Estimate the contact angle according to formula (3);

[0164] Calculate the transition rate:

[0165] (14);

[0166] where, is the transition rate;

[0167] Finally, estimate the dynamic behavior of CO2 molecules and their interaction with the solid surface based on the transition rate.

[0168] The present invention provides a method for calculating the contact angle of CO2 based on two-scale molecular simulation. Combining with the MDMC algorithm, it can improve the calculation efficiency while ensuring the calculation accuracy, and is widely applicable to the fields of EOR, CCUS and other micro-wetting behavior research.

[0169] To prove the feasibility and superiority of the present invention, the following embodiments are given.

[0170] Figure 2 A simplified "toy model" for state transition modeling is shown. Two typical state regions (state A and state B) are set in the figure, and two potential barrier regions (potential barrier 1 and potential barrier 2) are set to simulate the path of molecules transitioning from state A to state B. This figure clearly shows the transfer trajectory of particles between potential barriers and the typical multi-potential well potential distribution structure, which helps to understand the basic mechanism of state transitions in multi-state systems.

[0171] Figures 3 - 6 The state transition probability matrix and its corresponding state transition process are shown. Figure 3 is the original 100×100 state transition probability matrix, which is an original Markov transition matrix containing 100 microstates. This matrix shows an obvious diagonal block structure, indicating that particles mainly jump within local state clusters. Figure 4 Corresponds to the image of the original state trajectory of CO2 molecules evolving with time steps in the microstate space, and the irregular transition process of particles between multiple states can be observed. Figure 5 Shows the 10×10 state transition probability matrix after coarse-graining the original state transition probability matrix, effectively retaining the transition trend and reducing the matrix complexity. Figure 6 Is the time trajectory diagram of the coarse-grained state, showing a clearer transition mode and residence characteristics, which helps to reveal the transfer behavior of particles at the macroscopic scale and rare state transition events.

[0172] Figure 7 A schematic diagram of the contact angle configuration of CO2 molecules in the wetting behavior on the quartz surface is shown. The contact angle is calculated by fitting the interface contour curve in the figure , which can be used for quantitative analysis of the interface wetting performance and provides a basis for subsequent modeling of interfacial tension and adsorption behavior.

[0173] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by those skilled in the art within the essence of the present invention should also fall within the protection scope of the present invention.

Claims

1. A method for calculating the contact angle of carbon dioxide based on two-scale molecular simulation, characterized in that, It includes the following steps: Step 1: Perform all-atom molecular dynamics calculations; Step 2: Construct a coarse-grained molecular dynamics potential energy model, conduct coarse-grained simulations, and predict the wetting behavior of the CO2 interface; Step 3: Combine molecular dynamics and Monte Carlo methods to construct a molecular dynamics-Monte Carlo algorithm for two-scale molecular simulations. The specific process is as follows: Step 3.1: Data collection and preprocessing for molecular dynamics simulations. The data collection process is as follows: By parsing the trajectory file generated in Step 1.1.4, extract the three-dimensional coordinate information of all particles and the simulation boundary information at each time step in the molecular dynamics simulation. The preprocessing process is as follows: First, eliminate the non-numerical tags in the three-dimensional coordinate information, then perform format conversion, and finally store it in a three-dimensional data structure; Step 3.2: Construct a discretized grid and divide the states, conduct statistical analysis of the molecular dynamics simulation trajectories, and construct the original state statistical matrix; Step 3.3: Construct a coordinate expansion matrix; Step 3.4: Construct a state transition probability matrix to record the frequency of particles transitioning from one state to another; Step 3.5: Use the Monte Carlo method to reconstruct the microscopic trajectories of CO2 molecules to simulate their multi-step transition process in the discrete state space.

2. The method for calculating the contact angle of carbon dioxide based on two-scale molecular simulation according to claim 1, wherein The specific process of Step 1 is as follows: Step 1.1: Select the CO2 molecular force field of the phase equilibrium transferable potential model and construct the CO2 / solid interface system; Step 1.2: Perform all-atom scale molecular dynamics simulations to obtain the interface properties. The interface properties include the adsorption energy of CO2 molecules and the contact angle of CO2 molecules on the solid surface.

3. The carbon dioxide contact angle calculation method based on two-scale molecular simulation according to claim 2, characterized in that, The specific process of Step 1.1 is as follows: Step 1.1.1: The selected CO2 molecular force field includes the Lennard-Jones potential energy parameters and the relevant charge distribution information. The Lennard-Jones potential energy parameters include the potential well depth and the inter-particle distance; Step 1.1.2: Use the Lennard-Jones potential energy parameters to construct the CO2 intermolecular interaction potential energy function and calculate the inter-particle interaction potential energy; (1)； Among them, is the potential energy of the interaction between particles; is the distance between two particles; is the depth of the potential well; is the distance between particles; Step 1.1.3: Determine the crystal structure of the solid surface and establish an atomic model; Step 1.1.4: Use a large-scale atomic / molecular massively parallel simulator to perform all-atom scale structural modeling and energy minimization operations. After this operation, the CO2 / solid interface system is constructed. The large-scale atomic / molecular massively parallel simulator synchronously generates a trajectory file.

4. The method for calculating the contact angle of carbon dioxide based on two-scale molecular simulation according to claim 3, wherein In Step 1.2, the calculation formula for the adsorption energy of CO2 molecules is: (2)； Among them, is the adsorption energy of CO2 molecules; is the total energy of the entire simulation system; is the independent energy of CO2 molecules; is the individual energy of the solid surface; The calculation formula for the contact angle of CO2 molecules on the solid surface is: (3)； wherein, is the contact angle of the CO2 molecule on the solid surface; is the solid-gas interfacial tension; is the solid-liquid interfacial tension; is the liquid-gas interfacial tension.

5. The method for calculating the contact angle of carbon dioxide based on two-scale molecular simulation according to claim 4, wherein The specific process of Step 2 is as follows: Step 2.1: Use the Boltzmann inversion method for potential energy fitting, construct a coarse-grained molecular dynamics potential energy model, and optimize it. The formula for the coarse-grained molecular dynamics potential energy model is: (4)； Among them, is the effective interaction potential energy function between coarse particles; is the Boltzmann constant; is the system temperature; is the radial distribution function obtained from all-atom simulation; Optimize the coarse-grained molecular dynamics potential energy model by the force matching method. During the optimization process, construct the following objective function to measure the difference in particle forces between the all-atom model and the coarse model: (5)； Among them, is the objective function; is the true force calculated for the -th particle in the all-atom model simulation; is the predicted force calculated for the -th particle under the current coarse-grained model parameters ; is the total number of simulated particles; Step 2.2: Conduct coarse-grained simulation to predict the wetting behavior of CO2 molecules at the interface; the wetting behavior of CO2 molecules at the interface specifically refers to the contact angle of CO2 molecules on the solid surface; calculate the contact angle formed by CO2 molecules on the solid surface according to formula (3). The Monte Carlo method is used to pre-optimize the initial spatial distribution of CO2 molecules. The specific process is as follows: Under the condition of fixing the effective interaction potential energy function between coarse particles, the positions of CO2 molecules are randomly perturbed through an acceptance-rejection strategy. Combining with the Metropolis criterion, the probability of accepting a new state is: (6)； Among them, is the probability of accepting the new state; is the exponential function with as the base; , are the current position of the particle and the candidate position after perturbation, respectively.

6. The method for calculating the contact angle of carbon dioxide based on two-scale molecular simulation according to claim 5, wherein The specific process of step 3.2 is as follows: Step 3.2.1: To discretize the molecular motion trajectories in three-dimensional space, a unified grid system needs to be constructed first and all particle states numbered. First, traverse the z-direction coordinate data of all particles in the trajectory file, and use the minimum value function and the maximum value function to determine the lowest position and the highest position that the particles have appeared at during the entire simulation process. Then, perform a unified translation on the z-direction coordinate data of all particles, that is, subtract from each z-direction coordinate, so as to translate the entire coordinate interval to a positive number interval starting from zero . Step 3.2.2: Discretely divide the entire three-dimensional space according to the size of the simulation box and the set grid size. After division, several states are obtained, and the state matrix is initialized accordingly; specifically, the grid division method is used to convert the CO2 molecular trajectory into discrete states; meanwhile, calculate the probability of a particle transitioning from one state to another state, which is used to guide the statistics of the state transition frequency. The formula is: (7)； Among them, is the probability that the particle jumps from state to state ; is the number of times the particle jumps from state to state ; is the total number of states; Step 3.2.3: Read the molecular dynamics simulation trajectory data, traverse the spatial coordinates of each particle at each time step, and sequentially extract its position data in the x, y, and z directions; perform periodic boundary condition processing and normalized coordinate transformation on the atomic coordinates. Then, calculate the spatial state number where the particle is located according to the grid size in each direction, store this number in the molecular trajectory array for marking, and at the same time accumulate the number of occurrences in the statistical matrix under the corresponding state index. Finally, the original state statistical matrix is obtained. The original state statistical matrix is specifically defined as , indicating the cumulative number of occurrences of the discrete state numbered in the molecular dynamics simulation; The formula for processing atomic coordinates using periodic boundary conditions is: (8)； Among them, , , respectively represent the new atomic coordinates in the x, y, and z directions after processing; , , are the initial atomic coordinates in the x, y, and z directions respectively; , , are the lengths of the simulation box in the x, y, and z directions respectively; represents the rounding operation on floating-point numbers.

7. The method for calculating the contact angle of carbon dioxide based on two-scale molecular simulation according to claim 6, wherein The specific process of step 3.3 is as follows: based on the original state statistical matrix, three columns are added to store the spatial center coordinates of each state, and the coordinate-expanded matrix obtained after expansion is denoted as ; where , , are the x, y, and z coordinates of the grid cell corresponding to the discrete state numbered respectively; Perform coordinate reduction operations on each state number in sequence. Through the set boundary size of the simulation box and the grid cell size, call the state-coordinate mapping function to map the state number to its geometric center coordinates in the three-dimensional grid, and store the corresponding coordinate values in the x, y, and z directions into the coordinate extension matrix. The calculation formula of the state-coordinate mapping function is: (9)； Among them, is the state number of the particle; , , are the discrete indices of the grid cell where the current particle is located in the x, y, and z directions, respectively; , are the total number of cells in the grid in the x and y directions, respectively.

8. The method for calculating the contact angle of carbon dioxide based on two-scale molecular simulation according to claim 7, wherein The specific process of step 3.4 is as follows: Step 3.4.

1. Traverse the state changes of each atom between consecutive time steps, and count the number of transitions from state to state . Initialize the sparse transition matrix according to the total number of states divided in Step 3.2.

2. (10)； Among them, represents the total number of times transferred from state to state ; is the total number of time steps; represents the state number of the particle at time step ; is the state transition time interval; is the Kronecker function; Step 3.4.2: Traverse the state evolution trajectories of all particles at each time step, extract the adjacent two state numbers at a fixed time interval, and use them as the starting state and target state of the transition respectively. Calculate the probability of each particle transitioning from one state to another according to formula (5); every time it is recognized that a particle transitions from state to state , accumulate the value at the corresponding position in the matrix, indicating that this transition occurs once; this process is repeated for all particles and the entire time series range, and finally a complete state transition probability matrix is formed; Step 3.4.3, normalize the state transition probability matrix to describe state transitions in the form of probabilities; the specific process is as follows: perform normalization processing on the number of transitions in each row to obtain the normalized state transition probability matrix ; Step 3.4.4: Calculate the steady-state probability distribution of the normalized state transition probability matrix. (11)； wherein, is the steady-state probability distribution; represents the steady-state probability distribution in each state.

9. The method for calculating the contact angle of carbon dioxide based on two-scale molecular simulation according to claim 8, wherein The specific process of step 3.5 is as follows: Step 3.5.1: Initialize the spatial state of each particle at the initial time step and set a random number in the interval [0,1] as the transition criterion. Step 3.5.2: Progressively advance each particle in the time series. Use the normalized state transition probability matrix in step 3.4.3 to calculate the cumulative distribution function of the corresponding row of the current state; by comparing the random number with the cumulative distribution function, determine whether the particle transitions and its target state number; obtain the transition probability corresponding to the current state from the normalized state transition probability matrix. If the transition probability is zero, retain the previous state to avoid trajectory interruption. Step 3.5.3: Use the Metropolis-Hastings algorithm for probability sampling to simulate the transition process of CO2 molecules; the Metropolis-Hastings algorithm determines the state transition of molecules through sampling probability. The formula is: (12)； wherein, is the probability of accepting the new state; , are the transition probabilities of the new and old states, respectively; Calculate the residence time of CO2 molecules on the solid surface and statistically analyze the state transitions of particles through the time series. (13)； Among them, is the total number of time steps for the continuous occurrence of particles in a certain specific state; represents the time step and the state number of the particle at that time step; is the initial state number of the particle; is the Dirac function; Calculate the transition rate. (14)； Among them, is the transition rate; Finally, estimate the dynamic behavior of CO2 molecules and their interaction with the solid surface based on the transition rate.

Images