这里会显示出您选择的修订版和当前版本之间的差别。
两侧同时换到之前的修订记录前一修订版后一修订版 | 前一修订版上一修订版两侧同时换到之后的修订记录 | ||
atk:用分子动力学方法模拟液体的粘度 [2019/09/29 21:05] – [分子粘性] xie.congwei | atk:用分子动力学方法模拟液体的粘度 [2019/09/29 22:39] – [扩展结果] xie.congwei | ||
---|---|---|---|
行 35: | 行 35: | ||
==== 力场 ==== | ==== 力场 ==== | ||
+ | 经典力场的原理是用被设计为近似总能量的简单函数代替计算上消耗过大的电子结构计算。这种方法可以确保有效地模拟更大的系统,远远超出了电子结构计算所能实现的。该方法的主要挑战之一是为替换函数拟合参数,以便它们能够准确地预测能量。 | ||
+ | 一种方法是将分子分成许多化学类型。然后使用这些类型找到不同类型交互作用的特定参数。例如,乙烷中的 sp< | ||
+ | |||
+ | 在本教程中,您将使用 OPLS-AA 力场模拟甲醇< | ||
===== 计算步骤 ===== | ===== 计算步骤 ===== | ||
==== 构建初始分子构型 ==== | ==== 构建初始分子构型 ==== | ||
+ | |||
+ | |||
+ | 打开 **Builder** {{: | ||
+ | |||
+ | {{ : | ||
+ | |||
+ | |||
+ | 现在需要为每个原子指定化学类型。OPLS 甲醇有 4 种类型: | ||
+ | |||
+ | * ' | ||
+ | * ' | ||
+ | * ' | ||
+ | * ' | ||
+ | |||
+ | 选择附着在碳上的三个氢原子,通过 Selection Tools {{: | ||
+ | |||
+ | <WRAP center tip 100%> | ||
+ | === 提示 === | ||
+ | 选中原子时,按住 //Ctrl// 键可将该原子添加到当前的选择中,而按住 //Shift// 键可将其删除。 | ||
+ | </ | ||
+ | |||
+ | 接下来,选择 Builders {{: | ||
+ | |||
+ | {{ : | ||
+ | |||
+ | |||
+ | 将鼠标悬停在新结构中的原子上,检查标签是否反映了每个原子的类型。 | ||
+ | |||
+ | |||
+ | |||
+ | |||
+ | |||
+ | |||
+ | |||
+ | |||
==== 为计算添加力场 ==== | ==== 为计算添加力场 ==== | ||
+ | 如果您的原子标签指定没有问题,点击 **Send To** 按钮 {{: | ||
+ | |||
+ | 首先,在 **Scripter** 里添加一个 **ForceFieldCalculator** 模块。在此模块中需要定义一个力场,因此需要将新的力场定义添加到脚本中。打开 **ForceFieldCalculator** 并选择 New,在打开的 **Potential Editor** 对话框输入新的力场项。 | ||
+ | |||
+ | |||
+ | {{ : | ||
+ | |||
+ | 给新的势命名为 '' | ||
=== 原子标签 === | === 原子标签 === | ||
+ | |||
+ | 点击 Add 打开 Particle Type Editor。在这里添加类型和与之相关的非键合信息。 | ||
+ | |||
+ | OPLS-AA 力场采用以下公式表达非键合能量 $E_{nb}$: | ||
+ | |||
+ | $$E_{nb} = \sum_{ij} 4\epsilon_{ij}\left[\left(\frac{\sigma_{ij}}{r_ij}\right)^{12} | ||
+ | - \left(\frac{\sigma_{ij}}{r_{ij}}\right)^6\right] + \frac{q_iq_j}{4\pi\varepsilon_0 r_{ij}}$$ | ||
+ | |||
+ | $\epsilon_{ij}$ 和 $\sigma_{ij}$ 项由 | ||
+ | |||
+ | $$\epsilon_{ij} = \sqrt{\epsilon_i\epsilon_j}$$ | ||
+ | |||
+ | $$\sigma_{ij} = \frac{\sigma_i+\sigma_j}{2}$$ | ||
+ | |||
+ | 甲醇中 4 种原子类型的非键合参数如下: | ||
+ | |||
+ | ^ Type ^ $\epsilon_{i}$ / kcal mol-1 ^ $\sigma_{i}$ / Å ^ $q_{i}$ / e ^ | ||
+ | | CT | 0.066 | 3.50 | 0.145 | | ||
+ | | OH | 0.170 | 3.12 | -0.683 | ||
+ | | HC | 0.030 | 2.50 | 0.04 | | ||
+ | | HO | 0.000 | 0.00 | 0.418 | | ||
+ | |||
+ | |||
+ | |||
+ | <WRAP center important 100%> | ||
+ | === 注意 === | ||
+ | HO 型不是 OPLS-AA 势中的范德瓦尔斯位点,因此该原子的参数为零。但它有静电场,因此带电荷。 | ||
+ | </ | ||
+ | |||
+ | 从弹出的周期表中选择元素,然后从下拉菜单中选中标签。添加相应的非键合参数。例如,CT 型的输入信息如下所示。 | ||
+ | |||
+ | |||
+ | {{ : | ||
=== 键合 === | === 键合 === | ||
+ | |||
+ | 键伸缩能 $E_{bond}$ 可以用简谐波弹簧方程近似表达为: | ||
+ | |||
+ | $$E_{bond} = \sum_{ij} K_{ij}\left( r_{ij} - r_0 \right)^2$$ | ||
+ | |||
+ | 甲醇中有三种不同的键类型,这些键的参数如下表所示: | ||
+ | |||
+ | ^ Type I ^ **Type J** ^ $K_{ij}$ / kcal mol< | ||
+ | | CT | OH | 320.0 | ||
+ | | CT | HC | 340.0 | ||
+ | | OH | HO | 553.0 | ||
+ | |||
+ | |||
+ | 为添加键合势,在 //Potential Components// | ||
+ | |||
+ | |||
+ | {{ : | ||
+ | |||
+ | |||
+ | |||
== 键角 == | == 键角 == | ||
+ | |||
+ | 键角弯折能 $E_{angle}$ 可以用另一个简谐波弹簧方程近似表达为: | ||
+ | |||
+ | $$E_{angle} = \sum_{ijk} K_{ijk}\left( \theta_{ijk} - \theta_0 \right)^2$$ | ||
+ | |||
+ | 甲醇中有三种不同的键角类型,这三种键角参数如下表所示: | ||
+ | |||
+ | ^ Type I ^ Type J ^ Type K ^ $K_{ijk}$/ eV rad< | ||
+ | | CT | OH | HO | 2.3861 | ||
+ | | HC | CT | OH | 1.5184 | ||
+ | | HC | CT | HC | 1.4317 | ||
+ | |||
+ | |||
+ | 为添加键角势,在 //Potential Components// | ||
+ | |||
+ | {{ : | ||
+ | |||
== 扭转 == | == 扭转 == | ||
+ | |||
+ | |||
+ | 在甲醇中,只有一种键扭转项。这将延伸到 '' | ||
+ | |||
+ | $$E_{torsion} = K (1+\cos(3\phi))$$ | ||
+ | |||
+ | OPLS-AA 势中的总扭转项实际上是更为复杂的 3 级傅立叶展开。这个势不能直接在界面上添加,但可以通过少量修改脚本而轻松设置。创建势的步骤,在 //Potential Components// | ||
+ | |||
+ | |||
+ | |||
== 范德瓦尔斯 == | == 范德瓦尔斯 == | ||
+ | |||
+ | 虽然已在粒子类型定义中为单个原子类型添加了参数,但还需要为类型组合添加项。原子类型 //HO// 不是范德瓦尔斯位点,因此不需要添加涉及该原子的配对项。其余 3 种类型可以给出 6 种唯一的组合。 | ||
+ | |||
+ | 对于 //Potential Components// | ||
+ | |||
+ | * 可以通过按下 //Apply Combination Rules// 从独立原子类型参数中填写 Sigma 和 epsilon。 | ||
+ | * 设置 $r_{cut}$ 为 '' | ||
+ | * 设置 $r_{i}$ 为 '' | ||
+ | * 设置 //Bonded mode scaling// 为 '' | ||
+ | |||
+ | 一个 '' | ||
+ | |||
+ | |||
+ | {{ : | ||
+ | |||
+ | |||
+ | |||
行 65: | 行 209: | ||
- | === 弛豫起始结构 === | + | 在粒子类型的定义中将原子电荷添加到模型中,且仍然需要添加静电项的求和法。 |
+ | |||
+ | 在 //Potential Components// | ||
+ | |||
+ | 右侧的参数选择为: | ||
+ | |||
+ | * 设置 $r_{cut}$ 为 '' | ||
+ | * 设置精确度为 '' | ||
+ | |||
+ | 其他参数的默认值适用于本次计算。 | ||
+ | |||
+ | 已经完成添加这些势之后,点击 **Potential Editor** 和 **ForceFieldCalculator** 的 //OK//。 | ||
+ | |||
+ | |||
+ | |||
+ | |||
+ | |||
+ | |||
+ | === 弛豫起始结构 === | ||
== 优化结构 == | == 优化结构 == | ||
+ | |||
+ | 现在已经完成了力场的定义,可以将计算添加到脚本中。 | ||
+ | |||
+ | <WRAP center important 100%> | ||
+ | === 注意 === | ||
+ | 本教程遵循运行分子动力学模拟的一般方法。有关方法和设置细节的更多信息,请参阅以下教程: | ||
+ | |||
+ | * [[https:// | ||
+ | * [[https:// | ||
+ | |||
+ | </ | ||
+ | |||
+ | 第一步是几何优化以从初始构型中消除任何大的力。 | ||
+ | |||
+ | 在 **Scripter** 中添加 Optimization {{: | ||
+ | |||
+ | |||
+ | |||
== 获得合适的密度 == | == 获得合适的密度 == | ||
+ | |||
+ | 构建结构时,密度太低会使甲醇分子不能更轻易地添加到块体结构中。现在晶胞需要收缩到大致合适的密度,这个数据可以通过简单的分子动力学计算获得。 | ||
+ | |||
+ | 添加一个采用 //NPT Berendsen// 恒压器的分子动力学模块。这个恒压器要比 MTK 恒压器更稳定,却没有那么准确,可用于从初始猜想结构中创建弛豫结构。100 ps 的运行通常足以使系统冷凝成液体密度。 | ||
+ | |||
+ | |||
=== 平衡结构 === | === 平衡结构 === | ||
+ | |||
+ | 添加一个采用 //NPT Martyna Tobias Klein// 恒压器的分子动力学模块。在本次计算中,时间设置为 100 ps 可以保证模拟达到平衡。平衡时间可能会有很大差异,取决于所研究的系统。模拟还应该以之前计算的速度开始,可以通过设置 //Initial Velocity// 类型为 '' | ||
+ | |||
+ | |||
=== 设置生成计算 === | === 设置生成计算 === | ||
+ | |||
+ | 使用 //NPT Martyna Tobias Klein// 恒压器添加一个分子动力学块。将动力学模拟运行设置为 500 ps。生成模拟适当的长度在很大程度上取决于所研究的系统以及最终结果的期望精度。在这部分中应该要设置的两个特定项为: | ||
+ | |||
+ | |||
+ | * 在 //Save Trajectory// | ||
+ | * 在 //Log interval// 框设置 1000 | ||
+ | |||
+ | 这两个细节的设置都将用于本教程的分析部分。 | ||
+ | |||
+ | |||
+ | |||
行 81: | 行 282: | ||
== 添加扭转项 == | == 添加扭转项 == | ||
+ | |||
+ | 现在基本的计算已经创建完成,在 **Scripter** 设置输出文件为 '' | ||
+ | |||
+ | {{ : | ||
+ | |||
+ | |||
+ | 按下 **Send To** 按钮 {{: | ||
+ | |||
+ | <code python> | ||
+ | potentialSet.addPotential(_potential) | ||
+ | _potential = CosineTorsionPotential( | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | k = [0.0*eV, ], | ||
+ | n = [0, ], | ||
+ | delta = [0.0*Radians, | ||
+ | ) | ||
+ | </ | ||
+ | |||
+ | 使用以下代码块替换上面的部分 | ||
+ | |||
+ | <code python> | ||
+ | potentialSet.addPotential(_potential) | ||
+ | _potential = CosineTorsionPotential( | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | k = [0.176*kiloCaloriePerMol, | ||
+ | n = [3, ], | ||
+ | delta = [0.0*Radians, | ||
+ | ) | ||
+ | </ | ||
+ | |||
+ | 保存编辑过的脚本。 | ||
== 定义 hook 函数 == | == 定义 hook 函数 == | ||
+ | |||
+ | 为了计算粘度,需要得到模拟期间压力张量的值。由于模拟过程中压力变化迅速,因此数据的高精度是必要的。虽然将结构保存至轨迹时也保存了压力张量,但是使用短的时间步长保存轨迹文件会导致文件太大超出合理范围。我们所需要的是采用一种方法只保存每一步的压力信息。 | ||
+ | |||
+ | 在 **QuantumATK** 中,上述问题可以通过 hook 函数轻松解决。这些函数会在每次分子动力学整合之前或之后运行。hook 函数被设计用于控制模拟,或是在这种情况下,记录和分析模拟器件的不同量值。 | ||
+ | |||
+ | <WRAP center important 100%> | ||
+ | === 注意 === | ||
+ | 有关如何使用 hook 函数修改模拟的示例,请参阅教程:[[https:// | ||
+ | </ | ||
+ | |||
+ | 为使用 hook 函数,必须首先添加一个实现 ‘__call__’ 函数的 Python 对象。传递给这个函数的参数是: | ||
+ | |||
+ | * // | ||
+ | * // | ||
+ | * // | ||
+ | * // | ||
+ | * // | ||
+ | |||
+ | 利用这些数值可以创建许多自定义分析和数据记录方法。在当前的案例中,需要记录压力张量,这可以通过下面脚本中定义的 // | ||
+ | |||
+ | <code python> | ||
+ | # Define the class that does the analysis for the MD simulation | ||
+ | class PressureRecordingUtility: | ||
+ | ''' | ||
+ | |||
+ | def __init__(self): | ||
+ | ''' | ||
+ | self._time_list = [] | ||
+ | self._volume_list = [] | ||
+ | self._pressure_list = [] | ||
+ | |||
+ | def __call__(self, | ||
+ | ''' | ||
+ | |||
+ | # Add the current time onto the list | ||
+ | self._time_list.append(time) | ||
+ | |||
+ | # Get the cell volume and add that to the list too | ||
+ | volume = configuration.bravaisLattice().unitCellVolume() | ||
+ | self._volume_list.append(volume) | ||
+ | |||
+ | # Get the velocity component of the pressure tensor | ||
+ | velocity = configuration.velocities() | ||
+ | masses = configuration.atomicMasses() | ||
+ | velocity_tensor = ((velocity.reshape(-1, | ||
+ | |||
+ | # Construct the pressure tensor | ||
+ | pressure_tensor = (velocity_tensor/ | ||
+ | |||
+ | # Add the pressure tensor to the list | ||
+ | self._pressure_list.append(pressure_tensor) | ||
+ | |||
+ | def times(self): | ||
+ | ''' | ||
+ | return Units.PhysicalQuantity(self._time_list) | ||
+ | |||
+ | def volumes(self): | ||
+ | ''' | ||
+ | return Units.PhysicalQuantity(self._volume_list) | ||
+ | |||
+ | def pressure_tensors(self): | ||
+ | ''' | ||
+ | return Units.PhysicalQuantity(self._pressure_list) | ||
+ | </ | ||
+ | |||
+ | 将此代码粘贴到模拟脚本的顶部,然后保存编辑的脚本。 | ||
+ | |||
+ | |||
+ | |||
+ | |||
+ | |||
== 在分子动力学模拟中使用 hook 函数 == | == 在分子动力学模拟中使用 hook 函数 == | ||
+ | |||
+ | 现在已定义了适合的 hook 函数,接下来需要告诉 **QuantumATK** 如何在模拟过程中使用这个函数。 | ||
+ | |||
+ | 找到脚本中执行分子动力学计算的最终生成阶段的部分。在该部分的顶端添加以下行: | ||
+ | |||
+ | <code python> | ||
+ | pressureRecorder = PressureRecordingUtility() | ||
+ | </ | ||
+ | |||
+ | 该行代码声明了用于记录压力数据的类的一个实例,分配用于存储压力数据的存储空间。 | ||
+ | |||
+ | 现在需要更改 **MolecularDynamics** 的调用以添加 hook 函数。编辑此函数调用以包含参数: | ||
+ | |||
+ | <code python> | ||
+ | post_step_hook=pressureRecorder | ||
+ | </ | ||
+ | |||
+ | 这行代码告诉 **QuantumATK** 使用 // | ||
+ | |||
+ | |||
+ | <WRAP center important 100%> | ||
+ | == 注意 == | ||
+ | **MolecularDynamics** 函数支持参数 // | ||
+ | </ | ||
+ | |||
+ | 最后,需要将 // | ||
+ | |||
+ | <code python> | ||
+ | data_file = ' | ||
+ | temperature = 300 | ||
+ | |||
+ | nlsave(data_file, | ||
+ | nlsave(data_file, | ||
+ | nlsave(data_file, | ||
+ | nlsave(data_file, | ||
+ | nlsave(data_file, | ||
+ | nlsave(data_file, | ||
+ | </ | ||
+ | |||
+ | 保存编辑后的脚本。 | ||
+ | |||
+ | |||
+ | |||
+ | |||
+ | |||
+ | |||
=== 运行计算 === | === 运行计算 === | ||
+ | |||
+ | 脚本中的所有必要部分已准备就位。点击 {{: | ||
+ | |||
==== 结果分析 ==== | ==== 结果分析 ==== | ||
+ | |||
=== 观察压力变化 === | === 观察压力变化 === | ||
+ | |||
+ | 采用 NPT 系综时,检验模拟过程中压力的控制方式非常重要。理想情况下,压力应在模拟过程中最终达到目标值并趋于平衡。如果这种情况没有发生,则可能需要更长的平衡或模拟时间。 | ||
+ | |||
+ | 在 NPT 系综中,默认情况下会将压力写入 HDF5 文件。为了将压力可视化,可以使用脚本从 HDF5 文件中读取数据并绘图。为此,请点击 {{: | ||
+ | |||
+ | <code python> | ||
+ | import pylab | ||
+ | |||
+ | # Read in the existing trajectory to ge the da{{ : | ||
+ | trajectory_file " | ||
+ | trajectory = nlread(trajectory_file, | ||
+ | time = trajectory.times() | ||
+ | pressure = trajectory.pressures() | ||
+ | |||
+ | # Set the average and target pressures | ||
+ | p_tot_avg = numpy.cumsum(pressure.inUnitsOf(GPa)) * GPa / (numpy.arange(pressure.shape[0])+1) | ||
+ | p_target = numpy.ones(pressure.shape[0]) * Pa | ||
+ | |||
+ | # Plot the resulting pressure and average pressure | ||
+ | pylab.figure() | ||
+ | pylab.plot(time.inUnitsOf(picosecond), | ||
+ | pylab.plot(time.inUnitsOf(picosecond), | ||
+ | pylab.plot(time.inUnitsOf(picosecond), | ||
+ | pylab.xlabel(' | ||
+ | pylab.ylabel(' | ||
+ | pylab.legend() | ||
+ | pylab.show() | ||
+ | </ | ||
+ | |||
+ | 也可以直接在此处下载脚本 [[https:// | ||
+ | |||
+ | 将脚本发送到 **Job Manager** 运行。脚本应该显示与以下相似的图。 | ||
+ | |||
+ | <WRAP center tip 100%> | ||
+ | === 提示 === | ||
+ | 也可以使用 atkpython 命令从命令行运行像这样的分析脚本。 | ||
+ | </ | ||
+ | |||
+ | {{ : | ||
+ | |||
+ | 请注意,在动力学模拟期间,压力摆动约 0.2 GPa。这些大的压力变化是分子动力学模拟的典型特征。重要的是,模拟运行时平均压力收敛到目标 0.1 MPa。 | ||
+ | |||
+ | |||
+ | |||
+ | |||
+ | |||
+ | |||
=== 由 Einstein 关系式计算粘度 === | === 由 Einstein 关系式计算粘度 === | ||
+ | |||
+ | 在分子动力学模拟期间使用 hook 函数收集的数据,可用于计算和可视化粘度的估值。 | ||
+ | |||
+ | 以上可以通过使用 **ATKPython** 脚本来完成。在 **Editor** 中打开一个新脚本并粘贴以下脚本。 | ||
+ | |||
+ | <code python> | ||
+ | import pylab | ||
+ | import scipy | ||
+ | |||
+ | # Read the data that is stored in the HDF5 file from the simulation | ||
+ | data_file = sys.argv[1] | ||
+ | time = nlread(data_file, | ||
+ | pressure_tensor = nlread(data_file, | ||
+ | volume_avg = nlread(data_file, | ||
+ | temperature = nlread(data_file, | ||
+ | time_step = nlread(data_file, | ||
+ | |||
+ | # Set up the calculation to create 100 time based estimates of the viscosity | ||
+ | N = pressure_tensor.shape[0] | ||
+ | N_steps = 101 | ||
+ | skip = int(N/100) | ||
+ | time_skip = time[:: | ||
+ | |||
+ | # Calculate the off-diagonal elements of the pressure tensor | ||
+ | P_shear = numpy.zeros((5, | ||
+ | P_shear[0] = pressure_tensor[:, | ||
+ | P_shear[1] = pressure_tensor[:, | ||
+ | P_shear[2] = pressure_tensor[:, | ||
+ | P_shear[3] = (pressure_tensor[:, | ||
+ | P_shear[4] = (pressure_tensor[:, | ||
+ | |||
+ | # At increasing time lengths, calculate the viscosity based on that part of the simulatino | ||
+ | pressure_integral = numpy.zeros(N_steps, | ||
+ | for t in range(1, | ||
+ | total_step = t*skip | ||
+ | |||
+ | for i in range(5): | ||
+ | integral = scipy.integrate.trapz( | ||
+ | y = P_shear[i][: | ||
+ | dx=time_step.inUnitsOf(Second) | ||
+ | ) | ||
+ | integral *= Second * Joule / Meter**3 | ||
+ | pressure_integral[t] += integral**2 / 5 | ||
+ | |||
+ | # Finally calculate the overall viscosity | ||
+ | # Note that here the first step is skipped to avoid divide by zero issues | ||
+ | kbT = boltzmann_constant * temperature | ||
+ | viscosity = pressure_integral[1: | ||
+ | |||
+ | # Print the final viscosity | ||
+ | print(" | ||
+ | |||
+ | # Display the evolution of the viscosity in time | ||
+ | pylab.figure() | ||
+ | pylab.plot(time_skip[1: | ||
+ | pylab.xlabel(' | ||
+ | pylab.ylabel(' | ||
+ | pylab.legend() | ||
+ | pylab.show() | ||
+ | </ | ||
+ | |||
+ | 也可以直接在此处下载脚本 [[https:// | ||
+ | |||
+ | 脚本读取由 hook 函数在 HDF5 文件汇总保存的数据。然后利用在理论部分给出的 Einstein 关系式计算粘度。注意,积分在压力张量的 5 个独特的非对角线分量 $P_{xy}$、$P_{yz}$、$P_{xz}$ $\frac{P_{xx}-P_{yy}}{2}$、$\frac{P_{yy}-P_{zz}}{2}$ 上取平均值。 | ||
+ | |||
+ | 然后,该脚本在模拟时间限制下输出最终粘度,并显示在模拟过程中粘度的估值如何变化。 | ||
+ | |||
+ | 以下是通过模拟计算的粘度的典型示例。为了进行比较,298K 下甲醇密度的实验值为 0.54 cP。刚性 OPLS-AA 模型仅包含非键合相互作用并保持每个分子的构型固定,粘度为 0.43 cP< | ||
+ | |||
+ | |||
+ | {{ : | ||
+ | |||
+ | |||
+ | |||
=== 由 Green-Kubo 关系式计算粘度 === | === 由 Green-Kubo 关系式计算粘度 === | ||
+ | 采用 Green-Kubo 关系可以获得类似的粘度结果。您可以在此处下载概述如何实现此目标的脚本 [[https:// | ||
+ | 我们鼓励 Python 的高级用户查看和运行脚本以了解其工作原理。该脚本还演示了 **ATKPython** 平台的一些强大功能和灵活性,不止运行,还可以分析模拟。 | ||
==== 扩展结果 ==== | ==== 扩展结果 ==== | ||
行 100: | 行 582: | ||
=== 运行更多独立轨迹 === | === 运行更多独立轨迹 === | ||
+ | 在上一节中,我们了解到为了提高估算粘度的准确性,需要进行更多的模拟。 | ||
- | === 提高准确性 === | + | 在 **QuantumATK** 中,通过将它们放在 Python 循环中,可以轻松地重复模拟。在甲醇的例子中,一种策略是循环模拟的平衡和产出阶段。随着速度在平衡阶段开始时的重新随机化,将会产生一系列独立的轨迹。由此,采用标准统计方法,平均粘度及其误差都可以轻易估算。 |
+ | 在 Python 中循环计算,可以将要重复的代码放在 //for// 循环中。为了实现粘度计算循环 30 次以产生不同的轨迹,需要将模拟代码放置在 //for// 循环中,例如: | ||
+ | |||
+ | <code python> | ||
+ | for loop in range(30): | ||
+ | # Molecular dynamics simulation commands | ||
+ | |||
+ | # Futher commands outside the loop | ||
+ | </ | ||
+ | |||
+ | <WRAP center important 100%> | ||
+ | === 注意 === | ||
+ | Python 的一个重要特性是它使用缩进来确定代码块的范围。这就确定了哪些行在循环内部或外部以及逻辑语句。上面示例中的缩进即展示了如何区分 Python 循环内部或外部的代码。 | ||
+ | </ | ||
+ | |||
+ | |||
+ | 创建多个轨迹的脚本中唯一必要的变化是修改输出的完成方式。在 HDF5 中,使用相同文本描述符写入数据的文件会覆盖以前的数据,因此需要唯一的描述符。这可以通过将循环变量(在上面的例子中变量// | ||
+ | |||
+ | <code python> | ||
+ | data_file = ' | ||
+ | temperature = 300 | ||
+ | |||
+ | nlsave(data_file, | ||
+ | nlsave(data_file, | ||
+ | nlsave(data_file, | ||
+ | nlsave(data_file, | ||
+ | nlsave(data_file, | ||
+ | nlsave(data_file, | ||
+ | |||
+ | 为读取每条轨迹的数据进行分析,需要以相同的方式利用 '' | ||
+ | |||
+ | |||
+ | |||
+ | |||
+ | |||
+ | === 提高准确性 === | ||
+ | 为了演示如何使用多次运行提高模拟的准确性,模拟重复了 30 次。从这些模拟结构得到的平均值为 0.41 cP,标准误差为 0.05 cP。下图显示了平均粘度的变化和标准误差与模拟次数的关系。 | ||
+ | {{ : | ||
+ | ===== 参考 ===== | ||