Molecular dynamics simulation mean square displacement accurate calculation method based on multiphase identification

Abstract

The invention provides a molecular dynamics simulation mean square displacement accurate calculation method based on multi-phase state recognition, and the method comprises the steps: precisely classifying target molecules in a system into four states: near-impurity target molecules, target molecules in hydrates, target molecules in clusters and target molecules in a liquid phase, and employing a multi-stage screening strategy in the classification process, and the independence of each type of molecules is ensured. Traditional MSD calculation only requires molecules to be in a target phase state in a start frame and an end frame, but a middle frame may undergo phase change, such as entering hydrate from a liquid phase, resulting in deviation of a calculation result from an intrinsic diffusion behavior. According to the method, the occurrence state of the molecules in the whole simulation process is tracked, the state needing to be calculated is accurately found out, it is ensured that the MSD only reflects the continuous diffusion process of the molecules in the pure liquid phase through a state continuous verification mechanism, phase change interference is eliminated, and targeted MSD calculation is conducted.

Description

Technical Field

[0001] This invention relates to the field of materials computation, and in particular to a method for accurate calculation of mean square displacement in molecular dynamics simulations based on multiphase identification. Background Technology

[0002] Mean square displacement (MSD) is an important indicator for characterizing molecular diffusion behavior and migration ability in molecular dynamics (MD) simulations. By statistically analyzing molecular trajectories within a system, MSD can reflect the diffusion coefficient and microscopic dynamic characteristics of a substance under different conditions. Existing MSD calculation methods are typically based on the time-series trajectories of the entire molecule and are widely used in molecular simulation studies of liquids, gases, solutions, and multiphase systems.

[0003] However, in multiphase systems, traditional MSD calculation methods suffer from the following problems: 1. Lack of phase identification capability: Existing calculation methods generally do not distinguish between molecules in the system and directly count their trajectories. This leads to the displacement information of target molecules when they undergo phase transitions during the simulation (such as moving from the liquid phase to the gas phase, or being adsorbed from the liquid phase to the solid surface), thus incorporating motion characteristics from different phases. 2. Lack of state persistence verification: Even if existing methods roughly distinguish phases by dividing spatial regions, it is difficult to guarantee the "state persistence" screening of target molecules. That is, it is impossible to extract only the continuous time period in the liquid phase for calculation. The results are often affected by frequent interphase movements of molecules, causing discontinuities and statistical errors in the MSD calculation process.

[0004] Therefore, existing MSD calculation methods are poorly applicable to complex multiphase systems, especially when it is necessary to accurately characterize the diffusion properties of liquid phase molecules, it is difficult to avoid the bias introduced by changes in molecular phase. Summary of the Invention

[0005] In view of this, the present invention aims to propose an accurate mean square displacement (MSD) calculation method for molecular dynamics simulation based on multiphase state identification. This method can identify the phase state of the target molecule in real time during the MSD calculation process and verify the state persistence of the target molecule in the liquid phase. Displacement statistics and MSD calculation are only performed during the period when the molecule is in the liquid phase, thereby effectively eliminating the interference caused by cross-phase motion and significantly improving the accuracy and applicability of MSD calculation.

[0006] To achieve the above objectives, the technical solution of the present invention is implemented as follows: A method for accurate calculation of mean square displacement in molecular dynamics simulations based on multiphase state identification, comprising the following steps: S1. Molecular classification: Near-impurity target molecules, hydrate crystal structures, target molecules within hydrates, and target molecules within clusters in the system are identified and labeled sequentially to form their respective ID lists; near-impurity target molecules, target molecules within hydrates, and target molecules within clusters are removed from the total target molecules, and the remaining molecules are the target molecules in the liquid phase, and an ID list is generated; S2, Displacement Vector Calculation: By matching target molecule IDs, the coordinates of target molecules in the liquid phase at frames t and t+1 are read, and the three-dimensional displacement vector of each target molecule in the liquid phase is calculated. , where Δ x corrected Δ y corrected Δ z corrected These represent the coordinate changes in the x, y, and z directions for frames t and t+1, respectively. S3 and MSD core computing: S31, State Continuity Verification: By scanning the state of each target molecule in each frame of the liquid phase, when the molecule is within the time interval... When all frames are in the liquid phase, the displacement of the molecule is included in the MSD calculation of time to ensure that the MSD only reflects the continuous diffusion process of the molecule in the liquid phase and eliminates phase transition interference. S32, Incremental Calculation: Scan all target molecules in the liquid phase and verify their position within the specified range using pre-stored state labels. State continuity; for molecules in the same continuous state within the same time interval, the square of the displacement is accumulated. Count the number of molecules in continuous states N ( ), Calculate MSD: ; in, t is the duration of the target molecules being in the same state; t is the initial time when the target molecules are in the same state; k is... At a certain moment within the time period; Δ x k Δ y k Δ z k These represent the changes in the target molecule's coordinates during the time interval from k to k+1; S( )for The cumulative squared displacement of molecules in the same continuous state over a time period; N( )for The number of molecules in the same continuous state within a time period.

[0007] Among them, impurities in the molecular classification (including non-crystals such as polymers and crystals such as silicon dioxide and metals) have a certain degree of influence on the movement of target molecules that need to be calculated nearby, and it is necessary to distinguish them from target molecules in other phases.

[0008] The hydrate is composed of a cage structure formed by water molecules and has specific geometric features, such as pentagonal and hexagonal faces. In the present invention, the hydrate cage is identified by analyzing the local structure of water molecules, and then the identification of the crystal structure of the hydrate is realized.

[0009] Further, the method for identifying target molecules near impurities is as follows: read the atomic coordinates of impurities and the coordinates of all target molecules in each frame of the trajectory, calculate the periodic boundary correction distance between each target molecule and all atoms of the impurities, and if the distance between the target molecule and any impurity atom is less than 1.0 nm, it is marked as a target molecule near impurities; the impurities are any one of polymers, silicon dioxide or metals.

[0010] Further, the method for identifying the crystal structure of the hydrate is as follows: calculate the dihedral angle between each oxygen atom of the water molecule and the oxygen atom of its neighboring water molecule, and the distance between the oxygen atom and its neighboring oxygen atom < 0.35 nm; based on the DFS algorithm, the oxygen atoms that satisfy the dihedral angle condition 0.819152 < cosθ < 1.0 are marked as the same cluster, the sizes of each cluster are counted, and the largest cluster is selected as the crystal structure of the hydrate.

[0011] Further, the specific method of the DFS algorithm is as follows: traverse each oxygen atom, and start DFS if it is not marked; the DFS recursively visits all adjacent oxygen atoms that satisfy the distance and angle conditions and marks them as the same cluster.

[0012] Further, the method for identifying target molecules inside the hydrate is as follows: exclude the marked target molecules near impurities, and read the coordinates of the oxygen atoms of the water molecules and the coordinates of all target molecules in the crystal structure of the hydrate, calculate the PBC correction distance between the target molecule and each oxygen atom of the hydrate, and if the distance between the target molecule and any oxygen atom of the hydrate is less than 0.5 nm, it is marked as a target molecule inside the hydrate.

[0013] Further, the method for identifying target molecules in the cluster is as follows: take the remaining target molecules that are not marked as target molecules near impurities or target molecules inside the hydrate as nodes, and connect them if the distance between the nodes is less than 0.5 nm; traverse all nodes, start DFS if not marked, recursively visit all adjacent nodes and mark them as the same connected domain, and select the connected domain with the largest number of nodes as the cluster, and its molecules are marked as target molecules in the cluster.

[0014] Further, before calculating the displacement vector in S2, dynamic box size correction is also required, and the correction method is as follows: Ⅰ. Take the average value of the box sizes in the x, y, and z directions for two adjacent frames respectively; In the x direction, ;in, L x,t and L x,t+1 These are the dimensions of the box boundary in the x-direction at time t and t+1, respectively; L x To calculate the average size of the box in the x-direction taken at frame t and frame t+1; In the y and z directions, the information corresponding to the x direction is obtained. L y , L z ; II. Calculate the original displacement components of the molecule in the x, y, and z directions in adjacent frames; In the x direction, ,in, x t , x t+1 Let x be the x-coordinate of the target molecule in frame t and frame t+1; Δ x The x-direction displacement of the target molecule in frame t and frame t+1; In the y and z directions, the information corresponding to the x direction is replaced to obtain Δ. y Δ z ; III. Correction Rules: like The actual displacement is The molecule crosses from the left boundary to the right boundary of the box; like The actual displacement is The molecule crosses from the right boundary to the left boundary; otherwise, constant; The correction rules in the y and z directions are the same as those in the x direction.

[0015] Furthermore, the specific method for continuous state verification is as follows: S311, State Marking: Mark the phase state of each target molecule in each frame, with 1 for liquid phase and 0 for others, forming a state matrix: frame × molecule; S312, Continuous liquid phase segment identification: For each molecule, scan its state sequence and identify all frame intervals where all frames are continuously in the liquid phase; S313, MSD calculation conditions: Only when the molecule is within the time interval The displacement of the molecule is only included when all frames are in the liquid phase. Calculation of MSD for time.

[0016] Compared with existing technologies, the method for accurate calculation of mean square displacement in molecular dynamics simulation based on multiphase state identification described in this invention has the following advantages: (1) The molecular dynamics simulation mean square displacement accurate calculation method based on multiphase state identification described in this invention addresses the phase mixing and contamination problem of traditional mean square displacement (MSD) calculation methods in complex multiphase systems. It accurately classifies the target molecules in the system into four states: near-impurity target molecules, target molecules in hydrates, target molecules in clusters, and target molecules in the liquid phase. The classification process adopts a multi-stage filtration strategy to ensure the independence of each type of molecule.

[0017] (2) Traditional MSD calculations only require the molecule to be in the target phase state at the start and end frames, but the intermediate frames may undergo phase transitions, such as from the liquid phase to the hydrate phase, causing the calculation results to deviate from the intrinsic diffusion behavior. This invention ensures that the MSD only reflects the continuous diffusion process of the molecule in the pure liquid phase through a state-based verification mechanism, eliminating phase transition interference. Attached Figure Description

[0018] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings: Figure 1 This is a flowchart of the method for accurate calculation of mean square displacement in molecular dynamics simulation based on multiphase state identification, as described in this invention. Figure 2 This is a schematic diagram of dynamic box size correction; (a) and (b) are the uncorrected cases, and (c) and (d) are the corrected cases; Figure 3 This is a schematic diagram of the continuous state verification mechanism; states (a) to (b) are included in the calculation, while state (c) is not included in the calculation; Figure 4 The model for the decomposition of PVP polymer hydrate crystals is shown in (a) and (b) is the initial model at 0 ns and the model at 100 ns. Figure 5 Comparison of MSD calculation results for liquid methane in a hydrate decomposition system in the presence of PVP using different MSD calculation methods; the red line represents the mean square displacement of liquid methane molecules calculated using the traditional MSD calculation method; the blue line represents the mean square displacement of methane molecules in the liquid phase calculated using this method. Detailed Implementation

[0019] It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other.

[0020] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0021] This invention proposes a method for accurate calculation of mean square displacement (MSD) in molecular dynamics simulations based on multiphase state identification, including molecular classification, displacement vector calculation, and MSD core calculation. The entire process employs a hybrid programming approach using C++ and Python, with C++ used for high-performance computing components such as molecular classification and displacement calculation, and Python used for process control and data analysis. This hybrid architecture fully leverages the execution efficiency of C++ and the flexibility of Python, enabling efficient processing of large-scale molecular dynamics simulation trajectories.

[0022] The process of this method is as follows: Figure 1 As shown, the specific steps are as follows: S1. Molecular classification: S11, Near-impurity target molecule recognition Impurities, such as polymers and other substances, can affect the motion of nearby target molecules that need to be calculated to a certain extent, and it is necessary to distinguish them from target molecules in other phases.

[0023] S111: Read the coordinates of impurity atoms and all target molecules in each frame of the trajectory.

[0024] S112. Calculate the periodic boundary correction distance (PBC correction) between each target molecule and all atoms of the impurities.

[0025] S113. If the distance between the target molecule and any impurity atom is less than 1.0 nm, it is marked as a near-impurity target molecule.

[0026] S114. Generate a list of near-impurity target molecule IDs for subsequent exclusion steps.

[0027] S12, Identification of hydrate crystal structure Hydrates are composed of cage-like structures formed by water molecules and have specific geometric features, such as pentagonal or hexagonal surfaces. This invention identifies hydrate cages by analyzing the local structure of water molecules.

[0028] S121. Dihedral angle feature extraction: Calculate the dihedral angle between each water molecule oxygen atom and its neighboring water molecule oxygen atoms (distance <0.35nm). The characteristic dihedral angle of the hydrate cage corresponds to the pentagonal and hexagonal surface structure, and the dihedral angle cosθ is between 0.819152 and 1.0. S122. Connectivity Analysis: Based on the depth-first search (DFS) algorithm, oxygen atoms that satisfy the dihedral angle condition are marked as belonging to the same cluster.

[0029] The specific process is as follows: traverse each oxygen atom, and if it is not marked, start DFS; DFS recursively visits all adjacent oxygen atoms that meet the distance and angle conditions and marks them as the same cluster.

[0030] S123. Largest Cluster Screening: Count the size of each cluster and select the largest cluster as the hydrate structure (excluding false positive signals of small size).

[0031] S124. Generate a list file of oxygen atom IDs for hydrates.

[0032] S13, Target molecule recognition within hydrates The target molecules located within the hydrate cage have restricted movement, and their diffusion behavior differs significantly from that of liquid-phase molecules.

[0033] S131. Read the coordinates of the oxygen atom in the hydrate crystal structure identified in S12 and the coordinates of all target molecules, and exclude the labeled near-impurity target molecules.

[0034] S132. Calculate the PBC correction distance between the target molecule and each hydrate oxygen atom.

[0035] S133. If the distance between the target molecule and any oxygen atom of the hydrate is less than 0.5 nm, it is marked as a target molecule within the hydrate.

[0036] S134. Generate a list of target molecule IDs within the hydrate.

[0037] S14, Target molecule recognition within clusters The target molecule cluster contains only the target molecule and can form an independent connected domain.

[0038] S141. Construct a connected graph: Use the remaining target molecules (unlabeled as near-impurity target molecules or target molecules within hydrates) as nodes, and connect nodes if the distance between nodes is less than 0.5 nm.

[0039] S142, DFS Marking Connected Components: Traverse all nodes. If a node is not marked, start DFS to recursively visit all adjacent nodes and mark them as the same connected component.

[0040] S143. Identification of the largest cluster: Select the connected component with the most nodes as the cluster (excluding small clusters).

[0041] S144. Generate a list of target molecule IDs within the cluster.

[0042] S15. Target molecule recognition in liquid phase The target molecule in the liquid phase is the core object of study on free diffusion behavior, and interference from molecules in other phases must be eliminated.

[0043] S151. Remove near-impurity target molecules, target molecules within hydrates, and target molecules within clusters from the total target molecules. The remaining molecules are the target molecules in the liquid phase.

[0044] S152. Generate the liquid phase target molecule ID and coordinate file for subsequent displacement calculations.

[0045] S2, Displacement Vector Calculation: This step calculates the displacement vector of the target molecules in the liquid phase between adjacent frames and solves the coordinate jump problem caused by periodic boundary conditions (PBC).

[0046] S21, Dynamic Box Size Correction In molecular dynamics simulations, the box size may vary (e.g., constant pressure simulations), and traditional PBC correction using a fixed box size will introduce errors.

[0047] S211. Take the average value of the box size for two adjacent frames: , In the y and z directions, the information corresponding to the x direction is obtained. L y , L z .

[0048] S212. Calculate the original displacement components of the molecule in adjacent frames: .

[0049] In the y and z directions, the information corresponding to the x direction is replaced to obtain Δ. y Δ z .

[0050] S213, Correction Rules: like The actual displacement is The molecule crosses from the left boundary of the box to the right boundary.

[0051] like The actual displacement is The molecule crosses from the right boundary to the left boundary.

[0052] otherwise, constant; The correction rules in the y and z directions are the same as those in the x direction.

[0053] Specific dynamic box size correction, such as Figure 2 As shown.

[0054] Figure 2 (a) and Figure 2 (b) is the uncorrected case. In this case, the absolute value of the displacement Δx generated by the target molecule in a frame time interval is greater than the actual displacement distance of the target molecule, which does not match the actual displacement. Figure 2 (c) and Figure 2 (d) shows the corrected case. The displacement Δx of the target molecule changes with the size of the box in each frame, which is consistent with the actual situation.

[0055] S22, Displacement Vector Generation S221. For each pair of adjacent frames (t and t+1): By matching the target molecule ID, the coordinates of the target molecule in the liquid phase at frame t and frame t+1 are read.

[0056] Applying the dynamic box-size correction described above, the three-dimensional displacement vector of each target molecule in the liquid phase is calculated. .

[0057] S222, Storing displacement vectors: The displacement vector of each molecule is stored in a frame sequence to generate a structured data file.

[0058] S3 and MSD core computing: S31. Status Continuity Verification: Traditional MSD calculations only require molecules to be in the target phase state at the start and end frames, but intermediate frames may undergo phase transitions (such as transitioning from the liquid phase to the hydrate phase), causing the calculation results to deviate from the intrinsic diffusion behavior.

[0059] S311, State Marking: Mark the phase state of each target molecule in each frame of the liquid phase (=1 in the liquid phase, 0 elsewhere), forming a state matrix (frame × molecule).

[0060] S312, Continuous liquid phase segment identification: For each molecule, scan its state sequence to identify all continuous frame intervals in the liquid phase.

[0061] S313, MSD calculation conditions: Only when the molecule is within the time interval The displacement of the molecule is only included when all frames are in the liquid phase. MSD calculation for time. Ensure that the MSD only reflects the continuous diffusion process of molecules in the liquid phase, eliminating phase transition interference.

[0062] The specific status continues to be verified, such as... Figure 3 As shown.

[0063] Figure 3 (a) From time t to Figure 3 (b) For a moment, the blue target molecule is in the liquid phase, but The molecule enters the cluster at time t; therefore, only the molecule at time t is considered. The time period of a moment is included in the calculation.

[0064] S32, Incremental Calculation: To reduce the computational burden, traditional MSD calculations require loading full trajectory displacement data, resulting in high memory consumption.

[0065] S321. Inter-frame displacement vector storage: Displacement vectors are stored independently for adjacent frames, rather than the entire trajectory matrix.

[0066] S322, By time interval Hierarchical calculation: outer loop traversal value: .

[0067] The inner loop iterates through the starting frame t: t=0,1,…,N- .

[0068] For each (t, )combination: 1) Scan all target molecules in the liquid phase and verify their position within the specified range using pre-stored state labels. State continuity; 2) For molecules in the same continuous state at the same time interval, sum the squares of the displacements: .

[0069] 3) Count the number of molecules in the continuous state N ( ), Calculate MSD: .

[0070] The following example demonstrates the calculation of the mean square displacement of methane molecules in a methane solution within a four-phase system containing PVP impurities, hydrate crystals, methane bubbles, and a methane solution. The target molecule is the methane molecule. A detailed explanation is provided below, using the specific model as follows: Figure 4 As shown, from left to right, the decomposition of hydrates and the diffusion of methane in the liquid phase are illustrated in the 0-100 ns system. The green spheres represent methane molecules; the red dots represent water molecules; the gray spheres represent carbon atoms in PVP; the red spheres represent oxygen atoms in PVP; the blue spheres represent nitrogen atoms in PVP; and the white spheres represent hydrogen atoms in PVP.

[0071] The specific steps are as follows: S1. Molecular classification: S11, near-PVP target molecule recognition S111: Read the PVP atom coordinates and all methane molecule coordinates in each frame trajectory.

[0072] S112. Calculate the periodic boundary correction distance (PBC correction) between each methane molecule and all atoms of PVP.

[0073] S113. If the distance between a methane molecule and any C atom in a PVP is less than 1.0 nm, it is labeled as a near-PVP methane molecule.

[0074] S114. Generate a list of near-PVP methane molecule IDs for subsequent exclusion steps.

[0075] S12, Identification of hydrate crystal structure S121. Dihedral angle feature extraction: Calculate the dihedral angle between each water molecule oxygen atom and its neighboring water molecule oxygen atoms (distance <0.35nm). The characteristic dihedral angle of the hydrate cage corresponds to the pentagonal and hexagonal surface structure, and the dihedral angle cosθ is between 0.819152 and 1.0. S122. Connectivity Analysis: Based on the depth-first search (DFS) algorithm, oxygen atoms that satisfy the dihedral angle condition are marked as belonging to the same cluster.

[0076] The specific process is as follows: traverse each oxygen atom, and if it is not marked, start DFS; DFS recursively visits all adjacent oxygen atoms that meet the distance and angle conditions and marks them as the same cluster.

[0077] S123. Largest Cluster Screening: Count the size of each cluster and select the largest cluster as the hydrate structure (excluding false positive signals of small size).

[0078] S124. Generate a list file of oxygen atom IDs for hydrates.

[0079] S13, methane molecule recognition within hydrates S131. Read the coordinates of oxygen atoms in hydrates and the coordinates of all methane molecules (excluding labeled near-PVP methane molecules).

[0080] S132. Calculate the PBC corrected distance between methane molecules and oxygen atoms of each hydrate.

[0081] S133. If the distance between a methane molecule and any oxygen atom of a hydrate is less than 0.5 nm, it is marked as a methane molecule within the hydrate.

[0082] S134. Generate a list of target molecule IDs within the hydrate.

[0083] S14, Methane Molecular Cluster Recognition Methane molecule clusters simulate methane bubbles, which contain only methane molecules and can form independent connected domains.

[0084] S141. Construct a connected graph: Use the remaining methane molecules (not labeled as near-PVP methane molecules or methane molecules within hydrates) as nodes, and connect nodes if the distance between nodes is less than 0.5 nm.

[0085] S142, DFS Marking Connected Components: Traverse all nodes. If a node is not marked, start DFS to recursively visit all adjacent nodes and mark them as the same connected component.

[0086] S143. Maximum bubble identification: Select the connected component with the most nodes as the cluster (excluding small bubbles).

[0087] S144. Generate a list of target molecule IDs within the bubble.

[0088] S15, Screening of methane molecules in liquid phase S151. Remove near-PVP methane molecules, hydrate methane molecules, and bubble methane molecules from the total methane molecules. The remaining molecules are the methane molecules in the liquid phase.

[0089] S152. Generate liquid methane molecule IDs and coordinate files for subsequent displacement calculations.

[0090] S2, Displacement Vector Calculation: S21, Dynamic Box Size Correction S211. Take the average value of the box size for two adjacent frames: , In the y and z directions, the information corresponding to the x direction is obtained. L y , L z .

[0091] S212. Calculate the original displacement components of the molecule in adjacent frames: , In the y and z directions, the information corresponding to the x direction is replaced to obtain Δ. y Δ z .

[0092] S213, Correction Rules: like The actual displacement is The molecule crosses from the left boundary of the box to the right boundary.

[0093] like The actual displacement is The molecule crosses from the right boundary to the left boundary.

[0094] otherwise, constant; The correction rules in the y and z directions are the same as those in the x direction.

[0095] S22, Displacement Vector Generation S221. For each pair of adjacent frames (t and t+1): By matching the methane molecule IDs, the coordinates of methane molecules in the liquid phase at frames t and t+1 are read.

[0096] Applying the above dynamic box-size correction, calculate the three-dimensional displacement vector of each methane molecule in the liquid phase. .

[0097] S222, Storing displacement vectors: The displacement vector of each molecule is stored in a frame sequence to generate a structured data file.

[0098] S3 and MSD core computing: S31. Status Continuity Verification: S311, State Marking: Mark the phase state of each methane molecule in each frame (1 in liquid phase, 0 elsewhere), forming a state matrix (frame × molecule).

[0099] S312. Continuous liquid phase segment identification: For each methane molecule, scan its state sequence to identify all continuous frame intervals in the liquid phase.

[0100] S313, MSD Calculation Conditions: Only when methane molecules are within the time interval The displacement of the methane molecule is only included when all frames are in the liquid phase. MSD calculation for time. Ensure that the MSD only reflects the continuous diffusion process of molecules in the liquid phase, eliminating phase transition interference.

[0101] S32, Incremental Calculation: S321. Inter-frame displacement vector storage: Displacement vectors are stored independently for adjacent frames, rather than the entire trajectory matrix.

[0102] S322, By time interval Hierarchical calculation: outer loop traversal value: .

[0103] The inner loop iterates through the starting frame t: t=0,1,…,N- .

[0104] For each (t, )combination: 1) Scan all molecules and verify within the interval using pre-stored state tags. State continuity; 2) For molecules in the same continuous state at the same time interval, sum the squares of the displacements: .

[0105] 3) Count the number of molecules in the continuous state N ( ), Calculate MSD: .

[0106] Figure 5 This comparison shows the MSD calculation results of liquid-phase methane in a hydrate decomposition system under PVP using different MSD calculation methods. The figure reveals that the results from traditional MSD calculation methods are unstable and exhibit some distortion, especially in the 80-100 ns range. This demonstrates that the MSD results calculated by the proposed method are more stable and accurate.

[0107] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for accurate calculation of mean square displacement in molecular dynamics simulations based on multiphase state identification, characterized in that: This method includes the following steps: S1. Molecular classification: Identify and label the near-impurity target molecules, hydrate crystal structures, target molecules within hydrates, and target molecules within clusters in the system in sequence, forming their respective ID lists; exclude the near-impurity target molecules, target molecules within hydrates, and target molecules within clusters from the total target molecules, and the remaining are the target molecules in the liquid phase, and generate an ID list; S2. Displacement vector calculation: By matching target molecule IDs, the coordinates of target molecules in the liquid phase at frames t and t+1 are read, and the three-dimensional displacement vector of each target molecule in the liquid phase is calculated. , where Δ x corrected Δ y corrected Δ z corrected These represent the coordinate changes in the x, y, and z directions for frames t and t+1, respectively. S3. Core MSD calculation: S31, State Continuity Verification: By scanning the state of each target molecule in each frame of the liquid phase, when the molecule is within the time interval... When all frames are in the liquid phase, the displacement of the molecule is included in the MSD calculation of time to ensure that the MSD only reflects the continuous diffusion process of the molecule in the liquid phase and eliminates phase transition interference. S32. Incremental calculation: Scan all target molecules in the liquid phase and verify their position within the specified range using pre-stored state labels. State continuity; for molecules in the same continuous state within the same time interval, the square of the displacement is accumulated. Calculate the number of molecules N in the continuous state. ), Calculate MSD: ; in, τ represents the duration of the target molecules being in the same state; t represents the initial moment when the target molecules are in the same state; k represents a moment within the time interval τ; Δ x k Δ y k Δ z k These represent the changes in the target molecule's coordinates during the time interval from k to k+1; S( )for The cumulative squared displacement of molecules in the same continuous state over a time period; N( )for The number of molecules in the same continuous state within a time period.

2. The method for accurate calculation of mean square displacement in molecular dynamics simulation based on multiphase state identification according to claim 1, characterized in that: The method for identifying near-impurity target molecules is as follows: Read the atomic coordinates of impurities and the coordinates of all target molecules in each frame of the trajectory, calculate the periodic boundary correction distance between each target molecule and all atoms of the impurity. If the distance between a target molecule and any impurity atom is less than 1.0 nm, it is labeled as a near-impurity target molecule; the impurity is any one of polymers, silica, or metals.

3. The method for accurate calculation of mean square displacement in molecular dynamics simulation based on multiphase state identification according to claim 1, characterized in that: The method for identifying the hydrate crystal structure is as follows: Calculate the dihedral angle between the oxygen atom of each water molecule and the oxygen atom of its neighboring water molecule, and the distance between this oxygen atom and its near oxygen atom < 0.35 nm; Based on the DFS algorithm, label the oxygen atoms that satisfy the dihedral angle condition 0.819152 < cosθ < 1.0 as the same cluster, count the size of each cluster, and select the largest cluster as the hydrate crystal structure.

4. The method for accurate calculation of mean square displacement in molecular dynamics simulation based on multiphase state identification according to claim 3, characterized in that: The specific method of the DFS algorithm is as follows: Traverse each oxygen atom, and start DFS if it is not labeled; DFS recursively visits all adjacent oxygen atoms that meet the distance and angle conditions and labels them as the same cluster.

5. The method for accurate calculation of mean square displacement in molecular dynamics simulation based on multiphase state identification according to claim 1, characterized in that: The method for identifying target molecules within hydrates is as follows: Exclude the already labeled near-impurity target molecules, and read the coordinates of the oxygen atoms of water molecules and the coordinates of all target molecules in the hydrate crystal structure, calculate the PBC correction distance between the target molecules and each oxygen atom of the hydrate. If the distance between a target molecule and any oxygen atom of the hydrate is less than 0.5 nm, it is labeled as a target molecule within the hydrate.

6. The method for accurate calculation of mean square displacement in molecular dynamics simulation based on multiphase state identification according to claim 1, characterized in that: The method for identifying target molecules within clusters is as follows: Take the remaining target molecules that are not labeled as near-impurity target molecules or target molecules within hydrates as nodes, and connect the nodes if the distance between them is less than 0.5 nm; Traverse all nodes, start DFS if not labeled, recursively visit all adjacent nodes and label them as the same connected domain, and select the connected domain with the largest number of nodes as the cluster, and its molecules are labeled as target molecules within the cluster.

7. The method for accurate calculation of mean square displacement in molecular dynamics simulation based on multiphase state identification according to claim 1, characterized in that: Before calculating the displacement vector in S2, dynamic box size correction is also required, and the correction method is as follows: Ⅰ. Take the average value of the box sizes in the x, y, and z directions for two adjacent frames respectively; In the x direction, ; in, L x,t and L x,t+1 These are the dimensions of the box boundary in the x-direction at time t and t+1, respectively; L x To calculate the average size of the box in the x-direction taken at frame t and frame t+1; In the y and z directions, the information corresponding to the x direction is obtained. L y , L z ; Ⅱ. Calculate the original displacement components of the molecules in the x, y, and z directions for two adjacent frames; In the x direction, ,in, x t , x t+1 Let x be the x-coordinate of the target molecule in frame t and frame t+1; Δ x The x-direction displacement of the target molecule in frame t and frame t+1; In the y and z directions, the information corresponding to the x direction is replaced to obtain Δ. y Δ z ; Ⅲ. Correction rules: like The actual displacement is The molecule crosses from the left boundary to the right boundary of the box; like The actual displacement is The molecule crosses from the right boundary to the left boundary; otherwise, constant; The correction rules in the y and z directions are the same as those in the x direction.

8. The method for accurate calculation of mean square displacement in molecular dynamics simulation based on multiphase state identification according to claim 1, characterized in that: The specific method for continuous state verification is as follows: S311. State labeling: Label the phase state of each target molecule in each frame, liquid phase = 1, others = 0, and form a state matrix: frame × molecule; S312. Identification of continuous liquid phase segments: For each molecule, scan its state sequence to identify the frame intervals where it is continuously in the liquid phase in all frames; S313, MSD calculation conditions: Only when the molecule is within the time interval The displacement of the molecule is only included when all frames are in the liquid phase. Calculation of MSD for time.

Images