费米维基

用户工具

站点工具

您的足迹:
atk:用分子动力学方法模拟液体的粘度

差别

这里会显示出您选择的修订版和当前版本之间的差别。

到此差别页面的链接

两侧同时换到之前的修订记录前一修订版
后一修订版		前一修订版
atk:用分子动力学方法模拟液体的粘度 [2019/09/29 22:14] – [为计算添加力场] xie.congweiatk:用分子动力学方法模拟液体的粘度 [2019/09/29 22:43] (当前版本) – [参考] xie.congwei
行 448: 行 448:
  
 ==== 结果分析 ==== ==== 结果分析 ====
 +
 === 观察压力变化 === === 观察压力变化 ===
 +
 +采用 NPT 系综时，检验模拟过程中压力的控制方式非常重要。理想情况下，压力应在模拟过程中最终达到目标值并趋于平衡。如果这种情况没有发生，则可能需要更长的平衡或模拟时间。
 +
 +在 NPT 系综中，默认情况下会将压力写入 HDF5 文件。为了将压力可视化，可以使用脚本从 HDF5 文件中读取数据并绘图。为此，请点击 {{:atk:editor.png?direct&25|}} 打开 **Editor**，粘贴以下脚本。
 +
 +<code python>  
 +import pylab
 +
 +# Read in the existing trajectory to ge the da{{ :undefined:methanol_pressure_graph_new-20190929.png |}}ta
 +trajectory_file "Production_Trajectory.hdf5"
 +trajectory = nlread(trajectory_file, MDTrajectory)[-1]
 +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), pressure.inUnitsOf(GPa), label='Pressure (GPa)', linewidth=0.5)
 +pylab.plot(time.inUnitsOf(picosecond), p_tot_avg.inUnitsOf(GPa), label='P Average', linewidth=2.0)
 +pylab.plot(time.inUnitsOf(picosecond), p_target.inUnitsOf(GPa), label='P Target', linewidth=2.0)
 +pylab.xlabel('Time (ps)')
 +pylab.ylabel('Pressure (GPa)')
 +pylab.legend()
 +pylab.show()
 +</code>
 +
 +也可以直接在此处下载脚本 [[https://docs.quantumatk.com/_downloads/pressure_analysis.py|↓ pressure_analysis.py]]
 +
 +将脚本发送到 **Job Manager** 运行。脚本应该显示与以下相似的图。
 +
 +<WRAP center tip 100%>
 +=== 提示 ===
 +也可以使用 atkpython 命令从命令行运行像这样的分析脚本。
 +</WRAP>
 +
 +{{ :atk:methanol_pressure_graph_new-20190929.png |}}
 +
 +请注意，在动力学模拟期间，压力摆动约 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, object_id="Time_8")[-1]
 +pressure_tensor = nlread(data_file, object_id="Pressure_8")[-1]
 +volume_avg = nlread(data_file, object_id="Volume_Average_8")[-1]
 +temperature = nlread(data_file, object_id="Temperature_8")[-1]
 +time_step = nlread(data_file, object_id="Time_Step_8")[-1]
 +
 +# 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[::skip]
 +
 +# Calculate the off-diagonal elements of the pressure tensor
 +P_shear = numpy.zeros((5,N), dtype=float) * Joule / Meter**3
 +P_shear[0] = pressure_tensor[:,0,1]
 +P_shear[1] = pressure_tensor[:,0,2]
 +P_shear[2] = pressure_tensor[:,1,2]
 +P_shear[3] = (pressure_tensor[:,0,0] - pressure_tensor[:,1,1]) / 2
 +P_shear[4] = (pressure_tensor[:,1,1] - pressure_tensor[:,2,2]) / 2
 +
 +# At increasing time lengths, calculate the viscosity based on that part of the simulatino
 +pressure_integral = numpy.zeros(N_steps, dtype=numpy.float) * (Second * Joule / Meter**3)**2
 +for t in range(1,N_steps):
 + total_step = t*skip
 +
 + for i in range(5):
 + integral = scipy.integrate.trapz(
 + y = P_shear[i][:total_step].inUnitsOf(Joule/Meter**3),
 + 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:] * volume_avg / (2*kbT*time_skip[1:])
 +
 +# Print the final viscosity
 +print("Viscosity is {} cP".format(viscosity[-1].inUnitsOf(millisecond*Pa)))
 +
 +# Display the evolution of the viscosity in time
 +pylab.figure()
 +pylab.plot(time_skip[1:].inUnitsOf(ps), viscosity.inUnitsOf(millisecond*Pa), label='Viscosity')
 +pylab.xlabel('Time (ps)')
 +pylab.ylabel('Viscosity (cP)')
 +pylab.legend()
 +pylab.show()
 +</code>
 +
 +也可以直接在此处下载脚本 [[https://docs.quantumatk.com/_downloads/viscosity_analysis.py|↓ viscosity_analysis.py]]
 +
 +脚本读取由 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<sup>[3]</sup>。尽管该示例给出了类似的结果，但在计算的粘度中仍存在合理的噪声。理想地，该粘度应该收敛到一个不同的值。这表明需要更多的压力取样来提高模拟的准确性。
 +
 +
 +{{ :atk:methanol_einstein_progress-20190929.png?600 |}}
 +
 +
 +
 === 由 Green-Kubo 关系式计算粘度 === === 由 Green-Kubo 关系式计算粘度 ===
  
 +采用 Green-Kubo 关系可以获得类似的粘度结果。您可以在此处下载概述如何实现此目标的脚本 [[https://docs.quantumatk.com/_downloads/green_kubo_analysis.py|↓ green_kubo_analysis.py]]。
  
 +我们鼓励 Python 的高级用户查看和运行脚本以了解其工作原理。该脚本还演示了 **ATKPython** 平台的一些强大功能和灵活性，不止运行，还可以分析模拟。
  
 ==== 扩展结果 ==== ==== 扩展结果 ====
行 458: 行 582:
 === 运行更多独立轨迹 === === 运行更多独立轨迹 ===
  
 +在上一节中，我们了解到为了提高估算粘度的准确性，需要进行更多的模拟。
  
-=== 提高准确性 ===+在 **QuantumATK** 中，通过将它们放在 Python 循环中，可以轻松地重复模拟。在甲醇的例子中，一种策略是循环模拟的平衡和产出阶段。随着速度在平衡阶段开始时的重新随机化，将会产生一系列独立的轨迹。由此，采用标准统计方法，平均粘度及其误差都可以轻易估算。
  
 +在 Python 中循环计算，可以将要重复的代码放在 //for// 循环中。为了实现粘度计算循环 30 次以产生不同的轨迹，需要将模拟代码放置在 //for// 循环中，例如：
  
 +<code python>  
 +for loop in range(30):
 +    # Molecular dynamics simulation commands
 +
 +# Futher commands outside the loop
 +</code>
 +
 +<WRAP center important 100%>
 +=== 注意 ===
 +Python 的一个重要特性是它使用缩进来确定代码块的范围。这就确定了哪些行在循环内部或外部以及逻辑语句。上面示例中的缩进即展示了如何区分 Python 循环内部或外部的代码。
 +</WRAP>
 +
 +
 +创建多个轨迹的脚本中唯一必要的变化是修改输出的完成方式。在 HDF5 中，使用相同文本描述符写入数据的文件会覆盖以前的数据，因此需要唯一的描述符。这可以通过将循环变量（在上面的例子中变量//循环//）包含在名称中实现。更改文件写入代码为：
 +
 +<code python>  
 +data_file = 'md_methanol_data.hdf5'
 +temperature = 300
 +
 +nlsave(data_file, pressureRecorder.times(), "Time_"+str(loop))
 +nlsave(data_file, pressureRecorder.pressure_tensors(), "Pressure_"+str(loop))
 +nlsave(data_file, pressureRecorder.volumes(), "Volume_"+str(loop))
 +nlsave(data_file, pressureRecorder.volumes().mean(), "Volume_Average_"+str(loop))
 +nlsave(data_file, temperature, "Temperature_"+str(loop))
 +nlsave(data_file, md_trajectory.timeStep(), "Time_Step_"+str(loop))</code>
 +
 +为读取每条轨迹的数据进行分析，需要以相同的方式利用 ''nlread'' 命令更改名称。
 +
 +
 +
 +
 +
 +=== 提高准确性 ===
  
  
 +为了演示如何使用多次运行提高模拟的准确性，模拟重复了 30 次。从这些模拟结构得到的平均值为 0.41 cP，标准误差为 0.05 cP。下图显示了平均粘度的变化和标准误差与模拟次数的关系。
  
 +{{ :atk:methanol_results-20190929.png?600 |}}
  
 +===== 参考 =====
 +  * [1] B. Hess: Determining the shear viscosity of model liquids from molecular dynamics simulations. [[https://aip.scitation.org/doi/abs/10.1063/1.1421362|J. Chem. Phys. 116, 209 (2016)]]
 +  * [2] S. H. Jamali, R. Hartkamp, C. Bardas, J. Sohl, T. J. H. Vlugt, O. A. Moultos: Shear viscosity computed from the finite-size effects of self-diffusivity in equilibrium molecular dynamics. [[https://pubs.acs.org/doi/10.1021/acs.jctc.8b00625|J. Chem. Theory. Comput. 14, 5959 (2018)]]
 +  * [3] (1, 2) D. Gonzalez-Salgado, C. Vega: A new intermolecular potential for simulations of methanol: The OPLS/2016 model. [[https://aip.scitation.org/doi/10.1063/1.4958320|J. Chem. Phys. 145, 034508 (2016)]]
 +  * [4] W. L. Jorgensen, D. S. Maxwell, J. Tirado-Rives: Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. [[https://pubs.acs.org/doi/10.1021/ja9621760|J. Am. Chem. Soc. 118, 11225 (1996) 118]]
atk/用分子动力学方法模拟液体的粘度.txt · 最后更改: 2019/09/29 22:43 由 xie.congwei

页面工具

© 2014-2022 费米科技（京ICP备14023855号）