物质和材料的跨时间尺度计算机仿真模拟方法及装置
阅读说明：本技术 物质和材料的跨时间尺度计算机仿真模拟方法及装置 (Cross-time-scale computer simulation method and device for substances and materials ) 是由 万亮 梅青松 刘浩文 于 2021-07-09 设计创作，主要内容包括：本发明提供一种物质和材料的跨时间尺度计算机仿真模拟方法及装置,包括针对待进行计算机仿真模拟的物质和材料,构建多粒子模型；采用多粒子系统的四维扩展Langevin动力学方程组计算物质和材料系统的微观状态对应的所有粒子的坐标X与速度V的值随时间变化的情况,其中粒子的近邻偏心度量的值通过计算以该粒子为起点到其各最近邻粒子的三维欧氏空间向量之和再平均所得到的近邻偏心位移的3个坐标轴分量的绝对值之和来得到；根据计算所得每个离散时间步所有粒子的坐标X与速度V的值,对相应的物质和材料系统进行热力学与动力学分析；其中没有进行多粒子系统初始参考状态设定的步骤。本发明技术方案简化了实现流程,使用时更便捷,运行效率明显提升。(The invention provides a cross-time scale computer simulation method and device of substances and materials, which comprises the steps of constructing a multi-particle model aiming at the substances and materials to be subjected to computer simulation; calculating the time-varying conditions of the coordinate X and the velocity V of all particles corresponding to the microscopic state of a substance and material system by adopting a four-dimensional extended Langevin kinetic equation set of a multi-particle system, wherein the value of the neighbor eccentric measurement of the particle is obtained by calculating the sum of the absolute values of 3 coordinate axis components of neighbor eccentric displacement obtained by averaging the sum of three-dimensional Euclidean space vectors from the particle serving as a starting point to each nearest neighbor particle; performing thermodynamic and kinetic analysis on corresponding substances and material systems according to the calculated values of the coordinates X and the speed V of all particles in each discrete time step; wherein no step of initial reference state setting of the multi-particle system is performed. The technical scheme of the invention simplifies the implementation process, is more convenient to use and obviously improves the operation efficiency.)
技术领域
本发明属于凝聚态物质和材料学与计算机仿真工程研究领域，特别涉及一种基于多粒子凝聚态系统加速分子动力学模型的计算机数值仿真模拟的改进技术方案。
背景技术
分子动力学计算机数值仿真模拟方法是目前物质研究领域内的一种最为主要的计算模拟方法，因其能直接仿真模拟给出物质在原子和分子尺度层面上的实时动态过程，被广泛地用在材料、物理、化学、生命科学等物质研究领域。在目前普遍使用的传统分子动力学模拟方法中，其基本原理可描述如下。
对于含有N个粒子(原子或分子)的任意物质和材料的多粒子系统，分子动力学模拟方法一般用所有粒子在三维欧氏空间中的坐标X＝(x1,x2,…,xN)与相应的粒子动量P＝(p1,p2,…,pN)(pi＝mivi，mi为粒子i的质量，vi为粒子i的速度，i∈{1,2,…,N})作为动力学量并构成相空间来描述多粒子系统的任意微观状态，再通过一组特定的动力学方程组来描述多粒子系统的微观状态随时间演化的动力学过程。在传统的分子动力学模拟方法中，常用的动力学方程组包括(参见文献[1])：
A.牛顿动力学(1687年)
Fi＝miai
B.Langevin动力学(1908年)
C.Nose-Hoover扩展动力学(1984年)
上述公式中，各个符号所对应的物理量分别为：
F—粒子所受作用力(一般为保守力)；
m—粒子的质量；
v—粒子的速度(x为粒子在三维欧氏空间中的坐标)；
a—粒子的加速度(x为粒子在三维欧氏空间中的坐标)；
γ—粒子在Langevin动力学假想溶液里运动的粘滞系数；
T—多粒子系统的热浴温度；
R(t)—白噪声函数；
kB—Boltzmann常数；
ζ—多粒子系统的Nose动力学扩展变量；
Q—Nose动力学扩展变量ζ对应的质量；
N—系统中所有粒子(原子或分子)的数目。
所有物理量符号中的下标i代表该物理量属于多粒子系统中标号为i的任意某粒子，物理量符号上的单点与双点分别表示对时间求一阶导数和二阶导数。
在上述模型框架下，为了较为真实可信地反映物质中原子、分子的运动情况，其动力学方程组数值离散化处理后的时间步长Δt一般不超过数个飞秒(约为通常固体晶格原子振动周期的5％)，从而使得其在目前计算机的运算能力下，最多只能模拟纳秒量级时间尺度的物理过程。而通常人们较为关注的物理现象如薄膜生长、固体相变、材料形变、材料的热处理等过程中的微观结构演化过程都涉及毫秒以上甚至数天、数年(如材料蠕变)的时间跨度。例如，作为一种铁碳合金，钢铁材料在500℃下其碳原子在铁原子晶格间跳跃1次需要等待的时间约为1微秒。如此巨大(相差约3～13个数量级)的时间尺度限制使得人们难以采用传统的分子动力学模拟方法来研究通常条件下物质和材料中所发生的物理过程及其规律。
尽管目前已有一些开源软件(参见文献[2-4])和商业软件(参见文献[5-8])将分子动力学仿真模拟方法包含到了软件的功能之中，但该方法基本上还只是停留在学术界实验室范围内被相关研究人员有限度地使用，而上述时间尺度限制正是阻碍其走向工业与工程实践从而被更为广泛地应用的关键瓶颈所在。如果该瓶颈问题可以被解决与突破，可以预见将会为相关工业领域特别是材料工业的研发方式带来崭新的面貌。届时分子动力学计算机数值仿真模拟方法将会成为类似于电子显微镜的材料研究必备工具，同时由于高性能计算机相比于昂贵的实验设备在成本方面的巨大优势，新材料研发的时间周期与经济成本都将会大幅缩小，从而有效地促进各种技术的发展。
针对这个时间尺度上的瓶颈问题，世界范围内的大量研究团队为此进行了持续不懈的研究，到目前为止发展并提出了一些所谓加速分子动力学模拟方法的解决方案。这些解决方案从技术思路上大体可分为三类：
1○系统势能面(Potential Energy Surface)抬升模拟方法，包括Hyperdynamics(参见文献[9,10])、Metadynamics(参见文献[11,12])、Strain-boost moleculardynamics(参见文献[13])、Adaptive boost molecular dynamics(参见文献[14])等方法与模型，其特点是通过为多粒子(原子、分子)凝聚态系统构造一个偏置势(Biaspotential)ΔU(X)，将这个偏置势叠加到系统本来的原子间相互作用势函数U(X)上，系统在这个叠加势的基础上基于传统分子动力学方程进行时间演化，将会有更大的几率从一个初态越过中间态势垒而过渡到另一个状态，从而达到了时间尺度跨越的效果。这种模型与模拟方法的缺点在于，由于多粒子系统的势能面是一个3N维空间中的超曲面，随着粒子数N的增加这个超曲面的形貌会变得异常复杂，因此需要小心注意偏置势的叠加是否会违反化学反应动力学基本规律；同时，偏置势的构造没有固定规则可循，且一般较为复杂，仅能针对特定的动力学过程如原子扩散、位错形核与运动等进行计算模拟，从而限制了其广泛应用。
2○多复本计算模拟方法，主要包括Parallel Replica Dynamics(参见文献[10,15])、Parallel Trajectory Splicing(参见文献[16])等方法与模型，其特点是基于多粒子凝聚态系统动力学过程的随机与统计特性，通过构造大量的系统初态的复本，为每个复本赋予不同的随机指标值，各个复本再基于传统分子动力学方程进行时间演化，由于每个复本发生状态改变的几率相同，N个复本同时进行时间演化等价于将原有单个系统中发生微观结构转变事件的几率提高了N倍。可以看出，由于要为每个复本进行同样的分子动力学计算模拟，这个方法实质上是以CPU计算资源换取模拟时间尺度的放大，其缺点在于对CPU数目要求较高，因此只适合于粒子数较少的物质系统的模型，从而限制了其广泛应用。
3○高温-低温映射模拟方法，主要包括Temperature Accelerated Dynamics方法和模型(参见文献[10,17,18])，其特点是根据化学反应动力学规律，多粒子凝聚态系统在温度升高时，系统发生微观状态变化的几率会大幅增加，因此该方法是通过在一个较高温度下对系统进行分子动力学时间演化，将所有高温下发生的微观结构转变事件记录下，再计算这些事件在目标低温下发生的几率，选择几率较大的事件作为低温下动力学状态变化的过程，从而达到时间尺度跨越的效果。该方法的主要缺点在于，由于同一个系统在较高温度和较低温度下的微观状态变化特征会存在较为显著的差别，较高温度下的分子动力学模拟不一定能有效合理地给出较低温度下的微观状态变化情况，因此可能会产生动力学物理失真。
除了上述三种技术方案外，通过借鉴由Shuichi Nose发展Nose-Hoover动力学时所提出的动力学扩展变量的思想(参见文献[19,20])，同时基于为描述多粒子凝聚态系统中特定微观结构转变事件动力学过程而引入的多粒子动力学系统“集合变量(Collectivevariable)”的概念(参见文献[11,21])，Luca Maragliano和Eric Vanden-Eijnden于2006年发表文章(参见文献[22])提出，可以将多粒子动力学系统的集合变量对应的热力学与统计力学量设定为动力学扩展变量，再采用一组扩展Langevin动力学方程组来描述系统的微观状态演化的动力学过程，实现对与集合变量的定义表达式相对应的特定微观结构转变事件进行分子动力学模拟的加速效果。该方法的有效性和可靠性强烈依赖于集合变量的构造方式的合理性(参见文献[21,23])。尽管有不少学者和研究人员尝试在该方法框架内改进集合变量的构造方式(如人工智能方法的使用，参见文献[24,25])以便应用于加速分子动力动力学仿真模拟，然而到目前为止，还没有一个较为合理可靠、普遍通用、简便有效的技术方案，来支持实现对多粒子凝聚态系统的分子动力学计算机数值仿真模拟获得较好的时间尺度跨越效果。
申请人的研究团队于2020年8月提出了一个专利申请《一种用于物质和材料的计算机仿真模拟方法及装置》已获得授权，专利号为202010884514.3，包括针对待进行计算机仿真模拟的物质和材料，构建一个该物质和材料系统的多粒子模型，用所有粒子的坐标与动量描述系统的微观状态，所述粒子为原子或分子；设置所述多粒子模型中每个粒子的粒子曳步运动度量D，用于描述模型中可以发生的任意微观结构转变过程的反应坐标，设定多粒子模型的初始参考状态；采用多粒子系统的四维扩展Langevin动力学方程组计算物质和材料系统的微观状态随时间演化的情况，该方程组为每个粒子赋予一个额外的动力学量s，s与D以谐振子的方式耦合；根据获取的每个离散时间步下模型中所有粒子的坐标与速度的值，对相应的物质和材料系统进行热力学与动力学分析。
但是，对于该方法所设置的用于描述模型中微观结构转变过程反应坐标的粒子曳步运动度量D，存在如下缺点：(1)对粒子曳步运动度量D进行计算时，需要事先设定多粒子模型的初始参考状态，同时需要对计算粒子曳步运动度量D所依据的参考状态每隔一定时间步数n用当前时刻系统的微观状态进行更新，n参量的具体数值在使用过程中需根据待进行计算机仿真模拟的物质和材料的实际情况来设定，这给该方法的实施和使用增加了一些复杂度；(2)粒子曳步运动度量D的计算公式较为繁琐，使得对于将额外的动力学量s与D以谐振子的方式耦合所得到的多粒子系统的四维扩展Langevin动力学方程组，在对其进行计算求解时，所需的计算机程序编制工作较繁琐，且计算量比较大，降低了该方法使用过程中的执行效率。
本发明相关文献：
[1]陈敏伯.计算化学：从理论化学到分子模拟.北京:科学出版社,2009.
[2]LAMMPS Molecular Dynamics Simulator.https://lammps.sandia.gov/.
[3]NAMD-Scalable Molecular Dynamics.http://www.ks.uiuc.edu/Research/namd/.
[4]GROMACS.http://www.gromacs.org/.
[5]BIOVIA Materials Studio.https://3dsbiovia.com/products/collaborative-science/biovia-materials-studio/.
[6]QuantumATK.https://www.synopsys.com/silicon/quantumatk.html.
[7]DL_POLY Molecular Simulation Package.https://www.scd.stfc.ac.uk/Pages/DL_POLY.aspx.
[8]The Vienna Ab initio Simulation Package.https://www.vasp.at/.
[9]VOTER A F.Hyperdynamics:Accelerated Molecular Dynamics ofInfrequent Events.Phys Rev Lett,1997,78(20):3908-11.
[10]VOTER AF,MONTALENTI F,GERMANN T C.Extending the timescale inatomistic simulation of materials.Annual Review of Materials Research,2002,32(1):321-46.
[11]LAIO A,PARRINELLO M.Escaping free-energy minima.Proceedings ofthe National Academy of Sciences of the United States of America,2002,99(20):12562-6.
[12]LAIO A,GERVASIO F L.Metadynamics:a method to simulate rare eventsand reconstruct the free energy in biophysics,chemistry and materialscience.Rep Prog Phys,2008,71(12):126601.
[13]HARA S,LI J.Adaptive strain-boost hyperdynamics simulations ofstress-driven atomic processes.Physical Review B,2010,82(18):184114.
[14]ISHII A,OGATA S,KIMIZUKA H,et al.Adaptive-boost moleculardynamics simulation of carbon diffusion in iron.Physical Review B,2012,85(6):064303.
[15]VOTER AF.Parallel replica method for dynamics of infrequentevents.Physical review B,Condensed matter,1998,57(22):R13985-8.
[16]PEREZ D,CUBUK E D,WATERLAND A,et al.Long-time dynamics throughparallel trajectory splicing.Journal of Chemical Theory&Computation,2015,acs.jctc.5b00916.
[17]SORENSEN M R,VOTER AF J.Temperature-accelerated dynamics forsimulation of infrequent events.Journal of Chemical Physics,2000,112(21):9599-606.
[18]ZAMORA R J,UBERUAGA B P,PEREZ D,et al.The Modern Temperature-Accelerated Dynamics Approach.Annu Rev Chem Biomol Eng,2016,7(1):87-110.
[19]NOSE S.A Unified Formulation of the Constant TemperatureMolecular-Dynamics Methods.Journal of Chemical Physics,1984,81(1):511-9.
[20]NOSE S.A molecular dynamics method for simulations in thecanonical ensemble.Molecular Physics,1984,52(2):255-68.
[21]FIORIN G,KLEIN M L.Using collective variables to drive moleculardynamics simulations.Molecular Physics,2013,111(22-23):3345-62.
[22]MARAGLIANO L,VANDEN-EIJNDEN E.A temperature accelerated methodfor sampling free energy and determining reaction pathways in rare eventssimulations.Chemical Physics Letters,2006,426(1-3):168-75.
[23]VALSSON O,TIWARY P,PARRINELLO M.Enhancing Important Fluctuations:Rare Events and Metadynamics from a Conceptual Viewpoint.Annual review ofphysical chemistry,2016,67(1):annurev-physchem-040215-112229.
[24]ZHANG J,CHEN M.Unfolding Hidden Barriers by Active EnhancedSampling.Physical Review Letters,2017,121(1):010601.
[25]CHEN W,TAN A R,FERGUSON A L.Collective variable discovery andenhanced sampling using autoencoders:Innovations in network architecture anderror function design.Journal of Chemical Physics,2018,149(7):072312.
[26]弗兰克(D.FRENKEL)等.分子模拟——从算法到应用.北京:化学工业出版社,2004.
[27]Interatomic Potentials Repository.https://www.ctcms.nist.gov/potentials/.
[28]BECQUART C S,RAULOT J M,BENCTEUX G,et al.Atomistic modeling of anFe system with a small concentration of C.Computational Materials Science,2007,40(1):119-29.
[29]VEIGA R G,PEREZ M,BECQUART C S,et al.Atomistic modeling of carbonCottrell atmospheres in bcc iron.Journal of physics Condensed matter,2013,25(2):025401.
[30]APOSTOL F,MISHIN Y.Angular-dependent interatomic potential forthe aluminum-hydrogen system.Physical Review B,2010,82(14):144115.
发明内容
本发明的目的在于，针对现有技术瓶颈，提出一种新的技术方案来进行在计算机上仿真模拟多粒子(原子或分子)物质和材料系统的微观状态随时间演化的情况，以实现对物质和材料系统进行高可靠、高效率、多用途的计算机仿真分析。
本发明的技术方案提供一种物质和材料的跨时间尺度计算机仿真模拟方法，包括以下步骤，
步骤1，针对待进行计算机仿真模拟的物质和材料，构建一个该物质和材料系统的多粒子模型；
步骤2，采用多粒子系统的四维扩展Langevin动力学方程组来计算物质和材料系统的微观状态对应的所有粒子的坐标X与速度V的值随时间变化的情况；
步骤3，根据物质和材料系统的多粒子模型在每个离散时间步对应的所有粒子的坐标X与速度V的值，对相应的物质和材料系统进行热力学与动力学分析，
在步骤1和步骤2之间，没有进行多粒子系统初始参考状态设定的步骤，
在步骤2中，进行以下操作，
所述多粒子系统的四维扩展Langevin方程组给多粒子系统中每个粒子赋予了一个额外的动力学量s，所有粒子的额外动力学量构成N维向量S＝(s1,s2,s3,…,sN)，所述多粒子系统中所有粒子的总数为N；
设定每个粒子的额外动力学量s为与该粒子的近邻偏心度量D相对应的热力学与统计力学量；
所述粒子的近邻偏心度量D＝D(X)的值通过计算以该粒子为起点到其各最近邻粒子的三维欧氏空间向量之和再平均所得到的近邻偏心位移R(X)的3个坐标轴分量的绝对值之和来得到，所有粒子的近邻偏心度量构成N维向量D＝(D1(X),D2(X),D3(X),…,DN(X))；所述粒子的近邻偏心度量D用于描述该粒子与其最近邻粒子之间的非协同运动的程度；
由于粒子的近邻偏心度量D与粒子在三维欧氏空间中坐标的坐标轴分量具有相同的物理量纲，多粒子系统被认为是在一个四维空间中运动，在这个四维空间中，将每个粒子在原三维欧氏空间中所定义的近邻偏心度量D＝D(X)与该粒子额外的动力学量s以谐振子的方式进行耦合，得到所有粒子在四维空间中粒子间相互作用的势能函数；在该势能函数的作用下，粒子在原三维欧氏空间中的动力学量X与粒子额外的动力学量S分别以相互独立的粘滞系数和热浴温度各自依据普通Langevin动力学方程进行状态时间演化；
对该四维扩展Langevin动力学方程组进行时间离散化数值处理，按时间步顺序自动化逐步运算求解多粒子模型的所有粒子的坐标X与速度V在每个离散时间步的数值，得到物质和材料系统的微观状态对应的所有粒子的坐标X与速度V的值随时间变化的情况。
而且，步骤2中，粒子的近邻偏心度量D的定义如下，
其中，
N—多粒子系统中的总粒子数；
X—所有粒子在三维欧氏空间中的坐标，X＝(x1,x2,…,xN)，
Di(X)—i粒子的近邻偏心度量，是所有粒子的坐标X的函数；
Ri(X)—i粒子的近邻偏心位移，是所有粒子的坐标X的函数；
—i粒子的近邻偏心位移Ri(X)在三维欧氏空间中三个坐标轴方向的分量，是所有粒子的坐标X的函数；
Ni—i粒子的以Rd为半径范围内所有近邻粒子的数目；
—i粒子的以Rd为半径范围内所有近邻粒子的集合；
xi—当前时刻下，i粒子在三维欧氏空间中的坐标；
xj—当前时刻下，j粒子在三维欧氏空间中的坐标。
而且，步骤2中，对于多粒子系统的所有粒子的坐标X与速度V的值如何在粒子间相互作用势能函数U(X)的作用下随时间变化，采用如下形式的四维扩展Langevin动力学方程组来进行计算：
其中，
N—多粒子系统中的总粒子数；
X—所有粒子在三维欧氏空间中的坐标，X＝(x1,x2,…,xN)，
S—多粒子系统额外的动力学量，S＝(s1,s2,s3,…,sN)；
—分别为i粒子在三维欧氏空间中α坐标轴方向的坐标、速度、加速度，
si，—分别为i粒子额外的动力学量、其对时间的一阶导数、其对时间的二阶导数；
mi—i粒子的质量；
γ—粒子在Langevin动力学假想溶液里运动的粘滞系数；
T—多粒子系统的主热浴温度；
kB—Boltzmann常数；
—与额外的动力学量S对应的、粒子在Langevin动力学假想溶液中运动的粘滞系数；
—与额外的动力学量S对应的、多粒子系统的副热浴温度；
Rx(t),Rs(t)—相互独立的白噪声函数；
U(X)—多粒子系统在三维欧氏空间中的粒子间相互作用的势能函数，是所有粒子的坐标X的函数；
sj—j粒子额外的动力学量；
Dj(X)—j粒子的近邻偏心度量，是所有粒子的坐标X的函数；
κ—多粒子系统额外的动力学量S与粒子的近邻偏心度量D之间的力耦合参量；
Uκ(X,S)—多粒子系统在四维空间中粒子间相互作用对应的势能函数，是所有粒子的坐标X与额外的动力学量S的函数。
而且，步骤2中，对所述四维扩展Langevin动力学方程组进行时间离散化数值处理，设置离散时间步长为Δt，离散时间步长Δt采用固定时间步长或非固定时间步长，自动化逐步运算将物质和材料的多粒子模型每个离散时间步对应的所有粒子的坐标X与速度V的值求解出，给出多粒子凝聚态物质和材料系统的微观状态随时间演化的情况。
而且，用于仿真分析铁碳合金中的间隙碳原子和铁原子空位在体心立方铁晶格中的扩散迁移动力学过程。
或者，用于仿真分析在溶解有一定浓度氢原子的金属铝双晶体中，氢原子在铝双晶体中扩散以及在晶界上偏聚的动力学过程。
本发明提供一种物质和材料的跨时间尺度计算机仿真模拟装置，用于实现如上所述物质和材料的跨时间尺度计算机仿真模拟方法。
本发明的优势在于：由于在本发明中所提出和定义的粒子的近邻偏心度量D＝(D1(X),D2(X),D3(X),…,DN(X))本质上反映了多粒子系统中每个粒子与其最近邻粒子之间的非协同运动(Shuffling，也称“曳步运动”)的程度，这种非协同运动的程度通常与多粒子凝聚态物质与材料系统中任意内部微观结构转变过程的反应坐标相对应，据此本发明巧妙地利用粒子的近邻偏心度量D来构造多粒子动力学系统的集合变量，再代入到LucaMaragliano和Eric Vanden-Eijnden所提出的(参见文献[22])扩展Langevin动力学方程组中，得到多粒子系统的四维扩展Langevin动力学方程组来描述其微观状态对应的所有粒子的坐标X与速度V的值随时间演化的情况，按照Luca Maragliano和Eric Vanden-Eijnden所给出的分析(参见文献[22])，只要合理地设置参数κ、γ与的值，再通过调节与额外的动力学量S对应的、多粒子系统的副热浴温度的值，就可以控制与增强多粒子凝聚态物质与材料系统中每个粒子的非协同运动程度的涨落幅度，实现对多粒子凝聚态物质和材料系统中任意内部微观结构转变过程的加速分子动力学计算机仿真模拟，在采用与传统的普通分子动力学模拟方法相当的数值模拟时间步长Δt的条件下，在有限的、与传统的普通分子动力学模拟方法相当的计算模拟时间步数和计算量内，对凝聚态物质和材料中真实物理时间跨度从1皮秒到数天范围不等的微观结构演化物理过程进行原子尺度计算机数值仿真模拟研究，从而可以将分子动力学模拟的时间尺度范围在传统的普通分子动力学模拟方法的基础上大幅拓展(增加约1～13个数量级)，在此基础上能有效地分析凝聚态物质和材料系统在一定宏观热力学条件下所发生的内部微观结构转变的过程及其动力学特征与规律，并能针对所研究的凝聚态物质和材料系统进行各种物理和化学特性的分析，具有重要的应用价值和广阔的商业前景。基于上述特征，在此将本发明所提出的通用加速分子动力学计算机仿真模拟的改进方法命名为“曳步运动加速分子动力学(Shuffling Accelerated MolecularDynamics(SAMD))”。
与申请人的研究团队于2020年8月所提出的专利号为202010884514.3的专利方法《一种用于物质和材料的计算机仿真模拟方法及装置》相比较，本发明在步骤2中所设置的用于描述任意粒子与其最近邻粒子之间非协同运动程度的粒子近邻偏心度量D，其作用与专利号为202010884514.3的专利方法中所设置的粒子曳步运动度量D具有相同的功能，即都可以用于描述模型中微观结构转变过程的反应坐标。但在本发明方法中所设置的粒子近邻偏心度量D的计算公式中，不需要进行多粒子系统参考状态的设定，因此本发明方法在实施使用时相对于专利号为202010884514.3的专利所给出的方法要更为便捷；而且本发明方法中所设置的粒子近邻偏心度量D的计算公式表达式相比于专利号为202010884514.3的专利方法中所设置的粒子曳步运动度量D的计算表达式要大幅简化，因此在实际的计算机仿真模拟执行中，本发明方法更易于进行计算机程序编制，在计算效率上也更为优越。
与申请人的研究团队于2021年5月所提交的专利申请《一种用于物质和材料的新型计算机仿真模拟方法及装置》相比较，本发明在步骤2中所设置的用于描述任意粒子与其最近邻粒子之间非协同运动程度的粒子近邻偏心度量D，其作用与申请名为《一种用于物质和材料的新型计算机仿真模拟方法及装置》的专利方法中所设置的粒子近邻偏心度量D具有相同的功能，但在计算表达式上，申请名为《一种用于物质和材料的新型计算机仿真模拟方法及装置》的方法中，使用以某个粒子为起点到其各最近邻粒子的三维欧氏空间向量之和再平均所得到的近邻偏心位移R(X)的向量积平方根来计算该粒子的粒子近邻偏心度量D，这种计算方式所得到的粒子近邻偏心度量D是粒子位移的非线性函数，而本发明则使用近邻偏心位移R(X)的3个坐标轴分量的绝对值之和来计算各粒子的粒子近邻偏心度量D，这种计算方式所得到的粒子近邻偏心度量D是粒子位移的线性函数，可以更为准确的描述和反映多粒子模型中微观结构转变过程的反应坐标，因而在对物质和材料实际进行计算机仿真模拟时，可以得到保真度更高的计算机仿真模拟数据结果。需要说明的是，同样作为对专利号为202010884514.3的专利方法《一种用于物质和材料的计算机仿真模拟方法及装置》的改进，相比于2021年5月所提交的申请名为《一种用于物质和材料的新型计算机仿真模拟方法及装置》的方法，本发明所提出的新的计算方法是申请人研究团队从另一个角度和思路出发，经过大量探索研究实践后得到的新的改进途径。
本发明方案实施简单方便，实用性强，可靠性高，解决了相关技术存在的实用性低、保真性不高以及实际应用不便的问题，能够提高用户体验，具有重要的市场价值。
附图说明
图1是本发明实施例中所提出的SAMD仿真模拟方法的主要执行步骤示意图。
图2是本发明具体应用例子1中所讨论的含有1个铁原子空位与3个间隙碳原子的铁晶体的原子模型示意图，子图(A)显示的是具有体心立方结构的铁晶体的晶格点阵，子图(B)显示的是铁晶体中铁原子空位与间隙碳原子的示意图。
图3显示的是基于图2(B)的模型，在给定系统的实际热力学温度为300K下，由普通分子动力学(Conventional MD)方法和本发明所提出的SAMD仿真模拟方法分别进行计算机数值模拟所给出的碳原子均方根位移随模拟时间步数变化的情况。
图4显示的是基于图2(B)的模型，采用本发明所提出的SAMD仿真模拟方法，在设定系统的副热浴温度为时进行计算机数值模拟，所得到的模型在图3中箭头所示不同时间步(nsteps)对应的时刻下间隙碳原子与铁原子空位所处位置与迁移情况的结果：(A)nsteps＝0；(B)nsteps＝5e+05；(C)nsteps＝1e+06；(D)nsteps＝1.5e+06；(E)nsteps＝2e+06。
图5是本发明具体应用例子2中所讨论的溶解有一定浓度氢原子的金属铝双晶体的原子模型示意图，子图(A)显示的是溶解有一定浓度氢原子的铝双晶体模型的投影图(大球为铝原子，小球为氢原子，氢原子分布未达到平衡态)，子图(B)显示的是将所有铝原子隐藏后仅显示氢原子初始分布状况的情况。
图6显示的是基于图5的模型，在给定系统的实际热力学温度为100K下，由普通分子动力学(Conventional MD)方法和本发明所提出的SAMD仿真模拟方法分别进行计算机数值模拟所给出的氢原子均方根位移随模拟时间步数变化的情况。
图7显示的是基于图5的模型，采用本发明所提出的SAMD仿真模拟方法，在设定系统的副热浴温度为时进行计算机数值模拟，所得到的模型在图6中箭头所示不同时间步(nsteps)对应的时刻下氢原子的分布与迁移情况的结果：(A)nsteps＝0；(B)nsteps＝5e+05；(C)nsteps＝1e+06；(D)nsteps＝1.5e+06；(E)nsteps＝2e+06。
图8显示的是基于图5的模型，采用本发明所提出的SAMD仿真模拟方法，在设定系统的副热浴温度分别为500K(子图(A))、1000K(子图(B))、1500K(子图(C))时进行计算机数值模拟，在总时间步数为2.0×106的模拟时间长度下，模拟最终得到的氢原子分布的结果；子图(D)为采用蒙特卡罗(MC)计算机数值模拟方法，同样在模型系统的实际热力学温度为100K下进行模拟所得到的氢原子分布的结果，反映出氢原子在该热力学温度下的平衡分布状态。
图9是现有技术中的一种用于物质和材料的计算机仿真模拟方法流程图。
具体实施方式
以下结合附图和实施例详细说明本发明技术方案。
本发明采用一种可称之为粒子的“近邻偏心度量(nearest neighbordecentration)”的物理量作为多粒子系统动力学的“集合变量(collective variable)”，并定义此集合变量对应的热力学和统计力学量为多粒子系统的动力学扩展变量，在此基础上提出了一种可以命名为“曳步运动加速分子动力学(Shuffling Accelerated MolecularDynamics(SAMD))”的通用加速分子动力学计算模拟的改进方法(以下简称SAMD仿真模拟方法)，通过在计算机上数值计算求解下文将要给出的多粒子系统的四维扩展Langevin动力学方程组，来仿真模拟多粒子凝聚态物质和材料系统的微观状态随时间演化的情况，并用于分析凝聚态物质和材料系统在一定宏观热力学条件下所发生的内部微观结构转变的过程及其动力学特征与规律。
需要注意的是，本文中所有涉及到“粒子”一词的使用，其英文对应为“particle”，按英文维基百科全书，指的是在物理和化学上尺寸较小、具有体积和质量的物体的局部，根据其尺寸范围可以涵盖从亚原子粒子如电子、质子、中子等，到微观粒子如原子、分子，再到介观颗粒如生物大分子、胶体粒子，甚至可以为宏观颗粒如粉粒和其它粒状物质等。本文给出的关于本发明技术与方法的阐述中所涉及的“粒子”则主要针对的是构成物质和材料的微观粒子即原子或分子，其主要特征之一是在物理理论中可以采用质点(mass point)的概念来描述。这并不影响本文所描述的发明技术和方法用在各种其它粒子如生物大分子、胶体粒子和其它粒状物质等作为研究对象的场合，因而其相应场合下的使用同样受本发明专利的保护。
参见图1，本发明实施例所提出的物质和材料的跨时间尺度计算机仿真模拟方法(SAMD仿真模拟方法)主要包括以下步骤：
步骤1：针对所研究和需要进行计算机仿真模拟的物质和材料，构建一个该物质和材料系统的多粒子模型，实现如下：
对于含有N个粒子(原子或分子)的任意多粒子凝聚态物质和材料系统，所有粒子构成的集合记为{1,2,…,N}，用所有粒子在三维欧氏空间(用笛卡尔坐标系表示，为三维欧氏空间的三个坐标轴)中的坐标X＝(x1,x2,…,xN)与相应的粒子动量P＝(p1,p2,…,pN)作为动力学量并构成相空间来描述多粒子系统的任意微观状态，任意某个粒子i的坐标为粒子i的动量为pi＝mivi，mi为粒子i的质量，vi为粒子i的速度，i∈{1,2,…,N}，所有粒子的速度为V＝(v1,v2,…,vN)，xi与vi和pi均为3维向量，X与V和P均为3N维向量。这些粒子间存在相互作用力，用U(X)表示系统所有粒子间的相互作用势能函数，则任意某个粒子i在α坐标轴方向所受的力可以表示为其中该粒子在三维欧氏空间中α坐标轴方向的坐标为
该物质和材料系统的多粒子模型的构建需要给出N,M＝(m1,m2,…,mN)的数值、U(X)函数的数据和X,V的初始数值。具体实施时，通过给定N,M的值(N,M一般为常数)、U(X)函数的数据和X,V的初值就构建了一个凝聚态物质和材料系统的多粒子模型，接下来的任务是确定该多粒子系统的X,V的值如何在粒子间相互作用势能函数U(X)的作用下随时间变化。
步骤2：通常研究多粒子凝聚态物质和材料系统，更主要关注的是其中所发生的各种内部微观结构转变的动力学过程，而不是每个粒子在其平衡位置附近的随机热振动。本发明提出，采用如下公式所定义的粒子的近邻偏心度量(nearest neighbordecentration)的物理量D＝D(X)来描述多粒子系统中每个粒子与其最近邻粒子之间的非协同运动(Shuffling，也称“曳步运动”)的程度：
其中，
N—多粒子系统中的总粒子数；
X—所有粒子在三维欧氏空间中的坐标，X＝(x1,x2,…,xN)，
Di(X)—i粒子的近邻偏心度量，是所有粒子的坐标X的函数；
Ri(X)—i粒子的近邻偏心位移，是所有粒子的坐标X的函数；
—i粒子的近邻偏心位移Ri(X)在三维欧氏空间中三个坐标轴方向的分量，是所有粒子的坐标X的函数；
Ni—i粒子的以Rd为半径范围内所有近邻粒子的数目；
—i粒子的以Rd为半径范围内所有近邻粒子的集合；
xi—当前时刻下，i粒子在三维欧氏空间中的坐标；
xj—当前时刻下，j粒子在三维欧氏空间中的坐标。
所有粒子的近邻偏心度量构成N维向量D＝(D1(X),D2(X),D3(X),…,DN(X))；由该公式所给出的粒子的非协同运动的程度通常与多粒子凝聚态物质与材料系统中任意内部微观结构转变过程的反应坐标相对应，因而由该公式所定义的粒子的近邻偏心度量D可以被用来描述多粒子凝聚态物质和材料系统中任意内部微观结构转变过程的反应坐标。
从上述粒子的近邻偏心度量D的定义表达式可以看出，只要给定了参量Rd的值，多粒子系统在任意时刻下每个粒子的近邻偏心度量D＝D(X)的值都可以基于当前时刻所有粒子的坐标X的值根据上述表达式直接求出。
在上述粒子的近邻偏心度量D的定义的基础上，对于多粒子系统的所有粒子的坐标X与速度V的值如何在粒子间相互作用势能函数U(X)的作用下随时间变化，本发明提出可以采用如下形式的四维扩展Langevin动力学方程组来进行计算：
该方程组表达式中各符号所代表的物理量或者其数学/物理意义分别为：
N—多粒子系统中的总粒子数；
X—所有粒子在三维欧氏空间中的坐标，X＝(x1,x2,…,xN)，
S—多粒子系统额外的动力学量，S＝(s1,s2,s3,…,sN)；
—分别为i粒子在三维欧氏空间中α坐标轴方向的坐标、速度(即坐标对时间的一阶导数)、加速度(即坐标对时间的二阶导数)，
si，—分别为i粒子额外的动力学量、其对时间的一阶导数、其对时间的二阶导数；
mi—i粒子的质量；
γ—粒子在Langevin动力学假想溶液里运动的粘滞系数；
T—多粒子系统的主热浴温度；
kB—Boltzmann常数；
—与额外的动力学量S对应的、粒子在Langevin动力学假想溶液中运动的粘滞系数；
—与额外的动力学量S对应的、多粒子系统的副热浴温度；
Rx(t),Rs(t)—相互独立的白噪声函数；
U(X)—多粒子系统在三维欧氏空间中的粒子间相互作用势能函数，是所有粒子的坐标X的函数；
sj—j粒子额外的动力学量；
Dj(X)—j粒子的近邻偏心度量，是所有粒子的坐标X的函数，由上文中的公式定义；
κ—多粒子系统额外的动力学量S与粒子的近邻偏心度量D之间的力耦合参量；
Uκ(X,S)—多粒子系统在四维空间(见下文)中粒子间相互作用对应的势能函数，是所有粒子的坐标X与额外的动力学量S的函数，由该方程组第三行表达式定义。
这里本发明为多粒子系统中每个粒子赋予了一个额外的动力学量s，所有粒子额外的动力学量s构成N维向量S＝(s1,s2,s3,…,sN)。S与多粒子系统中所有粒子的近邻偏心度量D＝(D1(X),D2(X),D3(X),…,DN(X))相对应，该对应关系的物理意义在于，如果以所有粒子的近邻偏心度量D＝(D1(X),D2(X),D3(X),…,DN(X))作为多粒子动力学系统的集合变量(N维)，则N维向量S＝(s1,s2,s3,…,sN)即为该集合变量对应的热力学与统计力学量。需要指出的是，与粒子额外的动力学量S对应的粘滞系数和热浴温度各自独立于多粒子系统原始动力学量(X,V)对应的粘滞系数γ与热浴温度T，而Rx(t)与Rs(t)也为互相独立的白噪声。表达式中的物理量T与在专利号为202010884514.3的专利方法中被分别称为“多粒子系统的热力学温度”和“与额外的动力学量s相关的虚拟温度”，而在本发明方法中则被分别称为“多粒子系统的主热浴温度”和“与额外的动力学量S对应的、多粒子系统的副热浴温度”，本发明方法中所赋予的名称更符合这两个物理量的本质物理属性。
该方程组的物理意义可以这样理解。粒子额外的动力学量s与作为原始多粒子系统集合变量的粒子近邻偏心度量D具有对应关系，而根据粒子的近邻偏心度量D的定义表达式，D具有与粒子在三维欧氏空间中坐标的坐标轴分量相同的物理量纲，因此动力学量s也与粒子的三维欧氏空间坐标的坐标轴分量具有相同的物理量纲，而赋予每个粒子额外的动力学量s相当于给每个粒子新增了一个额外的动力学自由度，于是可以认为多粒子系统是在一个四维空间中运动。在这个四维空间中，将所有粒子在原三维欧氏空间中的粒子近邻偏心度量D与粒子额外的动力学量S以谐振子的方式进行耦合而得到粒子在四维空间中的粒子间相互作用势能函数Uκ(X,S)。在势函数Uκ(X,S)的作用下，多粒子系统的动力学量X与S分别以相互独立的粘滞系数和热浴温度各自依据普通Langevin动力学方程进行状态时间演化。因此，本发明称上述方程组为四维空间下多粒子系统的扩展Langevin动力学方程组。
SAMD仿真模拟方法的主要任务即是对上述四维扩展Langevin动力学方程组在计算机上进行数值计算求解。可以看出该方程组是一个以时间为自变量的耦合二阶常微分方程组，通过简单的变量替换(如)即可以将其转化为耦合一阶常微分方程组，然后对时间进行数值离散化，设置离散时间步长固定为Δt(也可以设置为非固定的时间步长)，用差分格式替换微分(如)将其转化为普通代数方程组，根据步骤1所给定的(X,V)的初值、(N,M)的值和U(X)函数的数据，以及粒子的近邻偏心度量D的计算公式、参量Rd的值，再给定方程组中的各个参量κ、T、的值，即可以逐步将每个离散时间步(Δt,2Δt,3Δt,…,nΔt)对应的所有粒子的坐标X与速度V的值求解出，给出多粒子凝聚态物质和材料系统的微观状态随时间演化的情况。
对该方程组的具体的计算求解可以有不同的实现和数值处理细节，如Verlet算法、蛙跳(Leap frog)算法、Gear预测-矫正算法等等(参见文献[1,26])，这些具体的数值计算处理细节技术已经很成熟，不再赘述。需要说明的是，虽然上述四维扩展Langevin动力学方程组的具体数值求解可以有不同的处理方式和算法，但只要采用了该方程组按本说明书步骤进行物质和材料的计算机数值仿真模拟，其使用均受本发明专利的保护。
步骤3：根据步骤2所计算出来的多粒子系统在每个离散时间步(Δt,2Δt,3Δt,…,nΔt)对应的(X,V)的值，对相应的凝聚态物质和材料系统进行各种热力学与动力学分析，包括多粒子系统模型的空间构型及其演化分析、系统在一定宏观热力学条件下所发生的内部微观结构转变的过程及其动力学特征与规律，等等。由于这些分析方法、手段和技术可参考传统的普通分子动力学模拟方法所用的方法和手段，本发明不再赘述。
图9提供了专利号为202010884514.3的专利申请《一种用于物质和材料的计算机仿真模拟方法及装置》中的方法流程图，明显可见本发明方法与所引用专利的方法在流程上的不同之处，其中没有进行多粒子系统初始参考状态设定步骤。这是本发明从新的角度提出模型中微观结构转变过程的反应坐标表达方案后，在整个计算机仿真模拟流程上相应提出的要素省略改进。
为便于实施参考起见，以下提供两个具体应用例子。首先按照上述实施方式中所列步骤1-步骤2-步骤3的基本框架(见附图1)对SAMD仿真模拟方法(主要是步骤2中的四维扩展Langevin动力学方程组)进行数值离散化和计算机自动化运算流程处理。然后针对用户任意感兴趣的物质与材料系统，只需按步骤1准备该物质与材料系统的多粒子模型，包括(N,M)的值、系统的粒子(原子或分子)间相互作用势能函数U(X)的数据和系统的微观状态对应的(X,V)的初值，将其作为输入数据；同时提供SAMD仿真模拟方法所需的各种参量的数值，主要包括：
1)动力学方程时间离散化处理所用的时间步长Δt(推荐的取值范围为0.1～5.0fs)；
2)为计算粒子的近邻偏心度量D而构造每个粒子的最近邻粒子集合所需的截断半径Rd(一般可设为多粒子系统径向分布函数曲线的第一个谷底对应的位置，其取值应尽量使得多粒子系统内部每个粒子的最近邻粒子集合里的粒子数目在8～20之间)；
3)四维扩展Langevin动力学方程组参量κ、γ与(推荐的取值范围分别为： 2～20ps-1、1～100ps-1)；
4)四维扩展Langevin动力学方程组参量T(多粒子系统的主热浴温度，根据仿真模拟需要和实际情况来设定)；
5)四维扩展Langevin动力学方程组参量(与粒子额外的动力学量S对应的、多粒子系统的副热浴温度，根据仿真模拟需要和实际情况来设定，一般可取100～10000K范围内的值)；
6)总运行时间步数NRun(根据仿真模拟需要和实际情况来设定，一般取值为100～1×108范围内不等)。
再采用计算机自动化运算流程运行求解每个离散时间步(Δt,2Δt,3Δt,…,NRunΔt)对应的(X,V)的值，即可以给出该物质和材料系统的微观状态随时间演化的情况，实现对该物质和材料系统动力学过程的计算机仿真模拟，最后在此基础上再分析该物质和材料系统内部发生的微观结构转变事件和相应的动力学特征与规律。
在具体应用例子的流程实现过程中，借用了目前最为流行与广泛使用的分子动力学模拟软件LAMMPS(参见文献[2])，利用LAMMPS的开源特性、健壮的代码质量、丰富的功能以及其优良的模块化程序设计框架和便利的扩展接口，可以较为方便地将SAMD仿真模拟方法的计算机自动化运算流程编制成计算机程序代码并融合进LAMMPS软件中，并且为LAMMPS软件提供了两个新的、为SAMD仿真模拟方法特设的LAMMPS输入脚本命令(computedecentr/atom和fix samd/decentr/atom)。在使用时，采用这个经过修改和功能扩充后的LAMMPS软件，在计算机上对其进行源程序代码编译后得到可执行文件，再按照通常的使用LAMMPS软件进行普通分子动力学模拟的步骤，在计算机上执行仿真模拟任务(参见LAMMPS使用手册，文献[2])：设置与编制LAMMPS输入命令行脚本——提供各种必要的输入数据文件——提交LAMMPS可执行程序执行计算——分析输出结果数据文件，只需额外在LAMMPS输入命令行脚本文件中相应地添加和设置这两个为SAMD仿真模拟方法特设的新命令行(compute decentr/atom和fix samd/decentr/atom)，即可进行SAMD仿真模拟。
需要强调的是，虽然在下述具体应用例子中，借用了LAMMPS开源软件，将SAMD仿真模拟方法的计算机自动化运算流程编制成计算机程序代码并融合进LAMMPS软件中，再采用这个经过修改和功能扩充后的LAMMPS软件进行本发明所提出的SAMD仿真模拟方法的实施，但本发明的实施并不受LAMMPS软件的任何限制，任何其它可行的计算机程序编制方案及相应的实施方法都可以对本发明技术方案进行具体实施，所有采用这些计算机程序编制方案对本发明的实施与使用都受本发明专利的保护。
具体应用例子1：铁碳合金中的间隙碳原子和铁原子空位在体心立方铁晶格中的扩散迁移动力学过程
常用钢铁材料最基本的组成即为铁碳合金，由铁原子周期排列成体心立方晶格，而少量碳原子则作为间隙原子位于铁原子晶格的八面体间隙中，如附图2(A)所示。在通常温度(如室温)下，铁原子在其晶格阵点位置作热振动，会有一定概率脱离阵点位置变成自间隙原子，同时留下一个阵点的铁原子空位，如附图2(B)所示。在铁碳合金中，碳原子的扩散迁移与铁原子通过空位进行扩散迁移是重要动力学过程，在钢铁材料表面渗碳处理、钢的热处理引起组织结构转变和力学性能改变等材料工艺中发挥着重要作用，是影响这些工艺过程和材料性能的关键物理化学过程。
步骤1：针对附图2(B)的铁晶体的原子模型(含有1个铁原子空位与3个间隙碳原子)，进行该多粒子(原子)模型的构建和数据准备如下：
1)Fe-C二元合金系的原子间相互作用势能函数U(X)由专门的势函数文件提供(可从网上公开的势函数文件库获取，参见文献[27]，这里将该文件命名为“potential.file”)，势函数文件“potential.file”中包含Fe原子与C原子各自的质量的数据；
2)设置和编写一个LAMMPS输入命令脚本文件，在该文件中加入lattice、region、create_box、create_atoms、delete_atoms等命令设定在模型空间中填充铁原子(共计5400个铁原子)，并设定在铁晶体模型中不同的八面体间隙处添加3个碳原子，同时指定删除模型中1个铁原子使该原子所处位置成为空位，采用pair_style和pair_coeff命令设定模型中Fe-C二元合金系的原子间相互作用势函数文件“potential.file”，再采用minimize、timestep、velocity、fix npt、run等命令指定对初始模型进行弛豫，用write_restart命令设定最后输出结果到文件“model.file”；
3)该输入命令脚本文件生成后即可以传递给LAMMPS可执行程序并将其提交到计算服务器上执行运算，得到相应的输出结果文件“model.file”，该文件中包含模型中总原子数N的值和所有原子的坐标X与速度V的初始值的数据。
步骤2：以步骤1中所得到的输出结果文件“model.file”对应的模型状态为基础，采用LAMMPS软件对该模型进行SAMD仿真模拟，具体实施过程如下：
1)设置和编写一个新的LAMMPS输入命令脚本文件，在该文件中加入read_restart命令读取步骤1得到的输出结果文件“model.file”来导入模型，用pair_style和pair_coeff命令指定模型中Fe-C二元合金系的原子间相互作用势函数文件“potential.file”(与步骤1文件相同)，用timestep、compute decentr/atom、fix langevin、fix samd/decentr/atom等命令指定对模型进行SAMD仿真模拟的各个参量的取值如下表所示：
这里，截断半径Rd的值正好等于体心立方铁原子晶格径向分布函数的第一个谷底的位置，即每个铁原子正常有14个最近邻原子，用compute msd命令指定在模拟过程中计算碳原子的均方根位移，用fix nph命令指定使用Nose-Hoover动力学将模型的压力控制为0应力状态，用thermo命令指定每隔5个时间步输出模型的热力学状态信息(包括前述碳原子的均方根位移)到文件“thermo.file”，用restart命令指定每隔10000个时间步输出模型当前时间步下的各种动力学状态信息(包括所有原子的坐标X与速度V在当前时刻下的值)到文件“kinetics.file”，用于后续分析模型的空间构型随时间演化的情况，用run命令指定SAMD仿真模拟总运行时间步数NRun为2000000。
2)该输入命令脚本文件生成后即可以传递给LAMMPS可执行程序并将其提交到计算服务器上执行运算，得到相应的SAMD仿真模拟的输出结果数据文件“thermo.file”和“kinetics.file”。
步骤3：对步骤2中所得到的输出结果数据(热力学状态信息、动力学状态信息)文件“thermo.file”和“kinetics.file”进行分析。提取热力学状态信息输出文件“thermo.file”中的碳原子的均方根位移和相应的时间步序号信息，用绘图软件给出碳原子均方根位移随时间步数增加变化的情况(见附图3)；对动力学状态信息文件“kinetics.file”所给出的模型的空间构型状态进行原子的热振动噪声消除，再使用OVITO原子构型可视化程序对消除热振动噪声后的模型空间构型进行计算机可视化图像处理与分析，导出模型在不同时间步(nsteps)下的原子空间构型(包括碳原子与铁原子空位)的图像：(1)nsteps＝0(附图4(A))；(2)nsteps＝5e+05(附图4(B))；(3)nsteps＝1e+06(附图4(C))；(4)nsteps＝1.5e+06(附图4(D))；(5)nsteps＝2e+06(附图4(E))。
本具体应用例子的步骤1、2中，除了新添加的compute decentr/atom和fix samd/decentr/atom这两个为SAMD仿真模拟方法特设的命令外，所有涉及到的LAMMPS命令的使用与设定均按照LAMMPS使用手册的说明(参见文献[2])以常规方式进行，在此一些相关细节不再赘述；这些命令都可以在LAMMPS使用手册中找到详细的说明(参见文献[2])，在此也不再对其作进一步的解释。
需要说明的是，对该模型中Fe-C二元合金系的原子间相互作用势函数，采用的是C.S.Becquart等人开发的针对Fe-C二元系的嵌入原子方法型原子间相互作用势(参见文献[28,29])，经过验证表明该势函数可以比较真实地反映Fe-C合金的基本物理特性，具有较好的可靠性与合理性。经计算可得，该势函数所给出的间隙碳原子与铁原子空位扩散激活能垒分别为0.81eV/atom和0.63eV/atom，与实验测算结果比较接近。
从附图3中可以看出，采用本发明所提出的SAMD仿真模拟方法对附图2(B)所示的模型进行计算机数值仿真模拟，当系统的副热浴温度时，在总时间步数为2.0×106的模拟时间范围内，未发现碳原子或铁原子空位的迁移运动；从附图3和附图4(A)-(E)上可以看出，随着系统的副热浴温度从3000K增加到5000K，碳原子的迁移频率和迁移距离会逐渐增加，表明通过调节系统的副热浴温度参量的值可以有效地控制本发明所提出的SAMD仿真模拟方法在时间尺度上的扩展与跨越的程度。
经过简单估算可得，对于系统的副热浴温度分别为3000K、4000K、5000K的SAMD仿真模拟，在总时间步数为2.0×106的模拟时间范围内，模型中碳原子的平均迁移跳跃次数分别约为4.39、26.67、61.28。实验测算分析结果表明，热力学温度为300K下碳原子在相邻八面体间隙位置间平均每秒跳动次数约为2次。通过对比分析可得，在本具体应用例子中，对于系统的副热浴温度分别为3000K、4000K、5000K的SAMD仿真模拟，其相应的时间加速比分别约为1.1×109、6.7×109、1.5×1010。同时，分析也可以发现，在本具体应用例子所针对的模型中(附图2)，间隙碳原子与铁原子空位作为迁移激活能较为接近的两种不同点缺陷，采用本发明所提出的SAMD仿真模拟方法能较好地反映出两者的相对激活频率(见附图4)，表明本发明所提出的SAMD仿真模拟方法可以较好地保持物质和材料中不同结构转变事件和过程的动力学相对概率，具有较好的物理动力学保真性。
具体应用例子2：在溶解有一定浓度氢原子的金属铝双晶体中，氢原子在铝双晶体中扩散以及在晶界上偏聚的动力学过程
自然界与工程环境中会存在各种各样含氢元素的物质，如广泛存在的水，还有各种酸、碱、盐溶液(如海水)，金属结构件在服役过程中与这些含氢物质接触后，由于表面化学或电化学反应，氢原子会渗透进金属中，弱化金属原子间结合强度，从而产生所谓的“氢脆”，给金属结构件的服役和工程设备的使用带来严重的安全隐患。由于氢原子是尺寸最小的原子，其在金属材料中的扩散一般较为容易，同时实验研究表明，溶解在金属中的氢原子一般会在晶体缺陷如位错、晶界上偏聚，特别是在晶界上的偏聚往往会造成晶界的显著弱化，成为裂纹萌生的源头和扩展的择优路径。因此，深入理解与认识氢原子在铝晶体中的扩散与偏聚行为对于研究金属铝的氢脆行为有重要的参考意义。
步骤1：针对附图5所示的溶解有一定浓度氢原子的金属铝双晶体(所含晶界为∑5<110>{310}对称倾斜晶界)的原子模型，首先采用LAMMPS软件进行该多粒子(原子)模型的构建和数据准备如下：
1)Al-H二元合金系的原子间相互作用势能函数U(X)由专门的势函数文件提供(可从网上公开的势函数文件库获取，参见文献[27]，这里将该文件命名为“potential.file”)，势函数文件“potential.file”中包含Al原子与H原子各自的质量的数据；
2)设置和编写一个LAMMPS输入命令脚本文件，在该文件中加入lattice、region、create_box、create_atoms、delete_atoms等命令指定按晶格类型与晶粒取向在模型空间中填充铝原子与氢原子(总计5664个铝原子与200个氢原子)，并指定删除相互间距离太靠近的原子对中的其中一个原子，用pair_style和pair_coeff命令指定模型中Al-H二元合金系的原子间相互作用势函数文件“potential.file”，再采用minimize、timestep、velocity、fix npt、run等命令指定对初始模型进行弛豫，用write_restart命令指定最后输出结果到文件“model.file”；
3)该输入命令脚本文件生成后即可以传递给LAMMPS可执行程序并将其提交到计算服务器上执行运算，得到相应的输出结果文件“model.file”，该文件中包含模型中总原子数N的值和所有原子的坐标X与速度V的初始值的数据，最后模型中总计有5846个原子(5664个铝原子和182个氢原子)。
步骤2：以步骤1中所得到的输出结果文件“model.file”对应的模型状态为基础，采用LAMMPS软件对该模型进行SAMD仿真模拟，具体实施过程如下：
1)设置和编写一个新的LAMMPS输入命令脚本文件，在该文件中加入read_restart命令读取步骤1得到的输出结果文件“model.file”来导入模型，用pair_style和pair_coeff命令指定模型中Al-H二元合金系的原子间相互作用势函数文件“potential.file”(与步骤1文件相同)，用timestep、compute decentr/atom、fix langevin、fix samd/decentr/atom等命令指定对模型进行SAMD仿真模拟的各个参量的取值如下表所示：
这里，截断半径Rd的值略小于面心立方铝原子晶格径向分布函数的第一个谷底的位置，用compute msd命令指定在模拟过程中计算氢原子的均方根位移，用fix nph命令指定使用Nose-Hoover动力学将模型的压力控制为0应力状态，用thermo命令指定每隔5个时间步输出模型的热力学状态信息(包括前述氢原子的均方根位移)到文件“thermo.file”，用restart命令指定每隔10000个时间步输出模型当前时间步下的各种动力学状态信息(包括所有原子的坐标X与速度V在当前时刻下的值)到文件“kinetics.file”，用于后续分析模型的空间构型随时间演化的情况，用run命令指定SAMD仿真模拟总运行时间步数NRun为2000000。
2)该输入命令脚本文件生成后即可以传递给LAMMPS可执行程序并将其提交到计算服务器上执行运算，得到相应的SAMD仿真模拟的输出结果数据文件“thermo.file”和“kinetics.file”。
步骤3：对步骤2中所得到的输出结果数据(热力学状态信息、动力学状态信息)文件“thermo.file”和“kinetics.file”进行分析。提取热力学状态信息输出文件“thermo.file”中的氢原子的均方根位移和相应的时间步序号信息，用绘图软件给出氢原子均方根位移随时间步数增加变化的情况(见附图6)；对动力学状态信息文件“kinetics.file”所给出的模型的空间构型状态进行原子的热振动噪声消除，再使用OVITO原子构型可视化程序对消除热振动噪声后的模型空间构型进行计算机可视化图像处理与分析，导出模型在不同时间步(nsteps)下的原子空间构型(隐藏铝原子)的图像：(1)nsteps＝0(附图7(A))；(2)nsteps＝5e+05(附图7(B))；(3)nsteps＝1e+06(附图7(C))；(4)nsteps＝1.5e+06(附图7(D))；(5)nsteps＝2e+06(附图7(E)，附图8(A)-(C))。
本具体应用例子的步骤1、2中，除了新添加的compute decentr/atom和fix samd/decentr/atom这两个为SAMD仿真模拟方法特设的命令外，所有涉及到的LAMMPS命令的使用与设定均按照LAMMPS使用手册的说明(参见文献[2])以常规方式进行，在此一些相关细节不再赘述；这些命令都可以在LAMMPS使用手册中找到详细的说明(参见文献[2])，在此也不再对其作进一步的解释。
需要说明的是，对该模型中Al-H二元合金系的原子间相互作用势函数，采用的是F.Apostol和Y.Mishin所开发的针对Al-H二元系的键角相关型多体原子间相互作用势(参见文献[30])，经过验证表明该势函数可以比较真实地反映Al-H合金的基本物理特性，具有较好的可靠性与合理性。经计算可得，该势函数能正确地给出氢原子在铝晶格中的四面体间隙位置相对于八面体间隙位置的能量较低，所给出的间隙氢原子的扩散激活能垒为0.19eV/atom，与第一性原理计算和实验测算结果很接近，同时也表明氢原子在金属铝中扩散较为容易。
从附图6中可以看出，采用本发明所提出的SAMD仿真模拟方法对附图5所示的模型进行计算机数值仿真模拟，随着系统的副热浴温度从500K上升到1500K，氢原子的迁移频率会逐渐增加并产生相当程度的迁移扩散距离(如附图7(A)-(E)所示)，表明通过调节系统的副热浴温度参量的值可以有效地控制本发明所提出的SAMD仿真模拟方法在时间尺度上的放大与跨越的程度。从附图7(A)-(E)上也可以看出，随模拟时间的增加，氢原子发生扩散迁移，从最初的较为分散的分布状态(附图7(A))，逐渐在晶界上产生氢原子偏聚的结果(附图7(B)-(E))。附图8(A)-(C)所给出的结果也表明，只要氢原子的迁移频率足够高(如增加模拟时间步数或适当提高系统的副热浴温度值)，本发明所提出的SAMD仿真模拟方法都能给出氢原子在晶界上偏聚的结果，这与由蒙特卡罗(MC)方法模拟所给出的结果(附图8(D))比较一致。氢原子在晶界上产生偏聚的主要原因是氢原子在铝晶格中的扩散迁移与其在晶界上及晶界附近进行来回扩散迁移的激活能垒有所不同，而本发明所提出的SAMD仿真模拟方法能较好地给出氢原子在晶界上偏聚的结果，表明本发明所提出的SAMD仿真模拟方法很好地实现了针对具有不同激活能垒和反应系数的结构转变事件与过程，赋予其相应的、具有合适大小的激活频率，从而较好地保证了物质与材料中结构转变动力学的物理真实性。
经过简单估算可以得到，对于系统的副热浴温度分别为500K、1000K、1500K的SAMD仿真模拟，在总时间步数为2.0×106的模拟时间范围内，模型中氢原子的平均迁移跳跃次数分别约为13.0、58.3、138.5，而热力学温度为100K下氢原子在铝晶体中扩散系数实验测算值约为5.98×10-17m2/s，相当于氢原子在铝晶体四面体间隙中平均每秒跳动约8750次，从而可以估算出本具体应用例子中，相应的SAMD仿真模拟的时间加速比分别约为1.5×106、6.5×106、1.6×107。
为便于理解本发明和现有技术的效果差异，通过实验：(1)对比两者在进行技术实现时，同样借用LAMMPS开源软件，将SAMD仿真模拟方法的计算机自动化运算流程编制成计算机程序代码并融合进LAMMPS软件中，进行计算机程序代码编制所需的代码行数；(2)对比两者在同样利用某计算服务器，各自进行上述具体应用例子1(使用8个CPU核心并行计算)和具体应用例子2(使用16个CPU核心并行计算)的实施过程中，进行单次计算得出数据结果所需要消耗的时间；提供相应数据如下：
可见，相比于专利号为202010884514.3的专利方法，本发明的优点在于，由于在本发明方法中所设置的粒子近邻偏心度量D的计算不需要进行多粒子系统参考状态的设定，且其计算表达式相比于专利号为202010884514.3的专利方法中所设置的粒子曳步运动度量D的计算表达式要更为简化，因此在实际的计算机仿真模拟执行中，本发明方法更易于进行计算机程序编制，同时在计算效率上也大幅提高。另外，由于实施流程和步骤上的简化，本发明方法在实施使用时相对于专利号为202010884514.3的专利所给出的方法要更为便捷易用。
具体实施时，本发明技术方案提出的方法可由本领域技术人员采用计算机软件技术实现自动运行流程，实现方法的系统装置例如存储本发明技术方案相应计算机程序的计算机可读存储介质以及包括运行相应计算机程序的计算机设备，也应当在本发明的保护范围内。
在一些可能的实施例中，提供一种用于物质和材料的计算机仿真模拟装置，包括处理器和存储器，存储器用于存储程序指令和相关输入与输出数据，处理器用于调用处理器中存储的指令执行如上所述的一种用于物质和材料的计算机仿真模拟方法。
在一些可能的实施例中，提供一种用于物质和材料的计算机仿真模拟装置，包括可读存储介质，所述可读存储介质上存储有计算机程序，所述计算机程序执行时，实现如上所述的一种用于物质和材料的计算机仿真模拟方法。
以上所述仅为使用本发明所给出的具体实施例，并不用于限制本发明。凡在本发明的指导思想、精神与原则之内，所做的任何修改、改进等，均应包含在本发明专利的保护范围内。