熔体和流体是自然界中物质迁移的普遍载体。我们开发了一个开源软件包，用于分析此类体系的从头算分子动力学模拟。我们计算了其结构特性（成键、团簇化、化学组分）、输运特性（扩散、黏度）以及热力学特性（振动谱）。
本文内容
摘要
摘要
我们开发了一个基于 Python 的开源软件包，用于分析流体的从头算分子动力学模拟结果。该软件包最适合应用于天然体系，如硅酸盐和氧化物熔体、水基流体以及各种超临界流体。该软件包由一系列 Python 脚本组成，包含两个主要的库，分别用于处理文件格式和晶体学相关计算。所有脚本均在命令行中运行。我们提出了一种简化的格式，用于存储模拟中原子轨迹及相关热力学信息，该格式保存为 UMD 文件，即通用分子动力学（Universal Molecular Dynamics）文件。UMD 软件包能够计算一系列结构、输运和热力学性质。从径向分布函数出发，该软件可定义键长、构建原子间连接矩阵，并最终确定化学组分。通过确定化学物种的寿命，可进行完整的统计分析。随后，专用脚本可计算原子以及化学物种的均方位移。对原子速度进行自相关分析，可得到扩散系数和振动谱。将相同的分析应用于应力，则可获得黏度。该软件包可通过 GitHub 网站以及 ERC IMPACT 项目专属页面免费获取。
引言
流体和熔体是自然环境中活跃的化学与物理传输载体。较高的原子扩散速率有利于化学交换与反应，低黏度结合变化的浮力有利于大规模物质传输，而晶体与熔体之间的密度关系则有利于行星内部的分层。由于缺乏周期性晶格、达到熔融状态通常需要高温，以及淬火困难，使得密度、扩散和黏度等明显性质的实验测定极具挑战性。这些困难使得替代性的计算方法成为研究此类材料的有力且有用的工具。
随着计算能力的发展以及超级计算机的普及，目前研究非晶态原子体系动力学状态主要有两种重要的数值原子模拟技术：蒙特卡洛方法1和分子动力学（MD）1,2。在蒙特卡洛模拟中，构型空间通过随机采样获得；若所有采样观测彼此独立，则蒙特卡洛方法在并行化过程中表现出线性扩展性。结果的准确性取决于随机数生成器的质量以及采样的代表性。当采样相互独立时，蒙特卡洛方法在并行化中具有线性扩展性。在分子动力学（MD）中，构型空间通过依赖时间的原子轨迹进行采样。从给定的初始构型出发，通过积分牛顿运动方程来计算原子轨迹。原子间作用力可通过模型原子间势函数（在经典MD中）或基于第一性原理的方法（在从头算或第一性原理MD中）计算得到。结果的准确性取决于轨迹的长度及其避免陷入局部能量极小值的能力。
分子动力学模拟包含大量信息，均与系统的动力学行为相关。热力学平均性质（如内能、温度和压力）的计算较为常规，可直接从模拟的输出文件中提取并进行平均；而与原子运动及其相互关系直接相关的物理量，则需要在提取原子坐标和速度之后进一步计算。
因此，人们投入了大量精力用于结果的可视化，目前在不同平台上已有多种可视化软件包可供使用，包括开源和非开源的软件[Ovito3, VMD4, Vesta5, Travis6等]。所有这些可视化工具均能高效处理原子间距离，从而可高效计算径向分布函数和扩散系数。许多开展大规模分子动力学模拟的研究团队拥有专有软件，用于分析模拟产生的其他各类性质，这些软件有时以共享软件形式或以其他受限方式向学术界开放，有时则仅限于特定软件包的范围和用途。一些软件包中已开发并实现了复杂的算法，用于提取原子间成键、几何构型以及热力学等方面的信息3,4,5,6,7等。
本文提出UMD 软件包——一个用 Python 编写的开源软件包，用于分析分子动力学模拟的输出结果。UMD 软件包能够计算多种结构、动力学和热力学性质（图 1）。该软件包可通过 GitHub 网站（https://github.com/rcaracas/UMD_package）以及欧洲研究理事会（ERC）IMPACT 项目专属页面（http://moonimpact.eu/umd-package/）免费获取。
为了使其具有通用性并更易于处理，我们的方法是首先从实际分子动力学模拟的输出文件中提取所有与热力学状态和原子轨迹相关的信息。这些信息被存储在一个专用文件中，其格式独立于运行模拟的原始分子动力学软件包。我们将此类文件命名为“umd”文件，即通用分子动力学（Universal Molecular Dynamics）的缩写。通过这种方式，我们的 UMD 软件包可被任何使用任意软件的从头算（ab initio）研究小组轻松使用，且仅需最少的适配工作。使用本软件包的唯一要求是，将特定分子动力学软件的输出解析为 umd 文件格式，如果尚无相应解析器，则需自行编写。目前，我们已为 VASP8 和 QBox9 软件包提供了相应的解析器。
图1：UMD库的流程图。
物理性质用蓝色表示，主要的Python脚本及其选项用红色表示。请点击此处查看该图的放大版本。
umd 文件是 ASCII 文件，通常的扩展名为“umd.dat”，但并非强制要求。所有分析组件均可读取 umd 格式的 ASCII 文件，而不论其实际文件扩展名如何。然而，某些为在多个模拟上快速进行大规模统计而设计的自动化脚本，会专门查找扩展名为 umd.dat 的文件。每个物理量在文件中占一行，每行以一个关键字开头。因此，该格式具有高度可扩展性，允许向 umd 文件中添加新的物理量，同时在不同版本间保持良好的可读性。下文讨论中所使用的、针对 4.6 GPa 和 3000 K 条件下辉石岩模拟的 umd 文件的前 30 行内容如 图 2 所示。
图 2：描述在 4.6 GPa 和 3000 K 条件下液态榴辉岩模拟的 umd 文件开头部分。
文件头之后是每个快照的描述。每种物理性质单独占一行，包含物理性质的名称、数值和单位，各项之间以空格分隔。请点击此处查看此图的放大版本。
所有 umd 文件均包含一个描述模拟体系内容的文件头：原子数、电子数和原子类型，以及每个原子的详细信息，例如其类型、化学符号、价电子数和质量。一个空行标记了文件头的结束，并将其与 umd 文件的主体部分分隔开。
随后详细描述模拟的每一步。首先，逐行列出瞬时热力学参数，每个参数包括：(i) 参数名称，如能量、应力、等效静水压、密度、体积、晶格参数等；(ii) 参数值；以及 (iii) 单位。接下来是一个描述原子的表格。表头行列出各项物理量，如笛卡尔坐标、速度、电荷等，以及它们的单位。随后每一行详细描述一个原子。每三个数据为一组，对应 x, y, z 三个轴，依次为：约化坐标、折叠到模拟原胞内的笛卡尔坐标、真实的笛卡尔坐标（能够正确反映原子在模拟过程中可能跨越多个晶胞的情况）、原子速度和原子受力。最后两个条目为标量：电荷和磁矩。
两个主要的库确保了整个软件包的正常运行。umd_process.py 库负责处理 umd 文件，例如读取和输出。crystallography.py 库则处理与实际原子结构相关的所有信息。crystallography.py 库的核心理念是将晶格视为一个矢量空间。晶胞参数及其取向构成了该空间的基向量。该“空间”具有一系列标量属性（如比容、密度、温度和原子数密度）、热力学性质（如内能、压强、热容等）以及一系列张量性质（如应力和弹性）。原子分布于该空间之中。“Lattice”类定义了这一整体结构，并包含若干简单的计算功能，例如计算比容、密度、由正格子求倒格子等。“Atoms”类用于定义原子，其特征包括一系列标量属性（如名称、符号、质量、电子数等）和一系列矢量属性（如空间位置——可相对于 Lattice 类中定义的矢量基，也可相对于通用笛卡尔坐标系——以及速度、受力等）。除这两个类之外，crystallography.py 库还包含一系列用于执行多种测试和计算的函数，例如计算原子间距或晶胞倍增。元素周期表也以字典形式内置于库中。
umd 软件包的各个组件会生成多个输出文件。通常情况下，这些文件均为 ASCII 文件，所有条目均以制表符分隔，并尽可能做到自明性强。例如，文件中始终会明确标注物理量及其单位。umd.dat 文件完全符合此规范。
方案
1. 分子动力学模拟结果的分析
注意：该软件包可通过 GitHub 网站（https://github.com/rcaracas/UMD_package）以及 ERC IMPACT 项目专用页面（http://moonimpact.eu/umd-package/）获取，为开放获取软件包。
- 使用该软件包中的一个或多个专用 Python 脚本提取每组特定的物理性质。在命令行中运行所有脚本；这些脚本均采用一系列参数标志，且尽可能在不同脚本之间保持标志的一致性。各标志的含义及其默认值均在表 1中汇总。
| 标志 | 含义 | 使用该标志的脚本 | 默认值 |
| -h | 简要帮助信息 | 所有脚本 | |
| -f | UMD 文件名 | 所有脚本 | |
| -i | 需舍弃的热化步数 | 所有脚本 | 0 |
| -i | 包含原子间键的输入文件 | speciation | bonds.input |
| -s | 频率采样间隔 | msd, speciation | 1（每一步均被考虑） |
| -a | 原子或阴离子列表 | speciation | |
| -c | 阳离子列表 | speciation | |
| -l | 键长 | speciation | 2 |
| -t | 温度 | vibrations, rheology | |
| -v | 均方位移分析中轨迹采样窗口宽度的离散化分段数 | msd | 20 |
| -z | 均方位移分析中轨迹采样窗口起始位置的离散化分段数 | msd | 20 |
表1：UMD软件包中最常用的参数标志及其最常见的含义。
- 首先，将使用第一性原理计算程序（如 VASP）完成的分子动力学模拟结果进行转换8 或 QBox9，转换为 UMD 文件。
- 如果在 VASP 中完成了分子动力学模拟，则在命令行中输入：
VaspParser.py -f -i <初始步骤>
其中，–f 参数定义 VASP OUTCAR 文件的名称，–i 参数定义热化长度。
注意：初始步骤由 –i 定义，可用于舍弃模拟的最初阶段，该阶段对应于热化过程。在典型的分子动力学运行中，计算的初始部分代表热化，即系统中所有原子达到温度的类高斯分布，以及整个系统在温度、压力、能量等参数上围绕平衡值波动所需的时间。在分析流体的统计性质时，不应考虑模拟中的这一热化阶段。
- 如果在 VASP 中完成了分子动力学模拟，则在命令行中输入：
- 转化。umd 文件转换为.xyz 文件，以便在其他多种软件（如 VMD）中进行可视化4 或 Vesta5在命令行中输入：
umd2xyz.py -f -i <初始步骤> - s <采样频率>
其中 -f 定义了 . 的文件名umd 文件中，–i 定义了要舍弃的热化阶段时间，–s 定义了存储在 . 中的轨迹的采样频率。umd 文件。默认值为 –i 0 –s 1，即考虑模拟的所有步骤，不丢弃任何步骤。 - 逆转 umd 使用 umd2poscar.py 脚本将文件转换为 VASP 类型的 POSCAR 文件；可通过预设频率选取模拟的快照。在命令行中输入：
umd2poscar.py -f -i <初始步骤> - l <最后一步> - s <采样频率>
其中 -l 表示要转换为 POSCAR 文件的最后一步。默认值为 -i 0 -l 10000000 -s 1。该 -l 值足够大，可覆盖典型的完整轨迹。
2. 进行结构分析
- 运行 gofrs_umd.py 脚本，计算所有原子类型 A 和 B 之间的原子对分布函数（PDF）gᴀʙ(r)（图3）。输出结果写入一个以制表符分隔的 ASCII 文件中，文件扩展名为 gofrs.dat。在命令行中输入以下命令：
gofrs_umd.py -f -s < Sampling_Frequency > -d -i
注：默认参数为采样频率（轨迹采样的频率）= 1 步；离散化间隔（用于绘制 g(r)）= 0.01 Å；初始步数（轨迹开始时被舍弃的步数）= 0。径向对分布函数 gᴀʙ(r) 表示以类型 A 原子为中心、半径为 r、厚度为 dr 的球壳内，距离为 d_ᴀʙ 处类型 B 原子的平均数目（图3）：
其中 ρ 为原子密度，NA 和 NB 分别为类型 A 和 B 的原子数目，δ(r−rᴀʙ) 为 delta 函数，当原子 A 与 B 之间的距离位于 r 与 r+dr 之间时，其值为 1。gᴀʙ(r) 第一个极大值对应的横坐标给出了类型 A 与 B 原子之间最可能的键长，该值最接近我们所能确定的平均键长。第一个极小值界定了第一配位层的范围。因此，对 PDF 从原点积分至第一个极小值处可得到平均配位数。对所有原子类型对 A 和 B 的 gᴀʙ(r) 进行傅里叶变换并求和，可得到流体的衍射图样，这与使用衍射仪进行实验测得的结果一致。然而实际上，由于 gᴀʙ(r) 中常常缺少高阶配位层的信息，因此无法完整获得衍射图样。
图3：对分布函数的确定。
（a）对于某一物种的每个原子（例如红色），统计其周围配位物种（例如灰色和/或红色）的原子数量随距离的变化。（b）对每个快照得到的距离分布图（在此阶段仅为一系列δ函数）在所有原子和所有快照上进行平均，并以理想气体分布加权，从而生成（c）连续的对分布函数。g(r) 的第一个极小值对应第一配位层的半径，该值将在后续的物种分析中使用。请点击此处查看该图的放大版本。
- 提取平均原子间键距作为第一配位层的半径。为此，确定 gᴀʙ(r) 函数的第一个极大值位置：在电子表格软件中绘制 gofrs.dat 文件，并查找每对原子的极大值和极小值。
- 使用电子表格软件，将 PDF 函数 gᴀʙ(r) 的第一个极小值位置确定为第一配位层的半径。这是整个流体结构分析的基础；PDF 可提供流体中原子的平均成键状态。
- 提取第一个极小值的距离（即横坐标），并将其写入一个单独的文件，例如命名为 bonds.input。或者，运行 UMD 软件包中的某个 analyze_gofr 脚本，以识别 gᴀʙ(r) 函数的极大值和极小值。在命令行中输入：
analyze_gofr_semi_automatic.py - 在程序打开的图形中，点击 gᴀʙ(r) 函数曲线上的极大值和极小值位置。该脚本将自动扫描当前文件夹，识别所有 gofrs.dat 文件，并对每个文件执行分析。每次脚本需要合理的初始估计值时，请在窗口中再次点击相应的极大值和极小值位置。
- 打开并查看自动生成的 bonds.input 文件，该文件包含原子间的键距数据。
3. 进行形态分析
- 利用图论中的连接性概念，计算原子间化学键的拓扑结构：原子作为节点，原子间的化学键作为路径。speciation_umd.py 脚本需要在 bonds.input 文件中定义的原子间键长。
注意：在每个时间步都会构建连接性矩阵：当两个原子之间的距离小于其对应的第一配位层半径时，即被视为成键，即相互连接。通过将原子视为图中的节点，并根据该几何准则定义其连接关系，构建出多种原子网络。这些网络即为原子物种，其集合定义了该特定流体中的原子形态分布（图 4）。
图4：原子簇的识别。
配位多面体通过原子间距离来定义。所有距离小于指定半径的原子均被视为成键。此处阈值对应于第一配位层（浅红色圆圈），如图1中所定义。聚合结构及相应的化学物种由成键原子形成的网络确定。注意中心的Red1Grey2原子簇与其他原子隔离，其余原子则形成无限延伸的聚合物。请点击此处查看该图的放大版本。
- 运行物种分析脚本以获得连接性矩阵，并获取配位多面体或聚合信息。在命令行中输入：
speciation_umd.py -f -s -i -l -c -a -m -r
其中 -i 参数指定包含原子间键长的文件，例如上一步生成的文件。或者，可使用 -l 参数为所有键指定单一长度来运行脚本。
注意：-c 参数用于指定中心原子，-a 参数用于指定配体。中心原子和配体均可为多种类型；若存在多种类型，必须用逗号分隔。-m 参数定义某一物种在分析中被计入所需的最短存在时间。默认情况下该时间为零，即所有出现情况均计入最终分析。- 使用 –r 0 参数运行 speciation_umd.py 脚本，该参数在第一层级上采样连接性图以识别配位多面体。例如，一个标记为中心原子的阳离子可能被一个或多个阴离子包围（图4）。物种分析脚本会识别出每一个配位多面体。所有配位多面体的加权平均值即为配位数，该数值与PDF积分所得结果一致。在命令行中输入：
speciation_umd.py -f -i -c -a -r 0
注意：流体中的平均配位数通常为分数。这种分数特性源于配位数的平均性质。基于物种分析的定义能够更直观且信息丰富地表征流体结构，其中不同物种（即不同配位方式）的相对比例被量化。 - 使用 –r 1 参数运行 speciation_umd.py 脚本，该参数在所有深度层级上采样连接性图以获得聚合信息。原子图构成的网络具有一定深度，因为原子通过键进一步连接到更远的其他原子（例如交替排列的阳离子和阴离子序列）（图4）。
- 使用 –r 0 参数运行 speciation_umd.py 脚本，该参数在第一层级上采样连接性图以识别配位多面体。例如，一个标记为中心原子的阳离子可能被一个或多个阴离子包围（图4）。物种分析脚本会识别出每一个配位多面体。所有配位多面体的加权平均值即为配位数，该数值与PDF积分所得结果一致。在命令行中输入：
- 依次打开两个输出文件 .popul.dat 和 .stat.dat；这两个文件是物种分析脚本的输出结果。每个团簇占据一行，记录其化学式、形成时间、消失时间、寿命，以及构成该团簇的原子列表矩阵。根据 .popul.dat 文件中的数据绘制模拟中发现的所有原子团簇的寿命分布图（图5）。
- 根据 .stat.dat 文件中的数据绘制各物种的种群丰度分析图。该分析包括绝对丰度和相对丰度：对于 -r 0 情况，对应于配位多面体的实际统计结果；而对于 -r 1 的聚合情况，则需谨慎处理，可能需要对原子的相对数量进行归一化。丰度值对应于寿命的积分。.stat.dat 文件还列出了每个团簇的尺寸，即构成该团簇的原子数量。
4. 计算扩散系数
- 提取原子均方位移（MSD）随时间的变化关系，以获得自扩散系数。MSD 的标准公式为：
其中前置因子为重正化项。利用均方位移（MSD）工具，可通过多种方法分析流体的动力学特性。
注意：T 为模拟的总时间 Nα 是某类原子的数量 α初始时间 t0 是任意的，跨越模拟的前半段。N启动 是初始次数。 τ 是计算均方位移（MSD）所采用的时间区间宽度；其最大值为模拟总时长的一半。在典型的MSD实现中，每个时间窗口紧接在前一个窗口结束处开始。但采用更稀疏的采样方式可在不改变MSD斜率的前提下加快计算速度。为此，第i个窗口从时间点开始 t0(i），但第 (i+1) 个窗口从时间 t0(i) + τ + v，其中该值为 v 由用户自定义。类似地，窗口的宽度也以用户定义的离散步长逐步增加，如下所示： τ(i) = τ(i-1) + z。其数值为 z (“水平步长”和 v (“垂直步进”为正值或零；两者的默认值均为20。 - 使用该序列计算均方位移（MSD） msd_umd 脚本。其输出内容将打印在 .msd.dat 文件中，每种原子类型、单个原子或原子簇的均方位移（MSD）作为时间的函数，以列为单位输出。
- 计算每种原子类型的平均均方位移（MSD）。首先针对每个原子计算其MSD，然后对每种原子类型进行平均。输出文件中每种原子类型对应一列。在命令行中输入：
msd_umd.py -f -z <水平跳跃> - v <垂直跳跃> - b <弹道学> - 计算每个原子的均方位移（MSD）。对模拟中每个原子分别计算其MSD，然后按原子类型进行平均。输出文件包含模拟中每个原子对应的一列，以及每种原子类型对应的一列。该功能可用于识别在两种不同环境中扩散的原子，例如液体和气体，或两种不同的液体。在命令行中输入：
msd_all_umd.py -f -z <水平跳跃> - v <垂直跳跃> - b <弹道学> - 计算化学物种的均方位移（MSD）。利用通过物种分析脚本识别并输出到 . 文件中的簇群群体进行计算。popul.dat 文件。对每个独立的簇计算均方位移（MSD）。输出文件中每一列对应一个簇。为避免考虑大尺度聚合物，需对簇的大小设置限制；其默认值为20个原子。在命令行中输入：
msd_cluster_umd.py -f - p - s <采样频率> - b <弹道学> - c <最大簇大小>
注：默认值为 –b 100 –s 1 –c 20。
- 计算每种原子类型的平均均方位移（MSD）。首先针对每个原子计算其MSD，然后对每种原子类型进行平均。输出文件中每种原子类型对应一列。在命令行中输入：
- 使用基于电子表格的软件绘制均方位移（MSD）曲线图6）。在均方位移（MSD）与时间的双对数图中，识别斜率变化。将代表初始阶段（通常较短）的第一部分分离出来，该部分表示 弹道的 即碰撞后原子速度的守恒。第二部分较长，代表了 扩散的 即碰撞后原子速度的散射行为。
- 根据均方位移（MSD）的斜率计算扩散系数，公式如下：
其中，Z 为自由度数目（平面扩散时 Z = 2，空间扩散时 Z = 3），t 为时间步长。
5. 时间关联函数
- 使用以下通用公式计算时间关联函数，作为衡量系统惯性的指标：
A 可以是多种随时间变化的变量，例如原子位置、原子速度、应力、极化等，每种变量通过格林-久保关系12,13 可得到不同的物理性质，有时还需进一步变换。 - 分析原子速度，以获得液体的振动谱以及原子自扩散系数的另一种表达形式。
- 运行 vibr_spectrum_umd.py 脚本，计算每种原子类型的原子速度-速度自相关（VAC）函数，并对其进行快速傅里叶变换。在命令行中输入：
vibr_spectrum_umd.py -f -t
其中 –t 是用户必须定义的温度。该脚本将输出两个文件：.vels.scf.dat 文件包含每种原子类型的 VAC 函数，.vibr.dat 文件包含按各原子种类分解的振动谱及总谱。 - 打开并读取 vels.scf.dat 文件。使用类似电子表格的软件绘制 vels.scf.dat 文件中的 VAC 函数。
- 保留傅里叶变换后 VAC 的实部。这将作为频率的函数给出振动谱：
其中 m 为原子质量。 - 使用类似电子表格的软件绘制 vibr.dat 文件中的振动谱（图7）。识别在 ω=0 处的有限值，该值对应流体的扩散特性，以及在有限频率处谱图上的多个峰。确定每种原子类型对振动谱的贡献。
注：按原子类型分解的结果显示，不同原子在 ω=0 处的贡献不同，对应其各自的扩散系数。与相应的固体相比，该谱图的整体形状更平滑，特征更少。 - 在终端读取对振动谱的积分，该积分可得到每种原子种类的扩散系数。
注：可通过振动谱积分获得热力学性质，但由于两个近似，结果应谨慎使用：积分在准谐近似下成立，而该近似在高温下不一定适用；此外，对应扩散的类气体部分谱需被舍弃。因此，积分应仅在类晶格部分的谱上进行。但这种分离通常需要多个额外的后处理步骤和计算14，当前 UMD 软件包未涵盖这些内容。
- 运行 vibr_spectrum_umd.py 脚本，计算每种原子类型的原子速度-速度自相关（VAC）函数，并对其进行快速傅里叶变换。在命令行中输入：
- 运行 viscosity_umd.py 脚本，分析应力张量各分量的自相关性，以估算熔体的黏度。在命令行中输入：
viscosity_umd.py -f -i -s -o -l
注：此功能尚处于探索阶段，任何结果均需谨慎对待。首先，应充分检查黏度随模拟时长的收敛性。- 根据应力张量的自相关性推导流体的黏度15：
其中 V 和 T 分别为体积和温度，κB 为玻尔兹曼常数，σij 为应力张量的 ij 非对角分量，以笛卡尔坐标表示。 - 采用更合适的拟合方法，以获得更稳健的黏度估计值15,16，并避免因模拟的有限尺寸和有限持续时间导致的应力张量自相关函数噪声。对于应力张量的自相关函数，使用以下函数形式15,16，其可获得良好结果：
其中 A、B、τ1、τ2 和 ω 为拟合参数。积分后，黏度的表达式为：
- 根据应力张量的自相关性推导流体的黏度15：
6. 模拟产生的热力学参数。
- 运行 averages.py 以从 umd 文件中提取压力、温度、密度和内能的平均值及其离散程度（以标准偏差表示）。在命令行中输入：
averages.py -f -s
其中 –s 0 为默认值。 - 使用分块法计算平均值的统计误差。
注：该方法有多种变体。根据 Allen 和 Tildesley2 的研究，通常做法是对时间块序列进行平均，时间块长度逐渐增加，并相对于算术平均值估算标准偏差17。当采样无关联时，在块数足够多且块长度足够长的极限情况下可达到收敛。然而，收敛的实际阈值通常需要手动选择。- 采用二分法18：从初始数据样本开始，在每一步 κ 中，通过将前一步 κ−1 中每两个相邻样本取平均，使样本数量减半：
- 运行 fullaverages.py 脚本以执行完整的统计分析，包括均值误差的估计。在命令行中输入：
fullaverages.py -s -u
注：该脚本已实现自动化，能够搜索当前目录下所有 .umd.dat 文件并对它们逐一进行分析。默认参数为 –s 0 –u 0。当 -u 0 时输出最简，当 -u 1 时输出完整，并打印多种可选单位。该脚本需要图形支持，因为它会生成图形图像以检查收敛性，从而评估均值的误差。
- 采用二分法18：从初始数据样本开始，在每一步 κ 中，通过将前一步 κ−1 中每两个相邻样本取平均，使样本数量减半：
结果
辉石岩（Pyrolite）是一种多组分硅酸盐熔体模型（0.5Na2O 2CaO 1.5Al2O3 4FeO 30MgO 24SiO2），其成分最接近于整体硅酸盐地球——即除铁核以外整个地球的地球化学平均组成19。早期地球以一系列大规模熔融事件为主导20，其中最后一次可能在原月盘凝聚之后席卷了整个行星21。辉石岩最能近似代表此类行星尺度岩浆洋的成分。因此，我们利用VASP程序中的从头算分子动力学模拟，系统研究了辉石岩熔体在3,000‒5,000 K温度范围和0‒150 GPa压力范围内的物理性质。这些热力学条件完整表征了地球历史上最极端的岩浆洋状态。本研究是成功应用UMD软件包对熔体进行深入分析的优秀范例22。我们计算了键长的分布与平均值，追踪了阳离子-氧配位数的变化，并将结果与以往关于不同成分非晶态硅酸盐的实验和计算研究进行了比较。深入分析帮助我们将标准配位数分解为其基本组分，揭示了熔体中存在非常规配位多面体，并提取了所有配位多面体的寿命。此外，分析还强调了模拟中采样的重要性，包括轨迹长度以及所建模体系中原子数量的影响。在后处理方面，UMD分析本身独立于这些因素，但在解释UMD软件包提供的结果时仍需考虑这些因素。本文展示了如何利用UMD软件包提取熔体若干特征性质的几个示例，具体应用于熔融态辉石岩。
由 gofrs_umd.py 脚本获得的 Si-O 对分布函数显示，在 T = 3000 K 和 P = 4.6 GPa 条件下，第一配位壳层的半径（即 g(r) 函数的第一个极小值）约为 2.5 埃。而 g(r) 函数的峰值位于 1.635 Å，这最接近键长的最佳近似值。长尾现象是由温度引起的。以该值作为 Si-O 键长界限，物种分析表明，SiO4 单元在熔体中占主导地位，其存在时间可长达数皮秒（图 5）。熔体中有相当一部分表现出部分聚合特征，这体现在存在如 Si2O7 类型的二聚体以及 Si3Ox 类型的三聚体单元。它们相应的寿命在皮秒量级。更高阶的聚合物寿命则显著更短。
图 5：Si-O 化学物种的寿命。
该物种鉴定于 4.6 GPa 和 3000 K 条件下的多组分熔体中。标签标示了 SiO3、SiO4 和 SiO5 单体以及各种 SixOy 多聚体。请点击此处查看此图的放大版本。
由上述 –z 和 –v 标志定义的垂直与水平步长的不同取值，会产生均方位移（MSD）的多种采样结果（图6）。即使 z 和 v 取较大值，也足以确定斜率，从而确定不同原子的扩散系数。在后处理中采用较大的 z 和 v 值时，时间上的节省效果显著。MSD 为模拟质量提供了非常强的验证标准。如果 MSD 的扩散部分不够长，则表明模拟时间过短，在统计意义上未能达到流体状态。MSD 扩散部分的最低要求在很大程度上取决于具体体系。可以要求熔体结构中所有原子至少更换一次位置，才能将其视为流体10。行星科学中的一个极佳实例是接近甚至低于其液相线的高压复杂硅酸盐熔体11。作为主要网络形成阳离子的 Si 原子，在超过几十皮秒后才发生位点交换。短于该阈值的模拟将显著欠采样可能的构型空间。然而，由于配位阴离子（即 O 原子）的运动速度比中心 Si 原子更快，它们可以在一定程度上补偿 Si 原子缓慢的迁移性。因此，整个体系的实际构型空间采样可能优于仅根据 Si 原子位移所推测的结果。
图6：均方位移（MSD）。
本图展示了多组分硅酸盐熔体中几种原子类型的均方位移。采用不同的水平和垂直步长（z 和 v）进行采样，结果具有一致性。实心圆点：-z 50 –v 50；空心圆点：-z 250 –v 500。请点击此处查看此图的放大版本。
最后，原子速度自相关函数（VAC）可得到熔体的振动光谱。图7展示了与上述相同的压强和温度条件下的光谱。我们分别表示了Mg、Si和O原子的贡献，以及总和值。在零频率处，光谱存在一个有限值，这对应于熔体的扩散特性。从振动光谱中提取热力学性质时，需要去除零频处类似气体的扩散特征，同时还需恰当地考虑其在较高频率下的衰减行为。
图7：辉石熔体的振动光谱。
原子速度-速度自相关函数的傅里叶变换的实部给出振动谱. 此处针对多组分硅酸盐熔体计算其谱图。流体在零频率下具有非零的类气体扩散特性。 请点击此处以查看此图的放大版本。
讨论
UMD 软件包的设计更适用于从头算模拟，这类模拟通常仅包含数万至数十万张快照，每个晶胞内含有数百个原子。只要后处理运行所用机器具备足够的活动内存资源，更大规模的模拟也可实现。该代码的独特之处在于其能够计算多种性质，并采用开源许可证。
用于在模拟过程中保持粒子数量不变的系综时，umd.dat 文件是适用的。UMD 程序包能够读取来自模拟盒子形状和体积发生变化的计算所产生的文件。这些文件涵盖了最常见的计算类型，例如 NVT 和 NPT 系综，其中粒子数 N、温度 T、体积 V 和/或压力 P 保持恒定。
由于时间关系，径向分布函数以及所有用于估算原子间距离的脚本（例如物种分布脚本）仅适用于正交晶胞，即立方、四方和正交晶系，其中晶轴之间的夹角为90°。
2.0 版本的主要开发方向是取消距离计算中的正交性限制，并为物种分析脚本增加更多功能：分析单个化学键、分析原子间键角，以及实现第二配位层的分析。在外部合作的支持下，我们正在将代码移植到 GPU 上，以实现对更大体系的快速分析。
披露
作者无任何利益冲突需要披露。
致谢
本工作由欧洲研究理事会（ERC）根据欧洲联盟“地平线2020”研究与创新计划（资助协议号681818 IMPACT，授予RC）、深部碳观测计划的极端物理与化学方向，以及挪威研究理事会通过其卓越研究中心资助计划（项目编号223272）提供支持。我们感谢通过eDARI计算资助系列stl2816获得GENCI超级计算机的使用权限，通过PRACE RA4947项目获得Irene AMD超级计算机的使用权限，以及通过UNINETT Sigma2 NN9697K项目获得Fram超级计算机的使用权限。FS由玛丽·斯克沃多夫斯卡-居里项目（资助协议ABISSE No.750901）提供支持。
材料
| 姓名 | 公司 | 目录编号 | 评论 |
|---|---|---|---|
| getopt 库 | 开源 | ||
| glob 库 | 开源 | ||
| matplotlib 库 | 开源 | ||
| numpy 库 | 开源 | ||
| os 库 | 开源 | ||
| Python 软件 | The Python Software Foundation | 版本 2 和 3 | 开源 |
| random 库 | 开源 | ||
| re 库 | 开源 | ||
| scipy 库 | 开源 | ||
| subprocess 库 | 开源 | ||
| sys 库 | 开源 |
参考文献
- Frenkel, D., Smit, B. Understanding Molecular Simulation. From Algorithms to Applications. , Elsevier. (2001).
- Allen, M. P., Tildesley, D. J., Allen, T. Computer Simulation of Liquids. , Oxford University Press. (1989).
- Zepeda-Ruiz, L. A., Stukowski, A., Oppelstrup, T., Bulatov, V. V. Probing the limits of metal plasticity with molecular-dynamics simulations. Nature Publishing Group. 550 (7677), 492-495 (2017).
- Humphrey, W., Dalke, A., Schulten, K. VMD: Visual molecular dynamics. Journal of Molecular Graphics & Modeling. 14 (1), 33-38 (1996).
- Momma, K., Izumi, F. VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data. Journal of Applied Crystallography. 44 (6), 1272-1276 (2011).
- Brehm, M., Kirchner, B. TRAVIS - A free Analyzer and Visualizer for Monte Carlo and Molecular Dynamics Trajectories. Journal of Chemical Information and Modeling. 51 (8), 2007-2023 (2011).
- Stixrude, L. Visualization-based analysis of structural and dynamical properties of simulated hydrous silicate melt. Physics and Chemistry of Minerals. 37 (2), 103-117 (2009).
- Kresse, G., Hafner, J. Ab initio Molecular-Dynamics for Liquid-Metals. Physical Review B. 47 (1), 558-561 (1993).
- Gygi, F. Architecture of Qbox: A scalable first-principles molecular dynamics code. IBM Journal of Research and Development. 52 (1-2), 137-144 (2008).
- Harvey, J. P., Asimow, P. D. Current limitations of molecular dynamic simulations as probes of thermo-physical behavior of silicate melts. American Mineralogist. 100 (8-9), 1866-1882 (2015).
- Caracas, R., Hirose, K., Nomura, R., Ballmer, M. D. Melt-crystal density crossover in a deep magma ocean. Earth and Planetary Science Letters. 516, 202-211 (2019).
- Green, M. S. Markoff Random Processes and the Statistical Mechanics of Time-Dependent Phenomena. II. Irreversible Processes in Fluids. The Journal of Chemical Physics. 22 (3), 398-413 (1954).
- Kubo, R. Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems. Journal of the Physical Society of Japan. 12 (6), 570-586 (1957).
- Lin, S. T., Blanco, M., Goddard, W. A. The two-phase model for calculating thermodynamic properties of liquids from molecular dynamics: Validation for the phase diagram of Lennard-Jones fluids. The Journal of Chemical Physics. 119 (22), 11792-11805 (2003).
- Meyer, E. R., Kress, J. D., Collins, L. A., Ticknor, C. Effect of correlation on viscosity and diffusion in molecular-dynamics simulations. Physical Review E. 90 (4), 1198-1212 (2014).
- Soubiran, F., Militzer, B., Driver, K. P., Zhang, S. Properties of hydrogen, helium, and silicon dioxide mixtures in giant planet interiors. Physics of Plasmas. 24 (4), 041401-041407 (2017).
- Flyvbjerg, H., Petersen, H. G. Error estimates on averages of correlated data. The Journal of Chemical Physics. 91 (1), 461-466 (1989).
- Tuckerman, M. E. Statistical mechanics: theory and molecular simulation. , Oxford University Press. (2010).
- McDonough, W. F., Sun, S. S. The composition of the Earth. Chemical Geology. 120, 223-253 (1995).
- Elkins-Taton, L. T. Magma oceans in the inner solar system. Annual Review of Earth and Planetary Sciences. 40, 113-139 (2012).
- Lock, S. J., et al. The origin of the Moon within a terrestrial synestia. J. Geophysical Research: Planets. 123, 910-951 (2018).
- Solomatova, N. V., Caracas, R. Pressure-induced coordination changes in a pyrolitic silicate melt from ab initio molecular dynamics simulations. Journal of Geophysical Research: Solid Earth. 124, 11232-11250 (2019).
重印与许可
