MS分子动力学模拟具体实施步骤

分子动力学模拟方法Molecular Dynamics Simulation Method分子动力学模拟方法是一种计算方法，可以预测原子和分子在不同温度和压力下的运动和力学行为。

该方法已被广泛应用于物理、化学、生物学和材料科学等领域，用于研究材料性质、生物分子结构和动态、相变等现象。

本文将介绍分子动力学模拟的基本原理、模拟过程以及如何用该方法研究材料或生物分子。

1. 基本原理分子动力学模拟基于牛顿力学原理，用原子和分子之间的势能函数描述系统内部的相互作用力。

根据牛顿第二定律 F=ma，通过求解系统中每个分子的运动方程来推导出分子的运动轨迹。

在计算中，采用的势能函数决定了分子之间的相互作用，包括范德华力、静电作用、键角等力。

基于这些相互作用力和分子的运动轨迹，可以计算出分子的位置、速度、加速度和能量等物理量。

2. 模拟过程分子动力学模拟的过程包括初始化、模拟和分析三个阶段。

2.1 初始化初始化阶段主要是为模拟设置一些参数，包括分子数、模拟时间、初速度、初位置和系统温度等。

初速度可以根据玻尔兹曼分布生成，初位置随机分布，系统温度也可以通过控制分子初速度实现。

模拟阶段分为两个步骤：计算分子运动和更新分子位置。

计算分子运动：在每个时间步中，使用牛顿运动方程计算每个分子的运动。

分子与其他分子之间的相互作用通过势能函数计算。

时间步长各不相同，一般为1-10飞秒。

更新分子位置：根据计算出的分子运动轨迹和速度，使用欧拉法更新分子位置。

在此过程中，通过周期性边界条件保证系统的连续性。

2.3 分析分析阶段主要是对模拟结果进行分析和处理，如计算能量、相变、速度相关的分布函数等。

有效的分析可以给出关键参数和物理量，如分子动力学能量、热力学性质和动力学行为。

3. 应用分子动力学模拟方法已经被广泛应用于物理、化学、生物学和材料科学等研究领域，尤其是材料和生物分子方面的研究具有广泛的前景。

3.1 材料科学分子动力学模拟可用于研究材料的力学、热力学和电学等性质。

分子动力学模拟流程Molecular dynamics simulation is a powerful tool in the field of computational chemistry and physics. It allows researchers to study the movement of atoms and molecules over time, providing valuable insights into the behavior of materials at the molecular level. By simulating the interactions between particles based on classical mechanics, scientists can explore various physical and chemical processes in great detail.分子动力学模拟是计算化学和物理领域中一种强大的工具。

它允许研究人员随着时间的推移研究原子和分子的运动，为在分子水平上材料的行为提供有价值的见解。

通过基于经典力学对粒子之间的相互作用进行模拟，科学家可以详细地探索各种物理和化学过程。

One of the key advantages of molecular dynamics simulation is its ability to capture the dynamics of complex systems that are difficult to study experimentally. By monitoring the trajectories of individual particles in a simulated environment, researchers can observe how macroscopic properties emerge from the interactions of atoms and molecules. This information is crucial for understanding the behaviorof materials under different conditions and for designing new materials with desired properties.分子动力学模拟的一个关键优势是它能够捕获实验难以研究的复杂系统的动态。

第3章铁基块体非晶合金-纳米晶转变的动力学模拟过程3.1 Discover模块3.1.1 原子力场的分配在使用Discover模块建立基于力场的计算中，涉及几个步骤。

主要有：选择力场、指定原子类型、计算或指定电荷、选择non-bond cutoffs。

在这些步骤中，指定原子类型和计算电荷一般是自动执行的。

然而，在某些情形下需要手动指定原子类型。

原子定型使用预定义的规则对结构中的每个原子指定原子类型。

在为特定的系统确定能量和力时，定型原子使工作者能使用正确的力场参数。

通常，原子定型由Discover使用定型引擎的基本规则来自动执行，所以不需要手动原子定型。

然而，在特殊情形下，人们不得不手动的定型原子，以确保它们被正确地设置。

图 3-11）计算并显示原子类型：点击Edit→Atom Selection，如图3-1所示图3-2弹出对话框，如图3-2所示从右边的…的元素周期表中选择Fe，再点Select，此时所建晶胞中所有Fe 原子都将被选中，原子被红色线圈住即表示原子被选中。

再编辑集合，点击Edit →Edit Sets，如图3-3、3-4所示。

图3-3图3-4弹出对话框见图3-4，点击New...，给原子集合设定一个名字。

这里设置为Fe，则3D视图中会显示“Fe”字样，再分配力场：在工具栏上点击Discover按钮，从下拉列表中选择Setup，显示Discover Setup对话框，选择Typing选项卡，见图3-5。

图3-5在Forcefield types里选择相应原子力场，再点Assign（分配）按钮进行原子力场分配。

注意原子力场中的价态要与Properties Project里的原子价态（Formalcharge）一致。

3.1.2力场的选择1）Energy，见图3-6。

图3-6力场的选择：力场是经典模拟计算的核心，因为它代表着结构中每种类型的原子与围绕着它的原子是如何相互作用的。

对系统中的每个原子，力场类型都被指定了，它描述了原子的局部环境。

材料计算设计基础实验指导书朱景川编哈尔滨工业大学2005年2月实验一、第一性原理计算1. 实验目的(1) 掌握第一性原理和密度泛涵的计算方法；(2) 学会使用Visualizer 的各种建模和可视化工具；(3) 熟悉CASTEP 模块的功能。

2. 实验原理CASTEP 是基于密度泛涵理论平面波赝势基础上的量子力学计算。

密度泛涵理论的基本思想是原子、分子和固体的基本物理性质可以用粒子密度函数进行描述。

可以归纳为两个基本定理：定理1：粒子数密度函数是一个决定系统基态物理性质的基本参量。

定理2：在粒子数不变的条件下能量对密度函数变分得到系统基态的能量。

不计自旋的全同费米子的哈密顿量为：H T U V =++其中动能项为：()()T dr r r ψψ+=∇∇⎰库仑作用项为：11'()(')()(')2'U drdr r r r r r r ψψψψ++=-⎰ V 为对所有粒子均相同的局域势u(r)表示的外场影响：()()()V dru r r r ψψ+=⎰粒子数密度函数为：()()()r r r ρψψ+=ΦΦ对于给定的()r υ，能量泛函[]E ρ定义为：[]()()E dr r r T U ρυρ=+Φ+Φ⎰；[]F T U ρ=Φ+Φ系统基态的能量：'''''[]''''[][]()()[][]()()[]E T U V GE F dr r r E G G F dr r r E G ρρυρφρυρρΦ=Φ+Φ+ΦΦ==+>⎰=+=⎰3. 实验内容实验 1. 材料的电子结构计算；实验2. 晶体材料的晶格[点阵]参数预报（要求材料体系为金属合金、化合物半导体或有机 高分子材料）；实验 3. 材料的弹性模量计算。

* 在三个实验内容中可以任选一个内容进行计算，有能力的同学也可以全做。

4. 实验设备和仪器(1)硬件：多台PC机和一台高性能计算服务器。

第3章铁基块体非晶合金-纳米晶转变的动力学模拟过程3.1 Discover模块3.1.1 原子力场的分配在使用Discover模块建立基于力场的计算中，涉及几个步骤。

主要有：选择力场、指定原子类型、计算或指定电荷、选择non-bond cutoffs。

在这些步骤中，指定原子类型和计算电荷一般是自动执行的。

然而，在某些情形下需要手动指定原子类型。

原子定型使用预定义的规则对结构中的每个原子指定原子类型。

在为特定的系统确定能量和力时，定型原子使工作者能使用正确的力场参数。

通常，原子定型由Discover使用定型引擎的基本规则来自动执行，所以不需要手动原子定型。

然而，在特殊情形下，人们不得不手动的定型原子，以确保它们被正确地设置。

图 3-11）计算并显示原子类型：点击Edit→Atom Selection，如图3-1所示图3-2弹出对话框，如图3-2所示从右边的…的元素周期表中选择Fe，再点Select，此时所建晶胞中所有Fe 原子都将被选中，原子被红色线圈住即表示原子被选中。

再编辑集合，点击Edit →Edit Sets，如图3-3、3-4所示。

图3-3图3-4弹出对话框见图3-4，点击New...，给原子集合设定一个名字。

这里设置为Fe，则3D视图中会显示“Fe”字样，再分配力场：在工具栏上点击Discover按钮，从下拉列表中选择Setup，显示Discover Setup对话框，选择Typing选项卡，见图3-5。

图3-5在Forcefield types里选择相应原子力场，再点Assign（分配）按钮进行原子力场分配。

注意原子力场中的价态要与Properties Project里的原子价态（Formalcharge）一致。

3.1.2力场的选择1）Energy，见图3-6。

图3-6力场的选择：力场是经典模拟计算的核心，因为它代表着结构中每种类型的原子与围绕着它的原子是如何相互作用的。

对系统中的每个原子，力场类型都被指定了，它描述了原子的局部环境。

分子动力学模拟的研究与应用分子动力学模拟是一种通过计算机模拟分子间相互作用以及它们随时间的演变而得出分子结构和行为的方法。

这种方法对于研究物质的性质和反应、设计新的材料和药物、优化现有材料和生产工艺等方面具有广泛的应用。

一、分子动力学模拟的原理与流程分子动力学模拟基于牛顿力学原理，将物质看作是由分子组成的，通过计算其分子之间的相互作用力和速度的数值解，来推断研究物质的宏观行为。

具体地说，分子动力学模拟流程包括以下几个步骤：1. 确定所研究物质的分子模型和初始结构。

通常情况下，这需要使用计算化学方法先预测分子的理论结构，并采取实验手段进行验证。

2. 确定模拟中所使用的计算力场。

力场是描述原子或分子间相互作用的一组公式和参数，决定了分子间的势能和力的大小和方向。

3. 设定模拟的温度、压力等条件。

模拟温度和压力等参数的设定需根据实验需求进行选择。

4. 计算分子间的相互作用力和速度。

分子之间的相互作用力通常采用范德华力、库伦力和键能等几种力。

计算相互作用力和速度的数值解可以通过使用微分方程和数值计算方法得到。

5. 根据计算得到的力和速度，重新计算分子的位置和速度。

这通过使用牛顿运动学公式得到。

6. 重复步骤4和5，直至模拟结束或达到设定的模拟时间，得到分子结构和其运动的相关信息。

二、分子动力学模拟的应用领域分子动力学模拟在许多领域中得到了广泛的应用，以下列举几个例子：1. 材料设计与性能优化。

通过在计算机中模拟、预测物质的性质和反应，可以用于设计新的材料，提高材料的性能和品质，减少生产成本，优化生产工艺等。

2. 药物研发与医学应用。

分子动力学模拟可用于模拟分子与蛋白质的相互作用，以研究药物在体内产生的效果和活性，从而设计更安全和有效的药物。

3. 天然资源的研究与利用。

运用分子动力学模拟对天然物质进行研究，如植物中的天然化合物，有助于推广利用天然资源提高生产效益。

4. 环境保护与能源研究。

利用分子动力学模拟计算大气、水体和岩石等的化学反应，有助于预测大气污染、水污染和化学品泄漏等环境问题的产生和发展趋势。

分子动力学模拟步骤英文回答：Molecular dynamics simulation is a computational method used to study the motion and behavior of atoms and molecules over time. It is widely used in various fields of science, such as chemistry, physics, and materials science. The simulation involves solving the equations of motion for each particle in the system, taking into account the forces acting on them.To perform a molecular dynamics simulation, several steps need to be followed. Firstly, the system under investigation needs to be defined. This includes specifying the number of particles, their initial positions, and velocities. Additionally, the boundary conditions and any external forces or constraints should be set.Once the system is defined, the next step is to choose an appropriate force field. A force field is a mathematicalmodel that describes the interactions between particles in the system. It includes terms for bonded and non-bonded interactions, such as bond stretching, angle bending, and van der Waals forces. The force field parameters need to be selected based on experimental data or theoretical calculations.After defining the system and force field, the simulation can be carried out. The equations of motion, typically Newton's equations, are numerically integrated to calculate the positions and velocities of the particles at each time step. This is done using algorithms such as the Verlet algorithm or the leapfrog algorithm. The time step size should be carefully chosen to ensure numericalstability and accuracy.During the simulation, various properties of the system can be calculated and analyzed. These include energy, temperature, pressure, and diffusion coefficients. Trajectories of individual particles can also be visualized to gain insights into their motion and interactions.Finally, the simulation results need to be analyzed and interpreted. This involves comparing the simulated properties with experimental data or theoretical predictions. It may also involve statistical analysis and data visualization techniques to identify trends and patterns in the data.In summary, the steps involved in a molecular dynamics simulation include defining the system, selecting a force field, performing the simulation, analyzing the results, and interpreting the findings. This computational method provides valuable insights into the behavior of atoms and molecules, complementing experimental and theoretical approaches.中文回答：分子动力学模拟是一种计算方法，用于研究原子和分子随时间的运动和行为。

分子动力学模拟方法与应用随着科技的不断发展，分子动力学模拟成为了研究和探索分子结构、运动和相互作用的一种重要手段。

和其他物理化学方法相比，分子动力学模拟具有高效、直观、精度高等优点，因此在各种化学、材料科学、生命科学等领域都有广泛应用。

下面我们来详细了解一下分子动力学模拟的方法和应用。

一、分子动力学模拟方法分子动力学模拟是通过计算机模拟分子间相互作用及其运动、运动轨迹、物理和化学性质等的一种方法。

通俗点说，就是通过自动模拟计算机上的分子交互运动，来预测分子在各种条件下的行为及其性质。

分子动力学模拟通常包括以下步骤：1. 构建模型首先，需要根据实验数据或理论推论来构建模型，包括分子的形状、大小、电荷等参数。

2. 选择力场力场是分子动力学模拟的一个重要组成部分，它是描述分子间相互作用的数学模型，包括键的形成和断裂、角度、二面角、静电作用、范德华力等。

目前常用的力场主要有CHARMM、AMBER、GROMOS等。

3. 定义模拟条件分子动力学模拟需要定义一系列条件来确定模拟的前提和结果。

比如温度、压力、模拟时间、模拟步长、边界条件等等。

4. 运行模拟运行模拟后，可以得到分子的各种性质和动态图像等结果，可以通过可视化软件实现对模拟结果的分析和展示。

二、分子动力学模拟的应用1. 生物医学领域在生物医学领域，分子动力学模拟可以用于研究蛋白质、核酸及其他生物分子的结构和功能。

它可以帮助科学家们预测蛋白质和核酸的折叠状态、解析分子间的相互作用、预测系统的稳定状态和反应路径等等。

2. 材料科学领域在材料科学领域，分子动力学模拟可以用于研究材料的物理和化学性质及其行为。

比如，研究纳米材料的表面形态、热力学性质和崩裂机理，预测热膨胀系数、熔点等组成。

此外，分子动力学模拟还可以辅助研究材料的加工和改变制造条件等。

3. 化学反应研究分子动力学模拟还可以用于反应研究。

它可以模拟化学反应中物质的结构和动态，提供反应路径和反应活化能等有用信息。

lammps msd 单位-回复LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) 是一种优秀的分子动力学模拟软件，可用于模拟材料的不同性质和行为。

在分子动力学模拟中，一个重要的参数是Mean Squared Displacement (MSD)（均方位移），它用于描述原子或分子在一定时间范围内的位移，进而提供有关材料物理特性的重要信息。

首先，让我们了解一下MSD 的单位。

在LAMMPS 中，默认长度单位是埃（Angstrom），时间单位是飞秒（femtosecond）。

所以，MSD 的单位是埃的平方（Angstrom squared，简写为"A^2"）。

接下来，我们将逐步介绍如何使用LAMMPS 计算和分析MSD。

第一步，准备LAMMPS 输入文件：要进行MSD 计算，首先需要准备一个输入文件（例如：input.in）来描述模拟系统和模拟参数。

该文件应包括分子结构的定义、原子的势能和相互作用以及模拟的时间范围和步长。

第二步，运行LAMMPS 模拟：打开命令行界面，导航到包含输入文件和相关数据的目录，并运行LAMMPS：lmp_serial -in input.in或者如果你的计算机支持并行计算：lmp_mpi -in input.inLAMMPS 将使用输入文件中的设置和参数自动运行模拟，并计算原子的位置和速度。

第三步，分析模拟结果：完成模拟后，将生成一个或多个轨迹文件，其中包含每个原子的位置和速度信息。

在LAMMPS 中，最常见的轨迹文件格式是随时间变化的XYZ 文件。

使用一种分析工具（如Python 或Matlab）打开轨迹文件，并计算MSD。

MSD 的计算公式如下：MSD(t) = 1/N * Σ[i=1,N] (r(i,t) - r(i,0))^2其中，MSD(t) 是时间t 的MSD，N 是原子数量，r(i,t) 是原子i 在时间t 的位置，r(i,0) 是原子i 的初始位置。

第3章 铁基块体非晶合金‐纳米晶转变的动力学模拟过程 3.1 Discover模块3.1.1 原子力场的分配在使用Discover模块建立基于力场的计算中，涉及几个步骤。

主要有：选择力场、指定原子类型、计算或指定电荷、选择non‐bond cutoffs。

在这些步骤中，指定原子类型和计算电荷一般是自动执行的。

然而，在某些情形下需要手动指定原子类型。

原子定型使用预定义的规则对结构中的每个原子指定原子类型。

在为特定的系统确定能量和力时，定型原子使工作者能使用正确的力场参数。

通常，原子定型由Discover使用定型引擎的基本规则来自动执行，所以不需要手动原子定型。

然而，在特殊情形下，人们不得不手动的定型原子，以确保它们被正确地设置。

图 3-11）计算并显示原子类型：点击Edit→Atom Selection，如图3‐1所示图3-2弹出对话框，如图3‐2所示从右边的…的元素周期表中选择Fe，再点Select，此时所建晶胞中所有Fe原子都将被选中，原子被红色线圈住即表示原子被选中。

再编辑集合，点击Edit →Edit Sets，如图3‐3、3‐4所示。

图3-3图3-4弹出对话框见图3‐4，点击New...，给原子集合设定一个名字。

这里设置为Fe，则3D视图中会显示“Fe”字样，再分配力场：在工具栏上点击Discover按钮，从下拉列表中选择Setup，显示Discover Setup对话框，选择Typing选项卡，见图3‐5。

图3-5在Forcefield types里选择相应原子力场，再点Assign（分配）按钮进行原子力场分配。

注意原子力场中的价态要与Properties Project里的原子价态（Formalcharge）一致。

3.1.2力场的选择1）Energy，见图3‐6。

图3-6力场的选择：力场是经典模拟计算的核心，因为它代表着结构中每种类型的原子与围绕着它的原子是如何相互作用的。

对系统中的每个原子，力场类型都被指定了，它描述了原子的局部环境。

一、实验名称：分子动力学模拟二、实验目的：利用分子动力学理论方法对研究的体系和问题进行理论数值模拟，为实验提供理论参考。

三、实验要求:1.学习并掌握分子动力学的基本理论内容2.掌握嵌入原子理论的主要内容3.对Linux操作系统有足够熟悉的使用经验4.学会分子动力学软件——Lammps的使用方法5.掌握Materials studio的模型搭建方法6.了解VMD的相关操作四、实验相关环境及软件：1.Linux操作系统：一种自由和开放源码的类操作系统，以高效性和灵活性著称，具有多任务、多用户的能力mmps分子模拟软件：即Large-scale Atomic/Molecular Massively Parallel Simulator，主要用于分子动力学相关的一些计算和模拟工作，一般来讲，分子动力学所涉及到的领域，LAMMPS代码也都涉及到了。

3.Materials studio：专门为材料科学领域研究者开发的一款可运行在PC 上的模拟软件，可以帮助解决化学、材料工业中的一系列重要问题。

4.VMD: 一个分子可视化程序，采用3D图形以及内置脚本来对大型生物分子系统进行显示、制成动画以及分析等操作。

五、实验原理：分子动力学（MD）模拟就是用计算机方法表示统计力学，作为实验的一个辅助手段，MD模拟用来研究不能用解析方法来解决的经典多体体系的平衡性质和传递性质，是搭接理论和实验的一个桥梁。

MD方法是一种定性方法，是按照该体系内部的内禀动力学规律来确定位移的转变，追踪系综中每个粒子的个体运动，然后根据统计物理规律给出微观量与宏观量的关系来研究材料性能的一种方法。

目前，该技术已成功地应用于研究晶格畸变，晶粒生长，拉应力-应变关系，蠕变行为，扩散，沉积等。

MD模拟与真实实验在许多方面非常相似，当做一个真实实验时，需要以下步骤: 准备试样; 将试样放在测试仪器中进行测量; 分析测量的结果。

而MD模拟应遵从同样方法: 首先准备试样，即建立一个由 N个粒子组成的模型体系; 将试样放入测试仪器中测量，解N个粒子组成的体系模型的牛顿运动方程，直至平衡; 平衡后，进行材料性质的计算; 最后，进行测量结果的分析，对模拟结果进行分析。

分子动力学的模拟程序设计一、课题名称：分子动力学的模拟程序设计二、班级和姓名：***三、主要内容：1.研究的内容和算法分子动力学模拟方法–确定性的演化过程。

按该体系内部的内禀动力学规律计算、确定位形的转变。

出发点：物理系统确定的微观描述。

具体的说就是描述系统的哈密顿、拉格朗日或者牛顿运动方程，用这些方程去驱动粒子的位置、速度和取向随时间的演化。

模拟过程：(1)通过对物理体系的微观数学描述建立一组方程组，直接求解每个分子的运动方程。

得到每个时刻、每个分子的坐标和动量。

(2)利用统计方法计算系统的静态和动态过程。

计算机模拟的问题1 观测时间是有限的对于某些问题，有限观测时间可以看成是无限长的。

如对一个分子系统的计算，它的观测时间远大于分子时间尺度。

2 观测体系是有限的。

引入周期性、全反射、漫反射等边界条件。

分子动力学模拟的计算对体系的分子动力学方程组进行求解时，需要将运动的连续性方程离散化变成有限差分方程，常用的方法有：欧拉法，龙格-库塔法、辛普生法等。

误差来源：(1) 动力学模型的近似程度；(2) 数值求解法的近似阶数；(3) 数值计算的舍入误差。

分子动力学的适用范围原则上，分子动力学方法所适用的微观物理体系并无任何限制。

这个方法适用的体系既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体既可以是分子，也可以是其它的微观粒子。

微正则系综(NVE)的分子动力学模拟步骤如下：(1) 给定初始空间位置 (2) 计算第n 步时粒子所受的力 (3) 计算第n+1步粒子的空间位置 (4) 计算第n 步的速度 ，根据速度得出系统的温度。

(5) 返回第(2)步，开始下一步计算。

2.源程序1) 简立方结构/*#include <cmath>#include <cstdlib>)}(),{()1()0(ii r r )()(n i n t F F i =m h F r r r n n n n ii i i /22)()1()()1(+-=-+h r r v n n n i i i 2/)()1()1()(-+-=#include <fstream>#include <iostream>#include <string>*/#include <iostream>#include <string.h>#include <math.h>#include <iomanip>#include <stdlib.h>#include <fstream>#include <cstdio>using namespace std;const int N = 64; // number of particlesdouble r[N][3]; // positionsdouble v[N][3]; // velocitiesdouble a[N][3]; // accelerationsdouble L = 10; // linear size of cubical volume，简单立方double vMax = 0.1; // maximum initial velocity componentvoid initialize(){// initialize positionsint n = int(ceil(pow(N, 1.0/3))); // number of atoms in each direction double a = L / n; // lattice spacing,晶格常数int p = 0; // particles placed so farfor (int x = 0; x < n; x++)for (int y = 0; y < n; y++)for (int z = 0; z < n; z++){if (p < N){r[p][0] = (x + 0.5) * a;r[p][1] = (y + 0.5) * a;r[p][2] = (z + 0.5) * a;}++p;}// initialize velocitiesfor (int p = 0; p < N; p++)for (int i = 0; i < 3; i++)v[p][i] = vMax * (2 * rand() / double(RAND_MAX) - 1);}void computeAccelerations(){for (int i = 0; i < N; i++) // set all accelerations to zero for (int k = 0; k < 3; k++)a[i][k] = 0;for (int i = 0; i < N-1; i++) // loop over all distinct pairs i,j for (int j = i+1; j < N; j++){double rij[3]; // position of i relative to jdouble rSqd = 0;for (int k = 0; k < 3; k++){rij[k] = r[i][k] - r[j][k];rSqd += rij[k] * rij[k];}double f = 24 * (2 * pow(rSqd, -7) - pow(rSqd, -4));for (int k = 0; k < 3; k++){a[i][k] += rij[k] * f;a[j][k] -= rij[k] * f;}}}void velocityVerlet(double dt){computeAccelerations();for (int i = 0; i < N; i++)for (int k = 0; k < 3; k++){r[i][k] += v[i][k] * dt + 0.5 * a[i][k] * dt * dt;v[i][k] += 0.5 * a[i][k] * dt;}computeAccelerations();for (int i = 0; i < N; i++)for (int k = 0; k < 3; k++)v[i][k] += 0.5 * a[i][k] * dt;}double instantaneousTemperature(){double sum = 0;for (int i = 0; i < N; i++)for (int k = 0; k < 3; k++)sum += v[i][k] * v[i][k];return sum / (3 * (N - 1));}int main(){initialize();double dt = 0.01;ofstream FileTemp("Temp.dat");for (int i = 0; i < 5000; i++){cout<<i<<'\n';velocityVerlet(dt);FileTemp<<fixed \<<setw(15)<<setprecision(5)<<i \<<setw(20)<<setprecision(10)<<instantaneousTemperature()<< '\n';}FileTemp.close();}2)面心立方结构/*#include <cmath>#include <cstdlib>#include <fstream>#include <iostream>#include <string>*/#include <iostream>#include <string.h>#include <math.h>#include <iomanip>#include <stdlib.h>#include <fstream>#include <cstdio>using namespace std;// simulation parametersint N = 64; // number of particlesdouble rho = 1.0; // density (number per unit volume)double T = 1.5; // temperature// function declarationsvoid initialize(); // allocates memory, calls following 2 functions void initPositions(); // places particles on an fcc lattice，面心立方结构void initVelocities(); // initial Maxwell-Boltzmann velocity distribution void rescaleVelocities(); // adjust the instanteous temperature to Tdouble gasdev(); // Gaussian distributed random numbersdouble **r; // positionsdouble **v; // velocitiesdouble **a; // accelerationsvoid initialize(){r = new double* [N];v = new double* [N];a = new double* [N];for (int i = 0; i < N; i++){r[i] = new double [3];v[i] = new double [3];a[i] = new double [3];}initPositions();initVelocities();}double L; // linear size of cubical volumevoid initPositions(){// compute side of cube from number of particles and number density L = pow(N / rho, 1.0/3);// find M large enough to fit N atoms on an fcc latticeint M = 1;while (4 * M * M * M < N)++M;double a = L / M; // lattice constant of conventional cell// 4 atomic positions in fcc unit celldouble xCell[4] = {0.25, 0.75, 0.75, 0.25};double yCell[4] = {0.25, 0.75, 0.25, 0.75};double zCell[4] = {0.25, 0.25, 0.75, 0.75};int n = 0; // atoms placed so farfor (int x = 0; x < M; x++)for (int y = 0; y < M; y++)for (int z = 0; z < M; z++)for (int k = 0; k < 4; k++)if (n < N){r[n][0] = (x + xCell[k]) * a;r[n][1] = (y + yCell[k]) * a;r[n][2] = (z + zCell[k]) * a;++n;}}double gasdev (){static bool available = false;static double gset;double fac, rsq, v1, v2;if (!available){do{v1 = 2.0 * rand() / double(RAND_MAX) - 1.0;v2 = 2.0 * rand() / double(RAND_MAX) - 1.0;rsq = v1 * v1 + v2 * v2;}while (rsq >= 1.0 || rsq == 0.0);fac = sqrt(-2.0 * log(rsq) / rsq);gset = v1 * fac;available = true;return v2*fac;}else{available = false;return gset;}}void initVelocities(){// Gaussian with unit variancefor (int n = 0; n < N; n++)for (int i = 0; i < 3; i++)v[n][i] = gasdev();// Adjust velocities so center-of-mass velocity is zerodouble vCM[3] = {0, 0, 0};for (int n = 0; n < N; n++)for (int i = 0; i < 3; i++)vCM[i] += v[n][i];for (int i = 0; i < 3; i++)vCM[i] /= N;for (int n = 0; n < N; n++)for (int i = 0; i < 3; i++)v[n][i] -= vCM[i];// Rescale velocities to get the desired instantaneous temperature rescaleVelocities();}void rescaleVelocities(){double vSqdSum = 0;for (int n = 0; n < N; n++)for (int i = 0; i < 3; i++)vSqdSum += v[n][i] * v[n][i];double lambda = sqrt( 3 * (N-1) * T / vSqdSum );for (int n = 0; n < N; n++)for (int i = 0; i < 3; i++)v[n][i] *= lambda;}void computeAccelerations(){for (int i = 0; i < N; i++) // set all accelerations to zero for (int k = 0; k < 3; k++)a[i][k] = 0;for (int i = 0; i < N-1; i++) // loop over all distinct pairs i,j for (int j = i+1; j < N; j++){double rij[3]; // position of i relative to jdouble rSqd = 0;for (int k = 0; k < 3; k++){rij[k] = r[i][k] - r[j][k];// closest image conventionif (abs(rij[k]) > 0.5 * L){if (rij[k] > 0)rij[k] -= L;elserij[k] += L;}rSqd += rij[k] * rij[k];}double f = 24 * (2 * pow(rSqd, -7) - pow(rSqd, -4));for (int k = 0; k < 3; k++){a[i][k] += rij[k] * f;a[j][k] -= rij[k] * f;}}}void velocityVerlet(double dt){computeAccelerations();for (int i = 0; i < N; i++)for (int k = 0; k < 3; k++){r[i][k] += v[i][k] * dt + 0.5 * a[i][k] * dt * dt;// use periodic boundary conditionsif (r[i][k] < 0)r[i][k] += L;if (r[i][k] >= L)r[i][k] -= L;v[i][k] += 0.5 * a[i][k] * dt;}computeAccelerations();for (int i = 0; i < N; i++)for (int k = 0; k < 3; k++)v[i][k] += 0.5 * a[i][k] * dt;}double instantaneousTemperature(){double sum = 0;for (int i = 0; i < N; i++)for (int k = 0; k < 3; k++)sum += v[i][k] * v[i][k];return sum / (3 * (N - 1));}int main(){initialize();double dt = 0.01;ofstream FileTemp("Temp2.dat");for (int i = 0; i < 4000; i++){cout<<i<<'\n';velocityVerlet(dt);FileTemp<<fixed \<<setw(15)<<setprecision(5)<<i \<<setw(20)<<setprecision(10)<<instantaneousTemperature()<< '\n';if (i % 200 == 0)rescaleVelocities();}FileTemp.close();}3.计算结果及具体分析讨论1.两种结构最终结果均收敛到某个值附近，面心立方结构收敛速度更快，并且收敛以后的方差明显更小。