Method for accurate calculation of mean square displacement based on polytype identification of molecular dynamics simulation

Abstract

The application provides a method for accurately calculating mean square displacement (MSD) of molecular dynamics simulation based on multi-phase state identification, which accurately classifies target molecules in the system into four states: near impurity target molecules, target molecules in hydrates, target molecules in clusters and target molecules in liquid phase. The classification process adopts a multi-stage screening strategy to ensure the independence of each type of molecule. Traditional MSD calculation only requires molecules to be in the target phase state at the start and end frames, but the intermediate frame may undergo phase transition, such as entering the hydrate from the liquid phase, resulting in a deviation of the calculation result from the intrinsic diffusion behavior. The application tracks the occurrence state of the molecule during the entire simulation process, accurately finds the state that needs to be calculated, and ensures that MSD only reflects the continuous diffusion process of the molecule in the pure liquid phase through a state persistence verification mechanism, excluding the interference of phase transition, and performing targeted MSD calculation.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of material calculation, and in particular to a molecular dynamics simulation mean square displacement accurate calculation method based on multi-phase state identification. BACKGROUND

[0002] Mean square displacement (MSD) is an important index to characterize the diffusion behavior and migration ability of molecules in molecular dynamics (MD) simulation. By statistically analyzing the trajectories of molecules in the system, MSD can reflect the diffusion coefficient and micro-dynamic characteristics of the substance under different conditions. The existing MSD calculation method is generally based on the time series trajectory of the whole molecule to calculate the mean square displacement, and is widely used in molecular simulation research of liquids, gases, solutions and multi-phase systems.

[0003] However, in a multi-phase system, the traditional MSD calculation method has the following problems: 1. Lack of phase state identification ability: the existing calculation method generally does not distinguish between molecules in the system and directly counts their trajectories. This results in the displacement information of the target molecule being included in the statistics when the target molecule undergoes a phase transition (such as entering the gas phase from the liquid phase, or being adsorbed onto the solid surface from the liquid phase) during simulation, thereby mixing the motion characteristics of different phases. 2. Lack of state persistence verification: even if the existing method roughly distinguishes the phase state by dividing the space, it is difficult to ensure the "state persistence" of the target molecule, that is, it is not possible to extract only the continuous period of the target molecule in the liquid phase for calculation, and the results are often affected by the frequent inter-phase motion of the molecule, causing discontinuity and statistical error in the MSD calculation process.

[0004] Therefore, the existing MSD calculation method has poor applicability in complex multi-phase systems, especially when it is necessary to accurately characterize the diffusion characteristics of liquid molecules, it is difficult to avoid the deviation introduced by the change of the phase state of the molecule. SUMMARY

[0005] Therefore, the present application aims to provide a molecular dynamics simulation mean square displacement accurate calculation method based on multi-phase state identification, which 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, only the displacement of the target molecule in the liquid phase is counted and the MSD is calculated, thereby effectively eliminating the interference caused by the inter-phase motion, and greatly improving the accuracy and applicability of the MSD calculation.

[0006] To achieve the above purpose, the technical scheme of the present application is as follows:

[0007] A molecular dynamics simulation mean square displacement accurate calculation method based on multi-phase state identification, the method comprising the following steps:

[0008] S1. Molecular classification:

[0009] 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;

[0010] S2, Displacement Vector Calculation:

[0011] 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.

[0012] S3 and MSD core computing:

[0013] 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.

[0014] S32, Incremental Calculation:

[0015] 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: ;

[0016] 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 Molecular accumulated displacement square of the same state in the time period; N ) is The number of molecules in the same state in the time period.

[0017] Wherein, the impurities (including non-crystals such as polymers and crystals such as silicon dioxide and metals) in the molecular classification have a certain degree of influence on the movement of the target molecules nearby which need to be calculated, and need to be distinguished from other phase target molecules.

[0018] The cage structure formed by water molecules has specific geometric characteristics, such as pentagonal and hexagonal faces. The application identifies the water cage by analyzing the local structure of water molecules, thereby realizing the identification of the water hydrate crystal structure.

[0019] Further, the identification method of the near-impurity target molecule is: reading the atomic coordinates of the impurities and the coordinates of all target molecules in each frame of trajectory, calculating the periodic boundary correction distance of each target molecule and all atoms of the impurities, 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; the impurities are any one of polymers, silicon dioxide or metal.

[0020] Further, the identification method of the water hydrate crystal structure is: calculating the dihedral angle of each oxygen atom of the water molecule and the oxygen atom of the adjacent water molecule, and the distance between the oxygen atom and its adjacent oxygen atom is <0.35 nm; based on the DFS algorithm, the oxygen atoms satisfying the dihedral angle condition 0.819152<cosθ<1.0 are marked as the same cluster, the size of each cluster is counted, and the largest cluster is selected as the water hydrate crystal structure.

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

[0022] Further, the identification method of the target molecule in the water hydrate is: excluding the near-impurity target molecule which has been marked, and reading the oxygen atom coordinates of the water molecules in the water hydrate crystal structure and the coordinates of all target molecules, calculating the PBC correction distance between the target molecule and each water hydrate oxygen atom, if the distance between the target molecule and any water hydrate oxygen atom is less than 0.5 nm, it is marked as a target molecule in the water hydrate.

[0023] Further, the identification method of the target molecule in the cluster is: taking the remaining target molecules which are not marked as near-impurity target molecules or target molecules in the water hydrate as nodes, and connecting if the distance between the nodes is less than 0.5 nm; traversing all nodes, if it is not marked, starting DFS, recursively visiting all adjacent nodes and marking them as the same connected domain, selecting the connected domain with the most nodes as the cluster, and marking the molecules in the cluster as the target molecules in the cluster.

[0024] Further, S2 also needs to make dynamic box size correction before calculating displacement vector, and the correction method is as follows:

[0025] I. The box sizes of adjacent two frames are averaged in x, y and z directions respectively;

[0026] In x direction, ; wherein, L x,t and L x,t+1 are the sizes of box boundary in x direction at time t and t+1 respectively; L x is the average size of box in x direction taken for calculating t frame and t+1 frame;

[0027] In y and z directions, the corresponding information of x direction is replaced to obtain L y , L z ;

[0028] II. The original displacement components of molecules in adjacent frames in x, y and z directions are calculated;

[0029] In x direction, , wherein, x t , x t+1 are the x coordinates of target molecules in t frame and t+1 frame respectively; x is the x direction displacement of target molecules in t frame and t+1 frame;

[0030] In y and z directions, the corresponding information of x direction is replaced to obtain y , z ;

[0031] III. Correction rule:

[0032] If , the actual displacement is , and the molecule crosses from the left boundary to the right boundary of the box;

[0033] If , the actual displacement is , and the molecule crosses from the right boundary to the left boundary of the box;

[0034] Otherwise, unchanged;

[0035] The correction rule in y and z directions is the same as that in x direction.

[0036] Further, the specific method of state persistence verification is as follows:

[0037] S311, state marking: mark the phase of each target molecule in each frame, liquid phase = 1, others = 0, and form a state matrix: frame x molecule;

[0038] S312, continuous liquid phase section identification: for each molecule, scan its state sequence, and identify all frame intervals that are continuously in the liquid phase;

[0039] S313, MSD calculation condition: only when the molecule is in the liquid phase in all frames in the time interval , the displacement of the molecule is included in the MSD calculation of time.

[0040] Compared with the prior art, the multi-phase state recognition-based accurate calculation method of molecular dynamics simulation mean square displacement of the application has the following advantages:

[0041] (1) The multi-phase state recognition-based accurate calculation method of molecular dynamics simulation mean square displacement of the application, aiming at the phase state mixed pollution problem existing in the traditional mean square displacement (MSD) calculation method in the complex multi-phase system, 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 liquid phase, and a multi-stage screening strategy is used in the classification process to ensure the independence of each type of molecule.

[0042] (2) The traditional MSD calculation only requires the molecule to be in the target phase at the starting frame and the ending frame, but the intermediate frame may experience a phase change, such as entering a hydrate from a liquid, which causes the calculation result to deviate from the intrinsic diffusion behavior. The application ensures that MSD only reflects the continuous diffusion process of molecules in pure liquid phase by a state persistence verification mechanism, excluding phase change interference. BRIEF DESCRIPTION OF DRAWINGS

[0043] The accompanying drawings, which form a part of this application, are used to provide a further understanding of the application, and the illustrative embodiments of the application and their description serve to explain the application. The accompanying drawings do not constitute an improper limitation on the application. In the drawings:

[0044] Figure 1 is a flowchart of the multi-phase state recognition-based accurate calculation method of molecular dynamics simulation mean square displacement of the application;

[0045] Figure 2 is a dynamic box size correction schematic diagram; (a) and (b) are uncorrected, and (c) and (d) are corrected;

[0046] Figure 3 is a persistent state verification mechanism schematic diagram; the state between (a) and (b) is counted in the calculation, and the state of (c) is not counted in the calculation;

[0047] Figure 4 ​The model of hydrate crystal decomposition in the presence of PVP polymer; (a) is the initial model at 0 ns, (b) is the model at 100 ns;

[0048] Figure 5 Comparison of the MSD calculation results of liquid methane in the hydrate decomposition system in the presence of PVP by different MSD calculation methods; the red line is the mean square displacement of the liquid methane molecule calculated by using the traditional MSD calculation method; the blue line is the mean square displacement of the methane molecule in the liquid calculated by the method. DETAILED DESCRIPTION

[0049] It should be noted that the embodiments in the present application and the features in the embodiments can be combined with each other without conflict.

[0050] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the protection scope of the present application.

[0051] The present application proposes a molecular dynamics simulation mean square displacement accurate calculation method based on multi-phase state recognition, including molecular classification, displacement vector calculation and MSD core calculation. The whole processing flow adopts C++ and Python hybrid programming, in which C++ is used for high-performance calculation part such as molecular classification and displacement calculation, and Python is used for process control and data analysis. This hybrid architecture makes full use of the execution efficiency of C++ and the flexibility of Python, and can efficiently process large-scale molecular dynamics simulation trajectories.

[0052] The flow of the method is as shown in Figure 1 The specific steps are as follows:

[0053] S1, molecular classification:

[0054] S11, near impurity target molecule recognition

[0055] Impurities such as polymers and other substances have a certain degree of influence on the movement of the target molecules nearby which need to be calculated, and need to be distinguished from other phase target molecules.

[0056] S111, read the coordinates of impurity atoms and all target molecules in each frame trajectory.

[0057] S112, calculate the periodic boundary correction distance (PBC correction) of each target molecule and all atoms of impurities.

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

[0059] S114, generate a near-impurity target molecule ID list file for subsequent steps.

[0060] S12, hydrate crystal structure identification

[0061] Hydrates are composed of cage-like structures formed by water molecules, which have specific geometric characteristics such as pentagonal and hexagonal faces. The present invention identifies hydrate cages by analyzing the local structure of water molecules.

[0062] S121, dihedral angle feature extraction: calculate the dihedral angle between each water molecule oxygen atom and its adjacent water molecule oxygen atom (distance < 0.35 nm). The characteristic dihedral angle of the hydrate cage corresponds to the pentagonal and hexagonal face structure, and the dihedral angle cosθ is between 0.819152 and 1.0.

[0063] S122, connectivity analysis: based on the depth-first search (DFS) algorithm, mark the oxygen atoms that meet the dihedral angle condition as the same cluster.

[0064] The specific process is: traverse each oxygen atom, if 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.

[0065] S123, maximum cluster screening: count the size of each cluster and select the largest cluster as the hydrate structure (exclude small size false positive signals).

[0066] S124, generate a hydrate oxygen atom ID list file.

[0067] S13, hydrate target molecule identification

[0068] Target molecules located inside the hydrate cage have limited motion and significantly different diffusion behavior from liquid-phase molecules.

[0069] S131, read the coordinates of the hydrate oxygen atoms in the hydrate crystal structure identified in S12 and all target molecule coordinates, excluding the marked near-impurity target molecules.

[0070] S132, calculate the PBC corrected distance between the target molecule and each hydrate oxygen atom.

[0071] S133, if the distance between the target molecule and any hydrate oxygen atom is less than 0.5 nm, mark it as a hydrate target molecule.

[0072] S134, generate a hydrate target molecule ID list.

[0073] S14, target molecule recognition within cluster

[0074] Only target molecules exist within the target molecule cluster, which can form an independent connected domain.

[0075] S141, construct connected graph: take the remaining target molecules (not marked as near impurity target molecules or target molecules in hydrates) as nodes, and connect if the distance between nodes is less than 0.5 nm.

[0076] S142, DFS mark connected domain: traverse all nodes, if not marked, start DFS, recursively visit all adjacent nodes and mark as the same connected domain.

[0077] S143, maximum cluster recognition: select the connected domain with the most nodes as the cluster (exclude small clusters).

[0078] S144, generate target molecule ID list within the cluster.

[0079] S15, target molecule recognition in liquid phase

[0080] Target molecules in liquid phase are the core objects of free diffusion behavior research, and need to exclude the interference of other phase molecules.

[0081] S151, remove near impurity target molecules, target molecules in hydrates and target molecules within clusters from total target molecules, and the remaining are liquid phase target molecules.

[0082] S152, generate liquid phase target molecule ID and coordinate file for subsequent displacement calculation.

[0083] S2, displacement vector calculation:

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

[0085] S21, dynamic box size correction

[0086] In molecular dynamics simulation, the box size may change (such as constant pressure simulation), and traditional PBC correction using fixed box size will introduce errors.

[0087] S211, take the average of the box sizes of the adjacent two frames:

[0088] ,

[0089] In the y and z directions, the corresponding information is replaced by the information in the x direction L y 、 L z .

[0090] S212, calculate the original displacement component of the molecule in the adjacent frame:

[0091] .

[0092] In the y, z direction, the corresponding alternative x direction information gets Δ y , Δ z .

[0093] S213, correction rule:

[0094] If , the actual displacement is , the molecule crosses from the left boundary to the right boundary of the box.

[0095] If , the actual displacement is , the molecule crosses from the right boundary to the left boundary.

[0096] Otherwise, unchanged;

[0097] The correction rule in the y, z direction is the same as that in the x direction.

[0098] The specific dynamic box size correction is shown in Figure 2 .

[0099] Figure 2 (a) and Figure 2 (b) are uncorrected, in which case the absolute value of the displacement Δx generated by the target molecule in the time interval of a frame is greater than the actual displacement distance of the target molecule, which does not conform to the actual situation; Figure 2 (c) and Figure 2 (d) are corrected, in which case the displacement Δx of the target molecule will change with the size of each frame box after correction, which conforms to the actual situation.

[0100] S22, displacement vector generation

[0101] S221, for each pair of adjacent frames (t and t+1):

[0102] Through target molecule ID matching, read the coordinates of the target molecule in the liquid phase in t frame and t+1 frame.

[0103] Apply the above dynamic box size correction to calculate the three-dimensional displacement vector of each target molecule in the liquid phase .

[0104] S222, store the displacement vector: the displacement vector of each molecule is stored in sequence according to the frame, generating a structured data file.

[0105] S3, MSD core calculation:

[0106] S31, State persistence verification:

[0107] Traditional MSD calculation only requires the molecule to be in the target phase at the start and end frames, but intermediate frames can undergo phase transition (e.g., from liquid to hydrate), leading to the deviation of the calculation result from the intrinsic diffusion behavior.

[0108] S311, State labeling: Label the phase of each target molecule in each frame of the liquid phase (1 in liquid phase, 0 otherwise), forming a state matrix (frame x molecule).

[0109] S312, Continuous liquid phase segment identification: For each molecule, scan its state sequence and identify all frame intervals that are continuously in the liquid phase.

[0110] S313, MSD calculation condition: Only when the molecule is in the liquid phase in all frames within the time interval , the displacement of this molecule is included in the MSD calculation of time. Ensure that MSD only reflects the continuous diffusion process of the molecule in the liquid phase, excluding the interference of phase transition.

[0111] The specific state persistence verification is shown as Figure 3 .

[0112] Figure 3 (a) Between t and Figure 3 (b) , the blue target molecule is in the liquid phase, but , the molecule enters the cluster, so only the time period of the molecule from t to is included in the calculation range.

[0113] S32, Incremental calculation:

[0114] Relieve the calculation burden, traditional MSD calculation needs to load the full trajectory displacement data, the memory consumption is large.

[0115] S321, Inter-frame displacement vector storage: Displacement vectors are stored independently by adjacent frames, rather than full trajectory matrices.

[0116] S322, Hierarchical calculation by time interval :

[0117] The outer loop traverses value: .

[0118] The inner loop traverses the start frame t: t = 0, 1, …, N- .

[0119] For each (t, ) combination:

[0120] 1) Scan all target molecules in liquid phase, verify the continuity of state in interval by pre-stored state label;

[0121] 2) Accumulate the square of displacement for the same state molecules in the same time interval:

[0122] .

[0123] 3) Count the number of molecules N with continuous state, calculate MSD:

[0124] .

[0125] The following is an example of calculating the mean square displacement of methane molecules in a methane solution in a four-phase system containing PVP impurities, hydrate crystals, methane gas bubbles, and methane solution, i.e., the target molecule is a methane molecule. The specific model is shown in Figure 4 , which shows the decomposition of hydrate and the change of methane diffusion in liquid phase in the system from 0 to 100 ns. The green spherical shape is the methane molecule; the red dot shape is the water molecule; the gray spherical shape is the carbon atom in PVP, the red spherical shape is the oxygen atom in PVP, the blue spherical shape is the nitrogen atom in PVP, and the white spherical shape is the hydrogen atom in PVP.

[0126] The specific steps are as follows:

[0127] S1, molecular classification:

[0128] S11, near PVP target molecule identification

[0129] S111, read the coordinates of PVP atoms and all methane molecules in each frame trajectory.

[0130] S112, calculate the periodic boundary correction distance (PBC correction) between each methane molecule and all PVP atoms.

[0131] S113, if the distance between the methane molecule and any C atom in PVP is less than 1.0 nm, mark it as a near PVP methane molecule.

[0132] S114, generate a near PVP methane molecule ID list file for subsequent steps.

[0133] S12, hydrate crystal structure identification

[0134] S121, dihedral angle feature extraction: calculate the dihedral angle between each water molecule oxygen atom and its adjacent water molecule oxygen atom (distance < 0.35 nm), the characteristic dihedral angle of the hydrate cage corresponds to the five-sided and six-sided face structure, and the dihedral angle cosθ is between 0.819152 and 1.0;

[0135] S122, Connectivity analysis: based on depth-first search (DFS) algorithm, label oxygen atoms that satisfy dihedral angle condition as the same cluster.

[0136] The specific process is: traverse each oxygen atom, if not labeled, start DFS; DFS recursively visits all adjacent oxygen atoms that satisfy distance and angle conditions, and is labeled as the same cluster.

[0137] S123, Maximum cluster screening: count the size of each cluster, and select the largest cluster as the hydrate structure (exclude small size false positive signal).

[0138] S124, Generate hydrate oxygen atom ID list file.

[0139] S13, Hydrate internal methane molecule recognition

[0140] S131, Read hydrate oxygen atom coordinates and all methane molecule coordinates (exclude near PVP methane molecules that have been labeled).

[0141] S132, Calculate the PBC corrected distance between methane molecules and each hydrate oxygen atom.

[0142] S133, If the distance between methane molecules and any hydrate oxygen atom is less than 0.5 nm, label it as hydrate internal methane molecule.

[0143] S134, Generate hydrate internal target molecule ID list.

[0144] S14, Methane molecule cluster recognition

[0145] Methane molecule cluster simulates methane bubbles, and only methane molecules exist in methane bubbles, which can form independent connected domains.

[0146] S141, Construct connectivity graph: take the remaining methane molecules (not labeled as near PVP methane molecules or hydrate internal methane molecules) as nodes, and nodes are connected if the distance between them is less than 0.5 nm.

[0147] S142, DFS label connected domain: traverse all nodes, if not labeled, start DFS, recursively visit all adjacent nodes and label as the same connected domain.

[0148] S143, Maximum bubble recognition: select the connected domain with the most nodes as the cluster (exclude small bubbles).

[0149] S144, Generate bubble internal target molecule ID list.

[0150] S15, Liquid phase methane molecule screening

[0151] S151, remove the near PVP methane molecules, the hydrate methane molecules and the bubble methane molecules from the total methane molecules, the rest is the liquid phase methane molecules.

[0152] S152, generate the liquid phase methane molecule ID and coordinate file, for subsequent displacement calculation.

[0153] S2, displacement vector calculation:

[0154] S21, dynamic box size correction

[0155] S211, take the average of the box sizes of the adjacent two frames:

[0156] ,

[0157] In the y and z directions, the corresponding information instead of the x direction is obtained L y 、 L z .

[0158] S212, calculate the original displacement component of the molecule in the adjacent frame:

[0159] ,

[0160] In the y and z directions, the corresponding information instead of the x direction is obtained y , Δ z .

[0161] S213, correction rule:

[0162] If , the actual displacement is , the molecule crosses from the left boundary to the right boundary of the box.

[0163] If , the actual displacement is , the molecule crosses from the right boundary to the left boundary.

[0164] Otherwise, unchanged;

[0165] The correction rule in the y and z directions is the same as that in the x direction.

[0166] S22, displacement vector generation

[0167] S221, for each pair of adjacent frames (t and t+1):

[0168] Through methane molecule ID matching, read the coordinates of the methane molecules in the liquid phase in the t frame and the t+1 frame.

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

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

[0171] S3 and MSD core computing:

[0172] S31. Status Continuity Verification:

[0173] 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).

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

[0175] 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.

[0176] S32, Incremental Calculation:

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

[0178] S322, By time interval Hierarchical calculation:

[0179] outer loop traversal value: .

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

[0181] For each (t, )combination:

[0182] 1) Scan all molecules and verify within the interval using pre-stored state tags. State continuity;

[0183] 2) For molecules in the same continuous state at the same time interval, sum the squares of the displacements:

[0184] .

[0185] 3) Count the number of molecules in the continuous state N ( ), calculate MSD:

[0186] .

[0187] Figure 5 The MSD of liquid methane in the hydrate decomposition system in the presence of PVP is calculated by different MSD calculation methods, and it can be seen from the figure that the result of the traditional MSD calculation method is unstable and has a certain distortion, especially in the process of 80-100 ns. The distortion is serious. It is proved that the MSD result calculated by the method is more stable and more accurate.

[0188] The above only describes the preferred embodiments of the present application and is not used to limit the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A method for accurate calculation of mean square displacement in molecular dynamics simulation based on polyphase state recognition, characterized in that: The method comprises the following steps: S1, molecular classification: The near-impurity target molecules, hydrate crystal structures, target molecules in hydrates, and target molecules in clusters in the system are sequentially identified and marked to form respective ID lists; the near-impurity target molecules, target molecules in hydrates, and target molecules in clusters are removed from the total target molecules, and the remaining target molecules are target molecules in liquid phase, and ID lists are generated; S2, displacement vector calculation: Through target molecule ID matching, the coordinates of the target molecule in the liquid phase at t frame and t+1 frame are read, and the three-dimensional displacement vector of each target molecule in the liquid phase is calculated wherein, Δ x corrected , Δ y corrected , Δ z corrected respectively, are the coordinate change amounts of t frame and t+1 frame in x, y and z directions S3, MSD core calculation: S31, state persistent verification: by scanning the state of each target molecule in each frame of liquid phase, when the molecule is in liquid phase in all frames in the time interval , the displacement of the molecule is included in the MSD calculation of time to ensure that MSD only reflects the continuous diffusion process of the molecule in the liquid phase, excluding phase transition interference; S32, incremental calculation: Scanning all target molecules in the liquid phase, verifying the continuity of the state in the interval using pre-stored state markers; for the same molecules in the same time interval, accumulating the square of the displacement , calculating the number of molecules in the continuous state N , calculating the 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 of claim 1, wherein the method is based on the identification of the polytype. The identification method of the near-impurity target molecules is as follows: the atomic coordinates of impurities and the coordinates of all target molecules in each frame of trajectory are read, the periodic boundary correction distance of each target molecule and all atoms of the impurities is calculated, if the distance between the target molecule and any atom of the impurities is less than 1.0 nm, the target molecule is marked as a near-impurity target molecule; the impurities are any one of polymers, silicon dioxide or metals.

3. The multi-phase identification based molecular dynamics simulation mean square displacement accurate calculation method of claim 1, wherein: The identification method of the hydrate crystal structure is as follows: the dihedral angle between the oxygen atom of each water molecule and the oxygen atom of its adjacent water molecule is calculated, and the distance between the oxygen atom and its adjacent oxygen atom is less than 0.35 nm; based on the DFS algorithm, the oxygen atoms satisfying the dihedral angle condition 0.819152<cosθ<1.0 are marked as the same cluster, the size of each cluster is counted, and the largest cluster is selected as the hydrate crystal structure.

4. The multi-phase identification based molecular dynamics simulation mean square displacement accurate calculation method of claim 3, wherein: The specific method of the DFS algorithm is as follows: each oxygen atom is traversed, and DFS is started if it is not marked; all adjacent oxygen atoms satisfying the distance and angle conditions are recursively accessed by DFS, and are marked as the same cluster.

5. The multi-phase identification based molecular dynamics simulation mean square displacement accurate calculation method of claim 1, wherein: The identification method of the target molecule in hydrate is as follows: the near-impurity target molecules that have been marked are excluded, the oxygen atom coordinates of the water molecules in the hydrate crystal structure and the coordinates of all target molecules are read, the PBC correction distance of the target molecules and the hydrate oxygen atoms is calculated, and if the distance between the target molecule and any hydrate oxygen atom is less than 0.5 nm, the target molecule is marked as a target molecule in hydrate.

6. The multi-phase identification based molecular dynamics simulation mean square displacement accurate calculation method of claim 1, wherein: The identification method of the target molecule in cluster is as follows: the remaining target molecules that are not marked as near-impurity target molecules or target molecules in hydrate are taken as nodes, and the nodes are connected if the distance between them is less than 0.5 nm; all nodes are traversed, DFS is started if a node is not marked, all adjacent nodes are recursively accessed by DFS and are marked as the same connected domain, and the connected domain with the largest number of nodes is selected as a cluster, and the molecules in the cluster are marked as target molecules in cluster.

7. The multi-phase identification based molecular dynamics simulation mean square displacement accurate calculation method of claim 1, wherein: S2 also needs to be corrected for dynamic box size before calculating the displacement vector, and the correction method is as follows: Ⅰ, the box sizes of adjacent two frames are averaged in x, y and z directions respectively; In the x direction, ; wherein, L x,t and L x,t+1 respectively the size of the box boundary in the x direction at time t and at time t+1 ; L x is the average of the size of the box taken in the x direction at time t and at time t+1. In the y, z direction, the information corresponding to the x direction is obtained L y , L z ; Ⅱ, the original displacement components of the molecules in adjacent frames in x, y and z directions are calculated; In the x direction, wherein, x t , x t+1 is the target molecule x coordinate in frame t and frame t+1; Δ x is the x direction displacement of the target molecule in frame t and frame t+1; In the y, z direction, the information corresponding to the x direction is replaced by Δ y , Δ z ; Ⅲ, the correction rule: If then the actual displacement is the molecule has crossed from the left to the right border of the box. If then the actual displacement is the molecule crosses from the right boundary to the left boundary; Otherwise, unchanged; The correction rules in y and z directions are the same as in x direction.

8. The multi-phase identification based molecular dynamics simulation mean square displacement accurate calculation method of claim 1, wherein: The specific method of the state continuous verification is as follows: S311, state marking: each target molecule in each frame is marked with a phase state, liquid phase = 1, and others = 0, to form a state matrix: frame x molecule; S312, continuous liquid phase section identification: for each molecule, the state sequence thereof is scanned, and all frame intervals continuously in liquid phase are identified; S313, MSD calculation condition: only when the molecule is in liquid phase in all frames within the time interval , the displacement of the molecule is included in the MSD calculation. time.

Images