使用 UMD 软件包分析 Ab Initio 分子动力学模拟中的熔体和流体

1. 分子动力学运行分析 注意：该软件包可通过 GitHub 网站（https://github.com/rcaracas/UMD_package）和 ERC IMPACT 项目的专用页面（http://moonimpact.eu/umd-package/）作为开放获取软件包获得。 使用包中的一个或多个专用 Python 脚本提取每组特定的物理属性。在命令行运行所有脚本;它们都使用一系列标志，这些标志在一个脚本到另一个脚本之间尽可能一致。 表 1 总结了标志、其含义和默认值。 旗 意义 使用它的脚本 默认值 -h 简短帮助 都 -f UMD 文件名 都 -i 要丢弃的热化步骤 都 0 -i 包含原子间键的输入文件 物种形成 债券.输入 -秒 频率采样 MSD， 物种形成 1（考虑每一步） -一 原子或阴离子列表 物种形成 -c 阳离子列表 物种形成 -l 键长 物种形成 2 -t 温度 振动、流变学 -v 轨迹采样窗口宽度的离散化，用于均方位移分析 MSD 20 -z 轨迹采样窗口开始的离散化，用于均方位移分析 MSD 20 表 1：UMD 包中使用的最常见标志及其最常见的意义。 首先将在第一原则代码（如 VASP8 或 QBox9）中执行的 MD 仿真输出转换为 UMD 文件。 如果 MD 模拟是在 VASP 中完成的，则在命令行中键入：VaspParser.py -f -i 其中 –f 标志定义 VASP OUTCAR 文件的名称，–i 定义热化长度。注意：由 –i 定义的初始步骤允许丢弃表示热化的模拟的第一步。在典型的分子动力学运行中，计算的第一部分表示热化，即系统所有原子描述温度的高斯式分布所需的时间，以及整个系统围绕平衡值表现出温度，压力，能量等的波动。在分析流体的统计特性时，不应考虑模拟的这一热化部分。 转换 .将 文件转换为 .xyz 文件，便于在各种其他软件包（如 VMD4 或 Vesta5）上进行可视化。在命令行类型处：umd2xyz.py -f -i -s 其中 –f 定义 的名称。umd 文件 –i 定义了要丢弃的热化周期，–s 定义了存储在 中的轨迹的采样频率。umd 文件。默认值为 –i 0 –s 1，即考虑模拟的所有步骤，而不会丢弃任何步骤。 使用 umd2poscar.py 脚本将 umd 文件反转为 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 -d -i 注：默认值为Sampling_Frequency（对轨迹进行采样的频率）= 1步;离散化间距（用于绘制 g（r）） = 0.01 Å;初始步骤（轨迹开始时丢弃的步骤数）= 0。径向PDF，gᴀʙ（r）是半径为r的球形壳内距离d_ᴀʙ的B型原子的平均数量，厚度dr以A型原子为中心（图3）：ρ为原子密度，NA和NB为A型和B型原子数，δ（r−rᴀʙ）如果原子A和B位于r和r + dr之间的距离，则delta函数等于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 脚本需要 bond.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）。 打开两个文件。人口.dat 和.统计.dat 连续;它们构成了物种形成脚本的输出。每个簇都写在一行上，指定其化学式，形成时间，死亡时间，寿命，矩阵，其中包含形成该簇的原子列表。绘制模拟中发现的所有化学物质的每个原子簇的寿命，如 .popul.dat 文件所示（图 5）。 绘制种群分析与每个物种的丰度，如 所示。统计.dat 文件。该分析，绝对和相对，对应于情况-r 0的配位多面体的实际统计量;对于聚合的情况，对于-r 1，这需要仔细处理，因为可能需要对相对数量的原子施加一些归一化。丰度对应于一生中的积分。这。stat.dat 文件还列出了每个簇的大小，即形成它的原子数量。 4. 计算扩散系数 提取原子的均方位移（MSD）作为时间的函数，以获得自扩散率。MSD的标准配方为：其中前置因子是重整化。使用MSD工具，有不同的方法来分析流体的动力学方面。注意：T 是模拟的总时间，N α是α类型的原子数。初始时间 t0 是任意的，跨越了模拟的前半部分。Ninit 是初始时间的次数。τ 是计算 MSD 的时间间隔的宽度;其最大值是模拟时间长度的一半。在典型的 MSD 实现中，每个窗口都从前一个窗口的末尾开始。但是，稀疏采样可以加快MSD的计算速度，而不会改变MSD的斜率。为此，第 i 个窗口从时间 t0（i） 开始，但 （i+1） -第 5 个窗口从时间 t0（i） + τ + v 开始，其中 v 的值是用户定义的。类似地，窗口的宽度以用户定义的离散步长增加，如下所示：τ（i） = τ（i-1） + z。z（”水平步长”）和v（”垂直步长”）的值为正或零;两者的默认值均为 20。 使用一系列 msd_umd 脚本计算 MSD。它们的输出打印在 .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（图 6）。在 MSD 与时间的对数-对数表示中，确定斜率变化。分开第一部分，通常很短，它代表 弹道 制度，即碰撞后原子速度的守恒。第二个较长的部分表示 扩散 状态，即碰撞后原子速度的散射。 从 MSD 的斜率计算扩散系数为：其中 Z 是自由度数（Z = 2 表示平面上的扩散，Z = 3 表示空间中的扩散），t 是时间步长。 5. 时间相关函数 使用一般公式计算时间相关性函数作为系统惯性的度量：A 可以是各种时间相关变量，例如原子位置，原子速度，应力，极化等，每个变量都通过绿色 – Kubo关系产生不同的物理性质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 是原子质量。 使用类似电子表格的软件从 振动.dat 文件中绘制振动频谱（图7）。确定 ω=0 处的有限值，该值对应于流体的扩散特性和有限频率下频谱的各种峰值。确定每种原子类型对振动谱的参与。注意：原子类型的分解表明，不同的原子具有不同的 ω= 0贡献，对应于它们的扩散系数。光谱的一般形状比相应的实体更平滑，特征更少。 在壳体上，读取振动光谱上的积分，该积分产生每种原子物种的扩散系数。注意：热力学性质可以通过从振动光谱积分获得，但应谨慎使用结果，因为有两个近似值：积分在准谐波近似中有效，这在高温下不一定成立;并且需要丢弃与扩散相对应的光谱的气体样部分。然后，集成应仅在频谱的晶格状部分上进行。但这种分离通常需要几个进一步的后处理步骤和计算14，而这些步骤和计算在本 UMD 包中未涵盖。 运行 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 和 ω 是拟合参数。积分后，粘度的表达式变为： 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 输出是完整的，并打印出几个替代单元。此脚本需要图形支持，因为它会创建一个图形图像，用于检查收敛以估计均值上的误差。