阅读说明：本技术 基于几何边界运算的分子动力学边界条件快速施加方法 (Method for quickly applying molecular dynamics boundary conditions based on geometric boundary operation ) 是由 黄志刚 葛露明 黄亚军 何斌 覃圣超 史纪童 于 2021-01-29 设计创作，主要内容包括：本发明公开了一种基于几何边界运算的分子动力学边界条件快速施加方法,包括：根据分子动力学仿真模型,提取其几何边界,通过几何边界运算找出边界粒子,获取边界域内每个边界粒子的位置信息；为每个边界粒子赋予合理的速度和作用力,使得边界粒子的宏观物理量与微观状态相一致,从而得到完整的分子动力学边界条件。本发明给出了在复杂分子动力学系统上快速施加边界条件的一般方案；施加了平滑过渡的速度边界条件和应力边界条件,边界粒子的速度分布和势能分布,与系统局部的微观动力学状态一致,将边界附近的状态波动限制在较低的水平,解决施加边界条件任务繁重的问题,同时极大减少边界振荡。(The invention discloses a method for quickly applying a molecular dynamics boundary condition based on geometric boundary operation, which comprises the following steps: extracting a geometric boundary of the molecular dynamics simulation model according to the molecular dynamics simulation model, finding out boundary particles through geometric boundary operation, and acquiring position information of each boundary particle in a boundary domain; and (3) endowing each boundary particle with reasonable speed and force, so that the macroscopic physical quantity of the boundary particle is consistent with the microscopic state, and the complete molecular dynamics boundary condition is obtained. The invention provides a general scheme for rapidly applying boundary conditions on a complex molecular dynamics system; the velocity boundary condition and the stress boundary condition of smooth transition are applied, the velocity distribution and the potential energy distribution of boundary particles are consistent with the local micro dynamic state of the system, the state fluctuation near the boundary is limited to a lower level, the problem of heavy task of applying the boundary condition is solved, and the boundary oscillation is greatly reduced.)

1. A molecular dynamics boundary condition rapid application method based on geometric boundary operation is characterized by comprising the following steps:

extracting a geometric boundary of the molecular dynamics simulation model according to the molecular dynamics simulation model, finding out boundary particles through geometric boundary operation, and acquiring position information of each boundary particle in a boundary domain;

and (3) endowing each boundary particle with reasonable speed and force, so that the macroscopic physical quantity of the boundary particle is consistent with the microscopic state, and the complete molecular dynamics boundary condition is obtained.

2. The method for rapidly applying boundary conditions of molecular dynamics based on geometric boundary operations according to claim 1, wherein the finding of boundary particles through geometric boundary operations to obtain the position information of each boundary particle in the boundary domain comprises:

constructing a solid model of the geometric boundary according to the geometric characteristics of the boundary;

acquiring the geometric information of the boundary abstraction according to the entity model, and performing discrete operation on the entity model according to the geometric information of the boundary model to obtain discrete particles; and judging the position relation between the discrete particles and the boundary by adopting an intersection method according to the geometric characteristics of the boundary so as to determine the boundary particles and the simulation domain particles and finally obtain the complete position information of the boundary particles.

3. The method for rapidly applying the molecular dynamics boundary condition based on the geometric boundary operation of claim 2, wherein in the process of constructing the solid model of the geometric boundary according to the geometric features of the boundary, the more regular geometric boundary is obtained by a bias method, and the more complex boundary is obtained by a Boolean operation.

4. The method for rapidly applying molecular dynamics boundary conditions based on geometric boundary operations according to claim 1, wherein the step of endowing each boundary particle with reasonable speed comprises the following steps:

for all the boundary particles, based on a speed correction mode, correcting the speed of the boundary particles by using a relaxation factor, so that the speed of the boundary particles is slowly changed towards the direction of an analog domain, wherein a speed correction formula is as follows:

V′_i＝V_i+h(x_i)(u-V_i) Formula 1

In formula 1, i represents the ith boundary particle; v'_iThe corrected speed of the boundary particles; v_iIs the initial velocity of the boundary particle; h (x)_i) Is a relaxation factor with a value of 0 to h₀Change within the interval; x is the number of_iA certain coordinate value which is changed from 0 to A along the normal direction of the analog domain for the particle; u is the expected speed of the boundary particles in the target state of the molecular dynamics simulation system, and A is the width of the boundary domain along the direction of the normal vector n.

5. The method of claim 4, wherein a normal vector n of the boundary domain facing the simulation domain is determined based on the geometric features of the boundary, and the velocity is sequentially assigned to all the boundary particles along the normal vector direction by the numerical calculation of the velocity correction formula according to the position information of the boundary particles.

6. The method for rapidly applying molecular dynamics boundary conditions based on geometric boundary operations according to claim 1, wherein the step of giving reasonable acting force to each boundary particle comprises the following steps:

in practical simulations, to apply boundary forces within the boundary domain, it is assumed that there are some external particles outside the boundary domain that exert an action potential on the boundary particle, the external particles being at the cutoff radius R_CEffective action potentials exist on the boundary particles within the range, and the sum of the action potentials forms effective boundary action force; assuming that any boundary particle is surrounded by a shell outside the boundary domain, the force of the single shell on any boundary particle within the boundary domain is first determined, followed by the total force exerted by the outer particle on any boundary particle within the boundary domain.

7. The method for rapidly applying molecular dynamics boundary conditions based on geometric boundary operations of claim 6, wherein the determining the acting force of the monolayer thin shell on any boundary particle in the boundary domain comprises:

the number of particles in the shell can be expressed as: rho g (r) S (r) dr, the acting force of the particles in the shell on the boundary particles exists, and along the normal direction of the simulation domain, the acting force of the particles in any thin shell on any particle in the boundary domain is as follows: - ρ g (r) T (r) V (r) dr; wherein rho is the particle density in the molecular dynamics simulation system; g (r) is a function of the radial distribution of the particles; s (r) is the area of the shell; t (r) is the projection of the shell area on the outer layer of the boundary domain; v (r) is the action potential of a single particle in the shell on any particle in the boundary domain; dr is the thickness of the shell.

8. The method for rapidly applying molecular dynamics boundary conditions based on geometric boundary operations of claim 7, wherein the total force applied by the external particle to any boundary particle in the boundary domain is expressed as:

in the above formula, Φ is the stress tensor; p is Viry pressure; t (r, x) is the projection of the thin shell area on the outer layer of the boundary domain when the boundary particle has the coordinate x.

9. The method of claim 8, wherein the boundary particle i is truncated at its truncation radius R_cWithin the range, the acting force of any external particle j outside the boundary area on the external particle j is F_ijThe effective acting force F applied to the boundary particle i by all the particles within the truncation radius range outside the boundary domain is calculated by using the formula 2 by assuming that the periphery of any boundary particle outside the boundary domain is wrapped by a thin shell_mThe final realization is to apply the appropriate boundary forces for all boundary particles.

Technical Field

The invention relates to the technical field of molecular dynamics computational simulation, in particular to a molecular dynamics boundary condition rapid application method based on geometric boundary operation.

Background

With the development of computer technology and numerical simulation theory, Molecular Dynamics (MD) methods are more and more concerned. The molecular dynamics is a numerical simulation method which obtains the phase trajectory of the system by numerically solving the mechanical motion equation of the molecular system and counts the structural characteristics and properties of the system. The method is widely applied to biological calculation, microfluid and micro-nano engineering simulation.

Molecular dynamics is a particle-based computational method, and the application of boundary conditions is a very difficult task. In a traditional grid-based computing method, such as a finite element method, boundary conditions are applied to a geometric boundary, and the traditional grid-based computing method can be conveniently operated in a graphical interaction mode and the like. However, molecular dynamics, as a discrete particle system, must impose boundary effects on all boundary particles. The selection and manipulation of a large number of boundary particles can be a very burdensome task. Especially for the computation region with complex shape, the fast operation of boundary particles still lacks an effective general processing method.

In addition, in a molecular dynamics system, the motion states of particles, including position, velocity, and potential energy, form a specific distribution. The motion state of the boundary particles must be adapted to the local distribution function, and the boundary effect cannot be simply averaged to the boundary particles, otherwise oscillation of the system state is caused, calculation error is increased, and even calculation non-convergence is caused.

At present, the molecular dynamics calculation scale is larger and larger, the calculation region shape is more and more complex, and the calculation requirements on unsteady and non-uniform systems are more and more. How to use a visualization method to quickly establish a calculation model and apply corresponding molecular dynamics boundary conditions is an important challenge facing modern molecular dynamics simulation. The following difficulties mainly need to be solved:

(1) the fast selection and manipulation of boundary particles, which can be as large as millions for large scale molecular dynamics calculations, is a very burdensome task.

(2) The proper speed and potential energy are given to each boundary particle, the action among molecular dynamic particles is very strong, if the motion state of the particles does not correspond to the distribution function, violent oscillation is possible to occur, instantaneous deviation can occur between the local condition and the target state, and boundary effects such as density fluctuation and the like are easy to occur.

Disclosure of Invention

The invention aims to provide a method for quickly applying a molecular dynamics boundary condition based on geometric boundary operation, which is used for solving the problem of heavy task of applying the boundary condition in the prior art.

In order to realize the task, the invention adopts the following technical scheme:

a molecular dynamics boundary condition rapid application method based on geometric boundary operation comprises the following steps:

extracting a geometric boundary of the molecular dynamics simulation model according to the molecular dynamics simulation model, finding out boundary particles through geometric boundary operation, and acquiring position information of each boundary particle in a boundary domain;

and (3) endowing each boundary particle with reasonable speed and force, so that the macroscopic physical quantity of the boundary particle is consistent with the microscopic state, and the complete molecular dynamics boundary condition is obtained.

Further, the finding out the boundary particle through the geometric boundary operation to obtain the position information of each boundary particle in the boundary domain includes:

constructing a solid model of the geometric boundary according to the geometric characteristics of the boundary;

acquiring the geometric information of the boundary abstraction according to the entity model, and performing discrete operation on the entity model according to the geometric information of the boundary model to obtain discrete particles; and judging the position relation between the discrete particles and the boundary by adopting an intersection method according to the geometric characteristics of the boundary so as to determine the boundary particles and the simulation domain particles and finally obtain the complete position information of the boundary particles.

Further, in the process of constructing the solid model of the geometric boundary according to the geometric features of the boundary, a bias method is adopted for obtaining the more regular geometric boundary, and Boolean operation is used for obtaining the more complex boundary.

Further, the endowing of each boundary particle with a reasonable speed includes:

for all the boundary particles, based on a speed correction mode, correcting the speed of the boundary particles by using a relaxation factor, so that the speed of the boundary particles is slowly changed towards the direction of an analog domain, wherein a speed correction formula is as follows:

V_i'＝V_i+h(x_i)(u-V_i) Formula 1

In formula 1, i represents the ith boundary particle; v_i' is the corrected velocity of the boundary particle; v_iIs the initial velocity of the boundary particle; h (x)_i) Is a relaxation factor with a value of 0 to h₀Change within the interval; x is the number of_iA certain coordinate value which is changed from 0 to A along the normal direction of the analog domain for the particle; u is the expected speed of the boundary particles in the target state of the molecular dynamics simulation system, and A is the width of the boundary domain along the direction of the normal vector n.

Furthermore, a normal vector n of the boundary domain towards the direction of the simulation domain is determined based on the geometric characteristics of the boundary, and the speed is sequentially given to all the boundary particles along the direction of the normal vector according to the position information of the boundary particles in a numerical calculation mode of the speed correction formula.

Further, the imparting of reasonable force to each boundary particle includes:

in practical simulations, to apply boundary forces within the boundary domain, it is assumed that there are some external particles outside the boundary domain that exert an action potential on the boundary particle, the external particles being at the cutoff radius R_CEffective action potentials exist on the boundary particles within the range, and the sum of the action potentials forms effective boundary action force; assuming that any boundary particle is surrounded by a shell outside the boundary domain, the force of the single shell on any boundary particle within the boundary domain is first determined, followed by the total force exerted by the outer particle on any boundary particle within the boundary domain.

Further, the determining the acting force of the monolayer thin shell on any boundary particle in the boundary domain comprises:

the number of particles in the shell can be expressed as: rho g (r) S (r) dr, the acting force of the particles in the shell on the boundary particles exists, and along the normal direction of the simulation domain, the acting force of the particles in any thin shell on any particle in the boundary domain is as follows:wherein rho is the particle density in the molecular dynamics simulation system;

g (r) is a function of the radial distribution of the particles; s (r) is the area of the shell; t (r) is the projection of the shell area on the outer layer of the boundary domain; v (r) is the action potential of a single particle in the shell on any particle in the boundary domain; dr is the thickness of the shell.

Further, the total force applied by the external particle to any boundary particle within the boundary domain is expressed as:

in the above formula, Φ is the stress tensor; p is Viry pressure; t (r, x) is the projection of the thin shell area on the outer layer of the boundary domain when the boundary particle has the coordinate x.

Further, for the boundary particle i, at its cutoff radius R_cWithin the range, the acting force of any external particle j outside the boundary area on the external particle j is F_ijThe effective acting force F applied to the boundary particle i by all the particles within the truncation radius range outside the boundary domain is calculated by using the formula 2 by assuming that the periphery of any boundary particle outside the boundary domain is wrapped by a thin shell_mThe final realization is to apply the appropriate boundary forces for all boundary particles.

Compared with the prior art, the invention has the following technical characteristics:

1. the method can quickly find out the boundary particles and apply reasonable boundary conditions, solves the problem of heavy task of applying the boundary conditions, and greatly reduces boundary oscillation.

2. The invention provides a general scheme for rapidly applying boundary conditions on a complex molecular dynamics system; and applying a speed boundary condition and a stress boundary condition of smooth transition, and the speed distribution and the potential energy distribution of boundary particles are consistent with the local micro dynamic state of the system, so that the state fluctuation near the boundary is limited to a lower level.

Drawings

FIG. 1 is a schematic flow diagram of the process of the present invention;

FIG. 2 is a schematic diagram of the boundary particle acquisition principle of the present invention;

FIG. 3 is a schematic diagram of the velocity imparted to boundary particles by the present invention;

FIG. 4 is a schematic diagram of the present invention applied boundary particle forces.

Detailed Description

Referring to fig. 1, the invention provides a method for rapidly applying a molecular dynamics boundary condition based on geometric boundary operation, which first finds boundary particles rapidly through geometric boundary operation on a molecular dynamics simulation model, and obtains position information of each boundary particle in a boundary domain; and secondly, endowing each boundary particle with reasonable speed and acting force, so that the macroscopic physical quantity of the boundary particle is consistent with the microscopic state, and the complete molecular dynamics boundary condition is obtained.

Specifically, the method for rapidly applying the molecular dynamics boundary condition of the invention comprises three processes, which are specifically described as follows:

1. boundary particle acquisition from given geometric boundaries

The step is the most basic step when boundary conditions are applied, and mainly obtains molecular dynamics boundary particle information through dispersion according to geometric boundaries in a molecular dynamics simulation model to construct a boundary model. The method comprises the following specific steps:

1.1 constructing solid models of geometric boundaries

Extracting a geometric boundary of the molecular dynamics simulation model according to the molecular dynamics simulation model, and constructing an entity model of the geometric boundary according to the geometric characteristics of the boundary; in the process, a bias method can be adopted to obtain a more regular geometric boundary, and a Boolean operation can be used to obtain a more complex boundary; wherein, for a solid model, a straight line in a two-dimensional model or a plane in a three-dimensional model can be considered as a more regular geometric boundary.

As shown in fig. 2 (a), in the example given in fig. 2, for the geometric features of the more regular molecular dynamics model boundary, the bias method may be directly adopted to bias the model b along the X direction and a along the Y direction, so as to obtain the boundary solid model; as shown in fig. 2 (b), when the boundary of the simulation model is irregular or is a complex curved surface, the boundary model may be obtained by performing boolean operations on two or more models.

1.2 building a geometric boundary particle model

After the step 1.1 is completed, acquiring the abstract geometric information of the boundary according to the entity model of the geometric boundary, and performing discrete operation on the entity model according to the geometric information of the boundary model; after the discrete particles are obtained, the position relation between the discrete particles and the boundary is judged by adopting an intersection method according to the geometric characteristics of the boundary, so that the boundary particles and the particles in the analog domain are determined, as shown in fig. 2, the position information of the complete boundary particles is finally obtained, and a boundary particle model is constructed.

2. Imparting reasonable velocities to boundary particles

The boundary conditions applied to molecular dynamics simulation systems are largely divided into two parts, velocity and force. In the step, the speed is given to the boundary particles, the speed needs to meet the Dirichlet boundary condition, and meanwhile, the microscopic state of the molecular dynamics simulation system particles is consistent with the macroscopic physical quantity. In order to avoid the large deviation of the particle speed in the boundary domain and the simulated domain, which causes the system oscillation, the particle speed in the boundary domain must be changed slowly, i.e. the speed is in 'soft transition'; the specific implementation manner of the step is as follows:

for all boundary particles, based on a speed correction mode, a relaxation factor in the form of a "bell-shaped function" is adopted to correct the speed of the boundary particles, so that the speed of the boundary particles changes slowly towards the direction of the analog domain, and the speed correction formula is as follows:

V_i'＝V_i+h(x_i)(u-V_i) Formula 1

In the formula 1, the reaction mixture is,i represents the ith boundary particle; v_i' is the corrected velocity of the boundary particle; v_iIs the initial velocity of the boundary particle; h (x)_i) Is a relaxation factor with a value of 0 to h₀Within a range of variation, h₀The specific value of (a) depends on the specific molecular dynamics simulation model, e.g. h in general₀The value of (d) can be 1; x is the number of_iA certain coordinate value which is the change of the particle from 0 to A (the width of the boundary domain along the direction of the normal vector n) along the normal direction of the analog domain; u is the expected velocity of the boundary particle in the molecular dynamics simulation system target state.

The boundary particles start from the boundary, the speed value towards the normal vector n direction of the simulation domain is slowly changed from u, the introduced relaxation factor function is required to satisfy h '(0) ═ 0, h' (A) ═ 0, and the relaxation factor is tightly supported in the boundary domain, so that the boundary action is limited in the boundary domain, and the introduction of the relaxation factor function depends on a specific molecular dynamics simulation model.

As shown in fig. 3, a normal vector n of the boundary domain toward the simulation domain is determined based on the geometric features of the boundary, and based on the position information of the boundary particles obtained in step 1.2, velocities are sequentially given to all the boundary particles along the normal vector direction by the numerical calculation method of formula 1.

After the relaxation factor is introduced and a speed correction formula is adopted, reasonable boundary particle speed can be obtained quickly, so that the boundary particle speed meets Dirichlet boundary conditions, the microscopic state of a simulation system is met, and system oscillation is reduced.

3. Applying appropriate forces to the boundary particles

In order to ensure that the boundary particles are consistent with the internal stress and the microscopic action potential of the simulation system in the microscopic state, the boundary action force is introduced to ensure the continuity of the internal stress of the system and adapt to the local particle distribution function. According to the local particle distribution function, virtual 'external particles' are introduced, the action of the external boundary on the internal boundary particles is reconstructed, and the additional acting force applied to the boundary particles is deduced based on the Viry theorem. The specific implementation mode is as follows:

in practical simulations, to apply boundary forces within a boundary domain, one can assume that edges are at handOutside the boundary region, there are some outer particles which exert an action potential on the boundary particles, the outer particles having a radius of truncation R_CEffective potentials exist within the range for the boundary particles, and the sum of these potentials constitutes the effective boundary force. The action potential among all particles in the molecular dynamics simulation model is always limited to the truncation radius R_CIn scope, therefore, it can be assumed that any boundary particle is surrounded by a shell outside the boundary region, the radius R of the shell being x (the distance from the particle in the boundary region to the boundary in the normal direction) to R_CWithin a range.

The number of particles in the shell can be expressed as: rho g (r) S (r) dr, the acting force of the particles in the shell on the boundary particles exists, and along the normal direction of the simulation domain, the acting force of the particles in any thin shell on any particle in the boundary domain is as follows:wherein rho is the particle density in the molecular dynamics simulation system;

g (r) is a function of the radial distribution of the particles; s (r) is the area of the shell; t (r) is the projection of the shell area on the outer layer of the boundary domain; v (r) is the action potential of a single particle in the shell on any particle in the boundary domain; dr is the thickness of the shell.

After the acting force of the single-layer thin shell on any boundary particle in the boundary domain is determined, the total acting force applied by the external particle on any boundary particle in the boundary domain can be obtained as follows:

in the above formula, Φ is the stress tensor; p is Viry pressure; t (r, x) is the projection of the thin shell area on the outer layer of the boundary domain when the boundary particle has the coordinate x.

As shown in FIG. 4, for the boundary particle i, at its cutoff radius R_cWithin the range, the acting force of any external particle j outside the boundary area on the external particle j is F_ijThe integral calculation of formula 2 is used to obtain the boundary particles to which all particles are applied within the cutoff radius range outside the boundary region by the thin shell methodEffective force F on i_mFinally, the proper boundary acting force is applied to all the boundary particles, so that the macroscopic physical quantity of the boundary particles is consistent with the microscopic state, and the complete molecular dynamics boundary condition is obtained.

The above embodiments are only used to illustrate the technical solutions of the present application, and not to limit the same; although the present application has been described in detail with reference to the foregoing embodiments, it should be understood by those of ordinary skill in the art that: the technical solutions described in the foregoing embodiments may still be modified, or some technical features may be equivalently replaced; such modifications and substitutions do not substantially depart from the spirit and scope of the embodiments of the present application and are intended to be included within the scope of the present application.