Method for constructing nuclear reaction coordinates based on molecular simulation and genetic algorithm

Abstract

The application provides a method for constructing nucleation reaction coordinates based on molecular simulation and a genetic algorithm, taking nucleation reaction transition state resolution as an optimization target, and fusing molecular topological information and the like, combining a data-driven method, so that high-fidelity characterization of nucleation details is obtained, mechanism distortion is avoided, and more comprehensive prediction and regulation capabilities are provided for synthesis of crystal materials; and nucleation reaction coordinates are automatically constructed, without relying on manual trial and error.

Description

Technical Field

[0001] This invention belongs to the field of materials computation, and in particular relates to a method for constructing nucleation reaction coordinates based on molecular simulation and genetic algorithm. Background Technology

[0002] The initial crystallization / nucleation process involves rapid fluctuations in the generation and growth of "crystal nuclei" on a nanometer / nanosecond scale, making them difficult to capture and quantify directly experimentally. Therefore, computational simulations of nucleation events are widely used to help quantify and capture microscopic reaction processes. However, the development of traditional reaction coordinates in computational simulations requires a labor-intensive and trial-and-error process, and relying solely on the ability to distinguish between standard crystals and liquids makes it difficult to guarantee that the coordinates can reasonably reflect the nucleation process, leading to the omission of crucial metastable states during the exploration of nucleation. Summary of the Invention

[0003] In view of this, the present invention aims to propose a method for constructing nuclear reaction coordinates based on molecular simulation and genetic algorithm, which can automatically construct nuclear reaction coordinates without relying on repeated manual trial and error.

[0004] To achieve the above objectives, the technical solution of the present invention is implemented as follows:

[0005] A method for constructing nuclear reaction coordinates based on molecular simulation and genetic algorithm, the method comprising the following steps:

[0006] S1. First, sample the nucleation transition state;

[0007] S11. Simulate the nucleation process using molecular simulation software and obtain the simulation trajectory; then, first use empirical order parameters to obtain preliminary unlabeled configuration data, and then use empirical order parameters to preliminarily determine the crystalline phase C and the liquid phase L. After repeated iterations, obtain the configuration data with labels C / L.

[0008] S12. Subsequently, local structure coding is constructed, and the empirical order parameters in S11 are replaced with local structure coding to obtain crystal nucleus information.

[0009] Local structure encoding transforms the molecular structure in Cartesian coordinates into a binary code that can be processed and iterated by a computer through topological ring structure recognition and topological bond structure definition. This allows for the identification of the crystalline phase C and the liquid phase L, and the acquisition of preliminary crystal nucleus information. Different codes correspond to different crystal nucleus information.

[0010] S2, then perform the genetic algorithm iteration;

[0011] Using a genetic algorithm and the crystal nucleus information obtained from S1, the nucleation reaction coordinates of the binary base generated by the local structure encoding are self-evolved and iterated to obtain an accurate encoding form;

[0012] Then, replace the empirical sequence parameter with the accurate encoding form and execute S11 again to obtain the accurate label C / L configuration data;

[0013] S3. Finally, output the crystal nucleus information;

[0014] The optimized encoding is encapsulated and visualized to process molecular simulation trajectory configurations in batches.

[0015] Furthermore, the specific method of S11 is as follows:

[0016] S111. First, initial sampling is performed, that is, the nucleation process is simulated using molecular simulation software and the simulation trajectory is obtained; empirical order parameters are used to determine whether the two phases CL are connected, to initially describe the nucleation process, and to obtain preliminary unlabeled configuration data.

[0017] S112. Secondly, targetless path sampling is adopted;

[0018] ① Randomly select a frame configuration in the trajectory, redistribute the velocities p of all particles to follow the Boltzmann distribution, and perform molecular dynamics simulations; use empirical order parameters to determine whether to enter the crystalline phase C or the liquid phase L;

[0019] ② Using the initial configuration in ①, assign velocity -p to all particles, perform molecular dynamics simulations, and use empirical order parameters to determine whether to enter the crystalline phase C or the liquid phase L;

[0020] ③ If the final configurations obtained through ① and ② fall into C and L respectively, then the newly generated trajectories are spliced ​​together to obtain a new trajectory, replacing the original trajectory; steps ① and ② are executed again, and the initial configuration is saved and marked as either configuration C or L, each with a 50% probability, determined based on the generated random number;

[0021] ④ If the crystal phase C is entered through both ① and ②, then mark the initial configuration as C, save the initial configuration, and execute steps ① and ② again;

[0022] ⑤ If the liquid phase L has been entered through both ① and ②, then mark the initial configuration as L, save the initial configuration, and execute steps ① and ② again;

[0023] ⑥ Iterate repeatedly to obtain a configuration file with tagged C / L.

[0024] Furthermore, the specific method for local structure encoding is as follows:

[0025] I. First, identify the minimum cycle topology;

[0026] ① Calculate and collect the index number (ID) of each particle's neighboring particles within the cutoff range, and define these neighboring particles as the neighboring particles of the central particle;

[0027] ②Summarize the set of IDs representing the connections between the central particle and its neighboring particles, i.e., {ID} i N, ID ij1 ID ij2 …ID ijn}, where N is ID i The total number of neighboring particles, ID ijn Indicates ID i The neighboring particles, n=1, 2, ...∈N;

[0028] Connection relationship: For example, if atom i is within the cutoff range of r and there is atom j, then atom i and atom j are considered to be connected by some kind of "bond", similar to two water molecules being connected by a hydrogen bond; therefore, the connection relationship here refers to the presence of n other atoms within this range of connection bonds for atom i.

[0029] ③ Utilize the set of connection relationship IDs of neighboring particles to obtain IDs using a recursive algorithm. i The ring topology information for the starting point;

[0030] ④ Eliminate invalid topological cycles to obtain the minimum cycle topological information;

[0031] ⑤ Obtain the minimum topological ring structure of each particle and perform statistical analysis to obtain a new set of topological information, i.e., {ID, N}. 3-ring N 4-ring , ...N r-ring}, r≥3, where N r-ring This represents the number of r-rings that the ID particle is in.

[0032] II. Then define the topological bond structure;

[0033] ①Based on the set of IDs and topological information of neighboring particles in I, define the connection bond between two adjacent particles as B. i-j ;

[0034] ② Define key topology information, B i-j ={a,b,c,…}, where:

[0035] a=N 3-ring (i)·N 4-ring (j)-N 4-ring (i)·N 3-ring (j);

[0036] b=N 4-ring (i)·N 5-ring (j)-N 5-ring (i)·N 4-ring (j);

[0037] c=N 5-ring(i)·N 6-ring (j)-N 6-ring (i)·N 5-ring (j);

[0038] By analogy, information from multiple dimensions can be obtained;

[0039] ③ The types of bond topological information during the statistical molecular simulation are identified and arranged from highest to lowest frequency. These are then randomly mapped to either "0" or "1", i.e., the mapping relationship f(B) i-j → 0 / 1;

[0040] III. Finally, the crystal nuclei and liquid-like particles are classified.

[0041] ① By utilizing the set of connection relationship IDs and topological information of neighboring particles, the bond topological information contained in each particle is obtained, and combined with the mapping relationship, the superimposed bond feature, i.e., B, is obtained. all =∑f(B i-j );

[0042] ②When particle B of i all When the value is ≥3, particle i is labeled as a crystalline particle; otherwise, it is labeled as a liquid-like particle.

[0043] ③ Connect all crystal types according to the ID set of the connection relationship between neighboring particles, obtain the crystal nucleus clusters connected together, and obtain the largest cluster, that is, obtain the preliminary crystal nucleus information.

[0044] Furthermore, in I, the calculation of the ID of adjacent particles also needs to consider periodic boundary conditions; the cutoff range is ≤0.35nm.

[0045] Furthermore, the specific steps of S2 are as follows:

[0046] S21. First, calculate the fitness.

[0047] ① The configuration file of the tag C / L obtained by sampling the nucleation transition state, as well as the local structure code, are used to calculate the crystal nucleus information, which includes the crystal nucleus size CV1 and morphology CV2;

[0048] ② Construct reaction coordinates z(CV1, CV2, ...), which are linear combinations of CVi, i.e.: z = a0 + a1CV1 + a2CV2 + ... + a i CVi, where a i These are coefficients to be determined;

[0049] ③ For the configuration labeled C, construct the probability function:

[0050] P B =(1+erf[z(CV1, CV2,…)]) / 2;

[0051] ④ For the configuration labeled L, construct the probability function:

[0052] P A =1-(1+erf[z(CV1, CV2,…)]) / 2;

[0053] ⑤ Construct the maximum likelihood estimate by multiplying all probabilities together, that is:

[0054] L(z) = ∏P B (a i )·∏P A (a i );

[0055] ⑥ Take the logarithm of L and find the maximum value of log(L(z)), that is, find the maximum value of a. i The maximum value of the function log(L(z)) of the independent variable;

[0056] ⑦ Obtain a i The specific value of the constructed reaction coordinate z(CV1, CV2, ...) is the maximum value of log(L(z)) corresponding to the constructed reaction coordinate z(CV1, CV2, ...). The larger the fitness, the better the performance of the constructed reaction coordinate.

[0057] S22, and then perform genetic iteration;

[0058] ① Perform batch operations on the binary encoding of random mapping in the local structure encoding, arrange the key topology information of different types from high to low according to the frequency of occurrence, and perform random mapping on each structure to obtain binary encoding, thus obtaining different binary encodings;

[0059] ②Calculate the fitness of each encoding separately;

[0060] ③ Randomly select 4 codes and compare their fitness, then pass on the two codes with the highest fitness to the next generation;

[0061] ④ For the selected codes, swap local segments in pairs;

[0062] ⑤ Perform random mutations on the selected codes, that is, mutate "0" to "1" or "1" to "0";

[0063] ⑥ Repeat steps ②-⑤ until the binary code with the highest fitness is obtained;

[0064] S23. Finally, perform sampling iterations;

[0065] ① Determine whether the new binary code obtained by genetic iteration converges with the previous generation's binary code, i.e., |log(L(z)| i+1 )-log(L(z) i)|<0.01;

[0066] ② If convergence is not achieved, use the newly iterated optimal reaction coordinates to perform S1 nucleation transition state sampling to obtain a more accurate labeled configuration dataset, and further perform S21 and S22 to iterate and obtain the final high-precision nucleation reaction coordinates;

[0067] ③ If convergence is achieved, execute the S3 crystal nucleus information output.

[0068] Furthermore, the specific steps of S3 are as follows:

[0069] S31. Encapsulate the iteratively optimized binary code;

[0070] S32. Input the binary code determined by the encapsulation into the local structure code, and output the crystal nucleus information, that is, the number of the crystal nucleus particles that make up the largest connectivity, and visualize the crystal nucleus with VMD or OVITIO.

[0071] S33, Calculate auxiliary thermodynamic and kinetic characteristics.

[0072] Furthermore, auxiliary thermodynamic and kinetic characteristics include crystal nucleus sphericity, the generalized diffusion coefficient of the critical nucleus, the adhesion rate of the nucleus, and the proportion of crystal structure.

[0073] Furthermore, S11 employs high-driving conditions or meta-dynamics to simulate the nucleation process.

[0074] Compared with existing technologies, the method for constructing nucleation reaction coordinates based on molecular simulation and genetic algorithms described in this invention has the following advantages:

[0075] This invention addresses the limitations of traditional nucleation reaction coordinates, which rely on repeated manual trial and error and have limited accuracy. It introduces a genetic algorithm to construct a self-evolving framework for developing nucleation reaction coordinates. Traditional methods often aim to distinguish between the initial liquid phase and the standard crystalline phase, emphasizing the separability of the final state—a "result-oriented" approach. This differs from the "process-oriented" nature of reaction coordinates, leading to insufficient accuracy. This invention, however, optimizes the resolution of the nucleation reaction transition state, integrates molecular topological information, and combines data-driven methods to achieve high-fidelity characterization of nucleation details, avoiding mechanistic distortion and providing more comprehensive prediction and control capabilities for the synthesis of crystalline materials. Attached Figure Description

[0076] 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:

[0077] Figure 1This is a flowchart of the method for constructing nucleation reaction coordinates based on molecular simulation and genetic algorithm as described in this invention;

[0078] Figure 2 A schematic diagram of topological keys and binary encoding;

[0079] Figure 3 A schematic diagram of a path sampling algorithm without a target;

[0080] Figure 4 For the construction of the ice-water mixing model;

[0081] Figure 5 This is a graph showing the relationship between crystal nucleus size and the probability of a forward reaction. Detailed Implementation

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

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

[0084] The entire processing flow of this invention employs a hybrid C++ and Perl programming approach. C++ is used to handle the extensive calculations and definitions of local structures and the maximum likelihood estimation algorithm, while Perl is used to interface the genetic algorithm with other data stream processing. Therefore, by combining the computational efficiency of C++ with the flexibility of the Perl scripting language, simulation results and reaction coordinate iterations can be processed more efficiently in batches.

[0085] The method for constructing nucleation reaction coordinates based on molecular simulation and genetic algorithm described in this invention is illustrated in the flowchart below. Figure 1 As shown, it includes the following steps:

[0086] S1, Nucleation Transition State Sampling

[0087] The main approach involves using molecular simulation software to fully sample the configurations of the nucleation process and automatically adding labels to these configurations (i.e., obtaining a "configuration-label (crystal C / liquid L)" dataset). Since using empirical order parameters introduces errors into the "configuration-label" dataset, the reaction coordinates obtained after optimization using the initial dataset are used to replace the empirical order parameters, and the nucleation transition state sampling module is run again to obtain a more accurate "configuration-label" dataset.

[0088] S11. Simulate the nucleation process using molecular simulation software and obtain the simulation trajectory; then, first use empirical order parameters to obtain preliminary unlabeled configuration data, and then use empirical order parameters to preliminarily determine the crystalline phase C and the liquid phase L. After repeated iterations, obtain the configuration data of labeled C / L.

[0089] The specific method is as follows:

[0090] S111, Initial Sampling

[0091] Since the path sampling method requires a coarse initial (o) trajectory configuration set (connecting the two phases of CL), a coarse (which may be distorted) nucleation process is simulated here using strongly driven conditions or original dynamics. Then, the sampling method of S112 will automatically evolve this distorted trajectory to fit the real reaction process.

[0092] ① The nucleation process was simulated using high-driving conditions or meta-dynamics, and the simulated trajectory (o) was obtained.

[0093] ② Use order parameters (empirical order parameters are used for the first sampling to determine whether the two phases CL are connected) to initially describe the nucleation process (i.e., obtain preliminary unlabeled configuration data).

[0094] S112, Path sampling without a target

[0095] Since crystal nucleation is a rare event, conventional, brute-force molecular simulations are insufficient to obtain statistically significant configuration data. Therefore, this step further performs targetless shooting sampling on the initial, limited configuration data obtained in S111, thereby evolving a richer (statistically significant) configuration dataset. These configuration datasets are then automatically labeled. A schematic diagram of the targetless path sampling algorithm is shown below. Figure 3 As shown.

[0096] ① Randomly select a frame configuration in the trajectory (o), redistribute the velocities p of all particles to follow the Boltzmann distribution, and perform molecular dynamics simulations; use empirical order parameters to determine whether to enter the crystalline phase C or the liquid phase L;

[0097] ② Using the initial configuration in ①, assign velocity -p to all particles, perform molecular dynamics simulations, and use empirical order parameters to determine whether to enter the crystalline phase C or the liquid phase L;

[0098] ③ If the final configurations obtained through ① and ② fall into C and L respectively, then the newly generated trajectories are spliced ​​together to obtain a new trajectory (n), which replaces the original trajectory (o); steps ① and ② are executed again, and the initial configuration is saved and marked as either configuration C or L, each with a 50% probability, and the determination is made based on the generated random number;

[0099] ④ If the crystal phase C is entered through both ① and ②, then mark the initial configuration as C, save the initial configuration, and execute steps ① and ② again;

[0100] ⑤ If the liquid phase L has been entered through both ① and ②, then mark the initial configuration as L, save the initial configuration, and execute steps ① and ② again;

[0101] ⑥ Iterate repeatedly to obtain a configuration file with tagged C / L.

[0102] S12. Subsequently, local structure coding is constructed, and the empirical order parameters in S11 are replaced with local structure coding to obtain crystal nucleus information.

[0103] Local structure encoding involves processing the molecular structure in the Cartesian coordinate system into a binary code that can be processed and iterated by a computer through topological ring structure identification and topological bond structure definition. This allows for the identification of the crystalline phase C and the liquid phase L, and the acquisition of preliminary crystal nucleus information. Different codes correspond to different crystal nucleus information.

[0104] The specific method for local structure encoding is as follows:

[0105] I. Identification of Minimum Cycle Topology

[0106] ① Calculate and collect the index number (ID) of each particle's neighboring particles within the cutoff range (e.g., 0.35nm), and define these neighboring particles as the neighboring particles of the central particle, taking into account periodic boundary conditions;

[0107] ②Summarize the set of IDs representing the connections between the central particle and its neighboring particles; that is, {ID} i N, ID ij1 ID ij2 …ID ijn}, where N is ID i The total number of neighboring particles, ID ijn Indicates ID i The neighboring particles, n=1, 2, ...∈N;

[0108] ③ Utilize the set of connection relationship IDs of neighboring particles to obtain IDs using a recursive algorithm. i The ring topology information of the starting point; for example: ID i -ID j -ID k -ID g -ID i ;

[0109] ④ Eliminate invalid topological cycles to obtain the minimum cycle topological information;

[0110] For example, a 7-membered ring (1-2-3-4-5-7-1) is actually composed of a 4-membered ring (1-2-3-4-1) and a 5-membered ring (4-5-6-7-1-4). Therefore, the 7-membered ring is not a minimal topological ring structure and will be eliminated.

[0111] ⑤ Obtain the minimum topological ring structure of each particle and perform statistical analysis to obtain a new set of topological information, i.e., {ID, N}. 3-ring N 4-ring , ...N r-ring}, r≥3, where N r-ring This represents the number of r-rings that the ID particle is in.

[0112] For example, the set can be {ID, N} 3-ring N 4-ring N 5-ring N 6-ring N 7-ring}

[0113] II. Definition of Topological Bond Structure

[0114] ①Based on the set of IDs and topological information of neighboring particles in I, define the connection bond between two adjacent particles as B. i-j ;

[0115] ② Define key topology information, B i-j ={a,b,c,…}, where:

[0116] a=N 3-ring (i)·N 4-ring (j)-N 4-ring (i)·N 3-ring (j);

[0117] b=N 4-ring (i)·N 5-ring (j)-N 5-ring (i)·N 4-ring (j);

[0118] c=N 5-ring (i)·N 6-ring (j)-N 6-ring (i)·N 5-ring (j);

[0119] By analogy, information from multiple dimensions can be obtained;

[0120] For example, according to the above-mentioned set of topological information, including 10 dimensions, B i-j ={a,b,cd,e,f,g,h,n,m}, where:

[0121] a=N3-ring (i)·N 4-ring (j)-N 4-ring (i)·N 3-ring (j);

[0122] b=N 4-ring (i)·N 5-ring (j)-N 5-ring (i)·N 4-ring (j);

[0123] c=N 5-ring (i)·N 6-ring (j)-N 6-ring (i)·N 5-ring (j);

[0124] d=N 6-ring (i)·N 7-ring (j)-N 7-ring (i)·N 6-ring (j);

[0125] e=N 3-ring (i)·N 5-ring (j)-N 5-ring (i)·N 3-ring (j);

[0126] f=N 3-ring (i)·N 6-ring (j)-N 6-ring (i)·N 3-ring (j);

[0127] g=N 3-ring (i)·N 7-ring (j)-N 7-ring (i)·N 3-ring (j);

[0128] h=N 4-ring (i)·N 6-ring (j)-N 6-ring (i)·N 4-ring (j);

[0129] n=N 4-ring (i)·N 7-ring (j)-N 7-ring (i)·N 4-ring (j);

[0130] m=N 5-ring (i)·N 7-ring (j)-N 7-ring (i)·N 5-ring (j);

[0131] ③ The types of bond topological information during the statistical molecular simulation are identified and arranged from highest to lowest frequency. These are then randomly mapped to either "0" or "1", i.e., the mapping relationship f(B) i-j → 0 / 1;

[0132] Specifically, such as Figure 2 Examples of the mapping relationship between some topological keys and binary codes are given.

[0133] III. Classification of Crystal Nuclei and Liquid Particles

[0134] ① By utilizing the set of connection relationship IDs and topological information of neighboring particles, the bond topological information contained in each particle is obtained, and combined with the mapping relationship, the superimposed bond feature, i.e., B, is obtained. all =∑f(B i-j );

[0135] ②When particle B of i all When the value is ≥3, particle i is labeled as a crystalline particle; otherwise, it is labeled as a liquid-like particle.

[0136] ③ Connect all types of crystals according to the set of connection relationship IDs of neighboring particles, obtain the crystal nucleus clusters connected together, and obtain the largest cluster (with the most connected particles), that is, obtain the preliminary crystal nucleus information.

[0137] S2, Genetic Algorithm Iteration

[0138] The genetic algorithm and the configuration-label (C / L) dataset obtained by the sampling module are mainly used to perform self-evolutionary iteration of the nucleation reaction coordinates of binary basis ("local structure encoding module").

[0139] Using a genetic algorithm and the crystal nucleus information obtained from S1, the nucleation reaction coordinates of the binary base generated by the local structure encoding are self-evolved and iterated to obtain an accurate encoding form;

[0140] Then, replace the empirical sequence parameter with the accurate encoding form and execute S11 again to obtain the accurate label C / L configuration data.

[0141] The specific steps are as follows:

[0142] S21, Fitness Calculation

[0143] ① The configuration file of the tag C / L obtained by sampling the nucleation transition state, as well as the local structure code, are used to calculate the crystal nucleus information, which includes the crystal nucleus size CV1 and morphology CV2;

[0144] ② Construct reaction coordinates z(CV1, CV2, ...), which are linear combinations of CVi, i.e.: z = a0 + a1CV1 + a2CV2 + ... + a iCVi, where a i These are coefficients to be determined;

[0145] ③ For the configuration labeled C, construct the probability function:

[0146] P B =(1+erf[z(CV1, CV2,…)]) / 2;

[0147] ④ For the configuration labeled L, construct the probability function:

[0148] P A =1-(1+erf[z(CV1, CV2,…)]) / 2;

[0149] ⑤ Construct the maximum likelihood estimate by multiplying all probabilities together, that is:

[0150] L(z) = ∏P B (a i )·∏P A (a i );

[0151] ⑥ Take the logarithm of L and find the maximum value of log(L(z)), that is, find the maximum value of a. i The maximum value of the function log(L(z)) of the independent variable;

[0152] ⑦ Obtain a i The specific value of the constructed reaction coordinate z(CV1, CV2, ...) is the maximum value of log(L(z)) corresponding to the constructed reaction coordinate z(CV1, CV2, ...). The larger the fitness, the better the performance of the constructed reaction coordinate.

[0153] S22, Genetic Iteration

[0154] ① Perform batch operations on the binary encoding of random mapping in the local structure encoding, arrange the key topology information of different types from high to low according to the frequency of occurrence, and perform random mapping on each structure to obtain binary encoding, thus obtaining different binary encodings;

[0155] ②Calculate the fitness of each encoding separately;

[0156] ③ Randomly select 4 codes and compare their fitness, then pass on the two codes with the highest fitness to the next generation;

[0157] ④ For the selected codes, swap local segments in pairs;

[0158] For example, the quotation marks in 001“001”010 and 110“111”100 are intersected, and the resulting segments are: 001“111”010 and 110“001”100;

[0159] ⑤ Perform random mutations on the selected codes, that is, mutate "0" to "1" or "1" to "0";

[0160] ⑥ Repeat steps ②-⑤ until the binary code with the highest fitness is obtained (convergence);

[0161] S23, Sampling Iteration

[0162] ① Determine whether the new binary code obtained by genetic iteration converges with the previous generation's binary code, i.e., |log(L(z)| i+1 )-log(L(z) i )|<0.01;

[0163] ② If convergence is not achieved, use the newly iterated optimal reaction coordinates to perform S1 nucleation transition state sampling to obtain a more accurate labeled configuration dataset, and further perform S21 and S22 to iterate and obtain the final high-precision nucleation reaction coordinates;

[0164] ③ If convergence is achieved, execute the S3 crystal nucleus information output.

[0165] S3, Crystal nucleus information output

[0166] It is mainly used to encapsulate and visualize the optimized encoding, enabling more efficient batch processing of molecular simulation trajectory configurations.

[0167] S31. Encapsulate the iteratively optimized binary code;

[0168] S32. Input the binary code determined by the encapsulation into the local structure code, and output the crystal nucleus information, that is, the number of the crystal nucleus particles that make up the largest connectivity, and visualize the crystal nucleus with VMD or OVITIO.

[0169] S33. Calculate auxiliary thermodynamic and kinetic characteristics, including crystal nucleus sphericity, generalized diffusion coefficient of critical crystal nuclei, crystal nucleus adhesion rate, and crystal structure proportion.

[0170] The following explanation uses the self-evolving reaction coordinates of the ice nucleation process.

[0171] By constructing a liquid block of 4000 water molecules, such as Figure 4 As shown, in the process of liquid water freezing in a supercooled environment, this embodiment uses molecular dynamics to simulate the ice formation process. It develops self-evolving reaction coordinates for the ice nucleation process to accurately identify / describe the crystal nucleus structure in the solution (distinguishing between solid and liquid phases), thereby accurately exploring the ice crystal nucleation and growth phenomena in this process.

[0172] S1, Nucleation Transition State Sampling

[0173] S11. Simulate the nucleation process using molecular simulation software and obtain the simulation trajectory; then, first use empirical order parameters to obtain preliminary unlabeled configuration data, and then use empirical order parameters to preliminarily determine the crystalline phase C and the liquid phase L. After repeated iterations, obtain the configuration data of labeled C / L.

[0174] The specific method is as follows:

[0175] S111, Initial Sampling

[0176] ① The ice nucleation process was simulated using metadynamics in conjunction with q6 and q3 sequence parameters, and the simulated trajectory (o) was obtained;

[0177] ② Use q3 and q6 to preliminarily describe the nucleation process.

[0178] S112, Path sampling without a target

[0179] ① Randomly select a frame configuration in the trajectory (o), redistribute the velocities p of all particles to follow the Boltzmann distribution, and perform molecular dynamics simulations; use empirical order parameters q3 and q6 to determine whether to enter the crystalline phase C or the liquid phase L;

[0180] ② Using the initial configuration in ①, assign velocity -p to all particles, perform molecular dynamics simulations, and use empirical order parameters q3 and q6 to determine whether to enter the crystalline phase C or the liquid phase L;

[0181] ③ If the final configurations obtained through ① and ② fall into C and L respectively, then the newly generated trajectories are spliced ​​together to obtain a new trajectory (n), which replaces the original trajectory (o); steps ① and ② are executed again, and the initial configuration is saved and marked as either configuration C or L, each with a 50% probability, and the determination is made based on the generated random number;

[0182] ④ If the crystal phase C is entered through both ① and ②, then mark the initial configuration as C, save the initial configuration, and execute steps ① and ② again;

[0183] ⑤ If the liquid phase L has been entered through both ① and ②, then mark the initial configuration as L, save the initial configuration, and execute steps ① and ② again;

[0184] ⑥ Iterate repeatedly to obtain a configuration file with tagged C / L.

[0185] S12. Subsequently, local structure coding is constructed, and the empirical order parameters in S11 are replaced with local structure coding to obtain crystal nucleus information.

[0186] The specific method for local structure encoding is as follows:

[0187] I. Identification of Minimum Cycle Topology

[0188] ① Calculate and collect the ID of each water molecule within the cutoff range (0.35nm, angle OOH<35°) of its neighboring water molecules, taking into account periodic boundary conditions;

[0189] ②Summarize the set of IDs representing the connections between the central water molecule and its neighboring water molecules; that is, {ID} i N, ID ij1 ID ij2 …ID ijn}, where N is ID i The total number of neighboring water molecules, ID ijn Indicates ID i The neighboring water molecules, n=1, 2, ...∈N;

[0190] ③ Utilize the set of IDs representing the connections between neighboring water molecules, and use a recursive algorithm to obtain the ID. i The ring topology information of the starting point; for example: ID i -ID j -ID k -ID g -ID i ;

[0191] ④ Eliminate invalid topological loops to obtain the minimum water loop topology information;

[0192] For example, a 7-membered ring (1-2-3-4-5-7-1) is actually composed of a 4-membered ring (1-2-3-4-1) and a 5-membered ring (4-5-6-7-1-4). Therefore, the 7-membered ring is not a minimal topological ring structure and will be eliminated.

[0193] A certain 7-membered ring can also be formed by (1-2-3-1) and (3-4-5-6-7-1-3), and this 7-membered ring will also be eliminated;

[0194] ⑤ Obtain the smallest topological ring structure containing each water molecule, and perform statistical analysis to obtain a new set of topological information, i.e., {ID, N}. 3-ring N 4-ring N 5-ring N 6-ring N 7-ring};

[0195] II. Definition of Topological Bond Structure

[0196] ①Based on the set of IDs and topological information of neighboring particles in I, define the connection bond between two adjacent particles as B. i-j ;

[0197] ② Define key topology information, including 10 dimensions, B i-j ={a,b,cd,e,f,g,h,n,m}, where:

[0198] a=N 3-ring (i)·N 4-ring (j)-N 4-ring (i)·N 3-ring (j);

[0199] b=N 4-ring (i)·N 5-ring (j)-N 5-ring (i)·N 4-ring (j);

[0200] c=N 5-ring (i)·N 6-ring (j)-N 6-ring (i)·N 5-ring (j);

[0201] d=N 6-ring (i)·N 7-ring (j)-N 7-ring (i)·N 6-ring (j);

[0202] e=N 3-ring (i)·N 5-ring (j)-N 5-ring (i)·N 3-ring (j);

[0203] f=N 3-ring (i)·N 6-ring (j)-N 6-ring (i)·N 3-ring (j);

[0204] g=N 3-ring (i)·N 7-ring (j)-N 7-ring (i)·N 3-ring (j);

[0205] h=N 4-ring (i)·N 6-ring (j)-N 6-ring (i)·N 4-ring (j);

[0206] n=N 4-ring (i)·N 7-ring (j)-N 7-ring (i)·N 4-ring (j);

[0207] m=N 5-ring (i)·N 7-ring (j)-N 7-ring (i)·N 5-ring (j);

[0208] ③ The types of bond topological information during the statistical molecular simulation are identified and arranged from highest to lowest frequency. These are then randomly mapped to either "0" or "1", i.e., f(B i-j → 0 / 1;

[0209] For example: (1,0,2,12,0,0,1,0,9,3)→0;

[0210] (0,10,3,7,0,1,1,0,0,0)→1;

[0211] (12,12,0,0,0,1,1,0,9,0)→1;

[0212] (0,0,0,0,0,12,12,0,0,0)→0; etc.

[0213] Therefore, by arranging these 10-dimensional tensors in a fixed order, a corresponding binary encoding mapping can be obtained: 0, 1, 1, 0, ...

[0214] III. Classification of Crystal Nuclei and Liquid Particles

[0215] ① By utilizing the set of connection relationship IDs and topological information of neighboring particles, the bond topological information contained in each water molecule is obtained. Combined with the mapping relationship, the superposition bond characteristics, i.e., B, are obtained. all =∑f(B i-j );

[0216] ②When the B of water molecules i all When the value is ≥3, the i-water molecules are labeled as crystalline particles; otherwise, they are labeled as liquid particles.

[0217] ③ Connect all types of crystals according to the ID set of the connection relationship between nearby water molecules, obtain the crystal nucleus clusters connected together, and obtain the largest cluster (with the most connected particles), that is, obtain the preliminary crystal nucleus information.

[0218] S2, Genetic Algorithm Iteration

[0219] S21, Fitness Calculation

[0220] ①Use the large number of labeled (C / L) configuration files obtained by the ice nucleation transition state sampling module, and the ice crystal nucleus identification algorithm defined in the local structure encoding module to calculate and obtain crystal nucleus information (crystal nucleus size CV1, sphericity CV2).

[0221] ② Construct reaction coordinates z(CV1, CV2, ...), which are linear combinations of CVi, i.e.: z = a0 + a1CV1 + a2CV2 + ... + a i CVi, where a iThese are coefficients to be determined;

[0222] ③ For the configuration labeled C, construct the probability function:

[0223] P B =(1+erf[z(CV1, CV2,…)]) / 2;

[0224] ④ For the configuration labeled L, construct the probability function:

[0225] P A =1-(1+erf[z(CV1, CV2,…)]) / 2;

[0226] ⑤ Construct the maximum likelihood estimate by multiplying all probabilities together, that is:

[0227] L(z) = ∏P B (a i )·∏P A (a i );

[0228] ⑥ Take the logarithm of L and find the maximum value of log(L(z)), that is, find the maximum value of a. i The maximum value of the function log(L(z)) of the independent variable;

[0229] ⑦ Obtain a i The specific value of the constructed reaction coordinate z(CV1, CV2, ...) is the maximum value of log(L(z)) corresponding to the constructed reaction coordinate z(CV1, CV2, ...). The larger the fitness, the better the performance of the constructed reaction coordinate.

[0230] S22, Genetic Iteration

[0231] ① Perform batch operations on the binary encoding of random mappings in the local structure encoding, arrange the key topology information of different types in descending order of frequency, and perform random mapping on each structure to obtain binary encoding, resulting in 80 different binary encodings, such as: (1,1,1,1,0,0,1,0,…), (1,0,1,1,1,0,1,0,…), (0,1,1,1,0,0,1,0,…), etc.

[0232] ②Calculate the fitness of each encoding separately;

[0233] ③ Randomly select 4 codes and compare their fitness, then pass on the two codes with the highest fitness to the next generation;

[0234] ④ For the selected codes, swap local segments in pairs;

[0235] For example, the quotation marks in 001“001”010 and 110“111”100 are intersected, and the resulting segments are: 001“111”010 and 110“001”100;

[0236] ⑤ Perform random mutations on the selected codes, that is, mutate "0" to "1" or "1" to "0";

[0237] ⑥ Repeat steps ②-⑤ until the binary code with the highest fitness is obtained (convergence);

[0238] S23, Sampling Iteration

[0239] ① Determine whether the new binary code obtained by genetic iteration converges with the previous generation's binary code, i.e., |log(L(z)| i+1 )-log(L(z) i )|<0.01;

[0240] ② If convergence is not achieved, use the newly iterated optimal reaction coordinates to perform S1 nucleation transition state sampling to obtain a more accurate labeled configuration dataset, and further perform S21 and S22 to iterate and obtain the final high-precision nucleation reaction coordinates;

[0241] ③ If convergence is achieved, execute the S3 crystal nucleus information output.

[0242] S3, Crystal nucleus information output

[0243] S31. Encapsulate the iteratively optimized binary code;

[0244] S32. Input the binary code determined by the encapsulation into the local structure code, and output the crystal nucleus information, that is, the number of the crystal nucleus particles that make up the largest connectivity, and visualize the crystal nucleus with VMD or OVITIO.

[0245] S33, Calculate auxiliary thermodynamic and kinetic characteristics.

[0246] Figure 4 By using optimized ice nucleation reaction coordinates, ice nuclei can be identified and extracted from liquid phase configurations that are difficult to discern with the naked eye. The following algorithm was used to evaluate the performance of the reaction coordinates:

[0247] ① Collect a large number of different configurations ( Figure 4 The CCP has collected over 100 configurations.

[0248] ② For each configuration, the nucleus size information is extracted using reaction coordinates (using the traditional CHILL algorithm and the optimized ice nucleation reaction coordinates, respectively).

[0249] ③ Perform 100+ molecular dynamics simulations for each configuration to obtain the probability of a positive reaction.

[0250] ④ Drawing Figure 5 The relationship between crystal nucleus size and the probability of forward reaction is shown in the figure.

[0251] from Figure 5 It can be seen that when the reaction coordinates and the forward reaction probability of the method of this invention fit well (i.e., they exhibit a clear positive correlation, and the larger the slope), it indicates that the selected nucleation reaction has excellent performance and can accurately describe the microscopic details of the nucleation process. Conversely, compared with the traditional CHILL algorithm, the defined nucleus size and the forward reaction probability are almost orthogonal, indicating a poor correlation with the forward reaction of ice nucleation, and its performance is far inferior to the reaction coordinates optimized by this invention.

[0252] 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 constructing nuclear reaction coordinates based on molecular simulation and genetic algorithm, characterized by, The method comprises the following steps: S1, first, nucleation transition state sampling is performed; S11, a molecular simulation software is used to simulate the nucleation process, and a simulation trajectory is obtained; then, an empirical order parameter is used to obtain preliminary unlabeled configuration data, and the crystal phase C and the liquid phase L are preliminarily determined, and through repeated iteration, configuration data with labels C / L is obtained; S12, then, local structure coding is constructed, and the local structure coding is used to replace the empirical order parameter in S11 to obtain crystal nucleus information; The local structure coding is processed into binary coding which can be iterated by a computer by topological ring structure identification and topological bond structure definition, so that the crystal phase C and the liquid phase L are identified, and preliminary crystal nucleus information is obtained; different coding corresponds to different crystal nucleus information; S2, secondly, a genetic algorithm iteration is performed; The nucleation reaction coordinates of the binary base generated by the local structure coding are iterated by using the genetic algorithm and the crystal nucleus information obtained in S1, so that accurate coding forms are obtained; Then, the accurate coding forms are used to replace the empirical order parameter to perform S11 again, and accurate configuration data with labels C / L is obtained; S3, finally, crystal nucleus information is output; The optimized coding is packaged and visualized, and a batch of molecular simulation trajectory configurations are processed.

2. The method of constructing nuclear reaction coordinates based on molecular simulation and genetic algorithm according to claim 1, wherein, The specific method of S11 is as follows: S111, first, initial sampling is performed, that is, a molecular simulation software is used to simulate the nucleation process, and a simulation trajectory is obtained; whether the two phases C-L are connected is determined by using an empirical order parameter, the nucleation process is preliminarily described, and preliminary unlabeled configuration data is obtained; S112, secondly, path sampling without target is adopted; ①a frame of configuration is randomly selected in the trajectory, all particle velocities p are re-distributed and subjected to a Boltzmann distribution, and molecular dynamics simulation is performed; whether the crystal phase C or the liquid phase L is entered is determined by using an empirical order parameter; ②the initial configuration in ① is used, velocities-p of all particles are distributed, molecular dynamics simulation is performed, and whether the crystal phase C or the liquid phase L is entered is determined by using an empirical order parameter; ③if the final configurations in ① and ② fall into C and L respectively, the newly generated trajectory is spliced to obtain a new trajectory, and the original trajectory is replaced; steps ① and ② are executed again, and the initial configuration is saved and marked as a C or L configuration, each having a probability of 50%, and the determination is made according to a generated random number; ④if ① and ② both enter the crystal phase C, the initial configuration is marked as C, and the initial configuration is saved; steps ① and ② are executed again; ⑤if ① and ② both enter the liquid phase L, the initial configuration is marked as L, and the initial configuration is saved; steps ① and ② are executed again; ⑥repeated iteration is performed to obtain a configuration file with labels C / L.

3. The method of constructing nuclear reaction coordinates based on molecular simulation and genetic algorithm according to claim 1, wherein, The specific method of the local structure coding is as follows: Ⅰ, first, minimum ring topological structure identification is performed; ①the index serial number corresponding to the adjacent particles of each particle within a cutoff range is calculated and collected, which is referred to as ID, and the adjacent particles are defined as the neighboring particles of the center particle; ② The connection relationship ID set of the center particle and its adjacent particles, that is, {ID i , N, ID ij1 , ID ij2 , …, ID ijn}, wherein N is the total number of adjacent particles of ID i , ID ijn represents the adjacent particle of ID i , n=1, 2, … ∈N; ③Using the connection relationship ID set of the adjacent particles, using the recursive algorithm to obtain ID i Ring topology information with the starting point; ④invalid topological rings are eliminated to obtain minimum ring topological information; ⑤ Obtain the minimum topological ring structure where each particle is located, and perform statistics to obtain a new topological information set, namely {ID, N 3-ring , N 4-ring , …N r-ring}, r≥3, wherein N r-ring represents the number of r-element rings through the ID particle is located; Ⅱ, then, topological bond structure is defined; ① According to the ID set and the topological information set of the connection relationship of adjacent particles in I, the connection bond between the two adjacent particles is defined B i-j ; (ii) defining key topological information, B i-j = {a, b, c,...}, where: a = N 3-ring (i) • N 4-ring (j) - N 4-ring (i) • N 3-ring (j); b = N 4-ring (i) • N 5-ring (j) - N 5-ring (i) • N 4-ring (j); c=N 5-ring (i)·N 6-ring (j)-N 6-ring (i)·N 5-ring (j); By analogy, multiple dimensional information is obtained; iii. Counting the types of bond topological information in the molecular simulation process, and arranging them in descending order of frequency of occurrence, and randomly mapping them as "0" or "1", that is, the mapping relationship f(B i-j )→0 / 1; Ⅲ, finally, crystal nucleus-like and liquid-like particles are classified; all =∑f(B i-j );​ (ii) when the B all ≥ 3, the i-particle is labeled as a quasi-crystalline particle, and vice versa. ③Connect all crystal particles according to the ID set of the connection relationship of adjacent particles, obtain the crystal nucleus cluster connected in one piece, and obtain the largest cluster, that is, obtain the preliminary crystal nucleus information.

4. The method of constructing nuclear reaction coordinates based on molecular simulation and genetic algorithm according to claim 3, wherein, In I, the ID of adjacent particles is also considered in the periodic boundary condition; the cutoff range is ≤0.35 nm.

5. The method of constructing nuclear reaction coordinates based on molecular simulation and genetic algorithm according to claim 1, wherein, The specific steps of S2 are: S21, first, fitness calculation is performed; ①Use the configuration file with labels C / L obtained in the nucleation transition state sampling and the local structure encoding to calculate and obtain the crystal nucleus information, including the crystal nucleus size CV1 and the morphology CV2; ii) constructing reaction coordinates z (CV1, CV2,...), which are linear combinations of CVi, i.e. z = a0+ a1CV1+ a2CV2+... a i CVi, where a i are undetermined coefficients; ③For the configurations labeled as C, construct the probability function: P B = (1 + erf[z(CV1, CV2,...)] ) / 2; ④For the configurations labeled as L, construct the probability function: P A =1-(1+erf[z(CV1, CV2,...)] / 2; ⑤Construct the maximum likelihood estimation to multiply all the probabilities, that is: L(z) = ΠP B (a i )·ΠP A (a i ) (6) Taking the logarithm of L, solve the maximum value of log(L(z)), that is, solve the maximum value of the function log(L(z)) with a i as the independent variable; Obtained a i The specific value of the maximum value of the log(L(z)) corresponding to the reaction coordinate z(CV1, CV2, …) is the fitness of the reaction coordinate; the greater the fitness, the better the performance of the constructed reaction coordinate; S22, then, genetic iteration is performed; ①Perform batch operations on the binary encoding of the random mapping in the local structure encoding, arrange the key topological information of different types according to the frequency from high to low, and obtain different binary encodings by randomly mapping each type of structure; ②Use each encoding to calculate their fitness respectively; ③Randomly select 4 encodings to compare their fitness, and pass the two encodings with the highest fitness to the next generation; ④Exchange local fragments between the selected encodings; ⑤Randomly mutate the selected encodings, that is, change "0" to "1" or "1" to "0"; ⑥Repeat ②-⑤ until the binary encoding with the highest fitness is obtained; S23, finally, sampling iteration is performed; i+1 whether the new binary code obtained by genetic iteration converges with the binary code of the last generation, i.e., |log(L(z i )) - log(L(z i ))| < 0.01; ②If it does not converge, use the best reaction coordinates of the new iteration to perform S1 nucleation transition state sampling to obtain a more accurate labeled configuration dataset, and further perform S21 and S22 to obtain the final high-precision nucleation reaction coordinates; ③If it converges, perform S3 crystal nucleus information output.

6. The method of constructing nuclear reaction coordinates based on molecular simulation and genetic algorithm according to claim 1, wherein, The specific steps of S3 are: S31, package the binary encoding optimized by iteration; S32, import the packaged binary encoding into the local structure encoding, and output the crystal nucleus information, that is, the number of the largest connected crystal nucleus particles, and visualize the crystal nucleus with VMD or OVITIO; S33, calculate auxiliary thermodynamic and kinetic characteristics.

7. The method of constructing nuclear reaction coordinates based on molecular simulation and genetic algorithm according to claim 6, wherein, The auxiliary thermodynamic and kinetic characteristics include the crystal nucleus sphericity, the generalized diffusion coefficient of the critical crystal nucleus, the adhesion rate of the crystal nucleus, and the proportion of the crystal structure.

8. The method of constructing nuclear reaction coordinates based on molecular simulation and genetic algorithm according to claim 1, wherein, S11 uses high driving conditions or meta-dynamics simulation to nucleate the process.

Images