费米维基

用户工具

站点工具

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

差别

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

到此差别页面的链接

两侧同时换到之前的修订记录前一修订版
后一修订版		前一修订版
atk:用分子动力学方法模拟液体的粘度 [2019/09/29 21:06] – [力场] xie.congweiatk:用分子动力学方法模拟液体的粘度 [2019/09/29 22:43] (当前版本) – [参考] xie.congwei
行 44: 行 44:
  
 ==== 构建初始分子构型 ==== ==== 构建初始分子构型 ====
 +
 +
 +打开 **Builder** {{:atk:builder.png?direct&25|}}，选择 Add {{:atk:arrow.png?direct&5|}} From Database。将数据库中甲醇的分子结构添加到 **Stash**，确保通过点击 Databases {{:atk:arrow.png?direct&5|}} Molecules 选择分子数据库。
 +
 +{{ :atk:methanol_database-20190929.png |}}
 +
 +
 +现在需要为每个原子指定化学类型。OPLS 甲醇有 4 种类型：
 +
 +  * 'CT' 为四面体碳
 +  * 'HC' 为附着在碳上的氢
 +  * 'OH' 为乙醇氧
 +  * 'HO' 为乙醇氢
 +
 +选择附着在碳上的三个氢原子，通过 Selection Tools {{:atk:arrow.png?direct&5|}} Tags 将其类型设置为 'HC'，并输入标签 'HC'。您现在应该可以使用标记选择这些原子。对其他三种原子类型重复此过程。
 +
 +<WRAP center tip 100%>
 +=== 提示 ===
 +选中原子时，按住 //Ctrl// 键可将该原子添加到当前的选择中，而按住 //Shift// 键可将其删除。
 +</WRAP>
 +
 +接下来，选择 Builders {{:atk:arrow.png?direct&5|}} Packmol 创建液态甲醇的块体模型。将甲醇分子从 stash 拖入分子类型栏，选择 250 个分子，晶胞长度设置为 28Å。点击 Create。
 +
 +{{ :atk:methanol_packmol-20190929.png?600 |}}
 +
 +
 +将鼠标悬停在新结构中的原子上，检查标签是否反映了每个原子的类型。
 +
 +
 +
 +
 +
 +
 +
 +
  
 ==== 为计算添加力场 ==== ==== 为计算添加力场 ====
  
 +如果您的原子标签指定没有问题，点击 **Send To** 按钮 {{:atk:sendto.png?direct&20|}} 将块体结构发送到 **Scripter** {{:atk:script_generator.png?direct&25|}}。
 +
 +首先，在 **Scripter** 里添加一个 **ForceFieldCalculator** 模块。在此模块中需要定义一个力场，因此需要将新的力场定义添加到脚本中。打开 **ForceFieldCalculator** 并选择 New，在打开的 **Potential Editor** 对话框输入新的力场项。
 +
 +
 +{{ :atk:methanol_potential_window-20190929.png |}}
 +
 +给新的势命名为 ''Methanol OPLS-AA''。打开 **Potential Editor** 时会给每个原子添加粒子类型，需要删除这些并替换粒子类型为 4 种原子标签。选中每个元素并单击 Remove，然后点击 Add 为每个原子添加标签。
  
  
  
 === 原子标签 === === 原子标签 ===
 +
 +点击 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}$ 项由  Lorentz-Berthelot 组合规则计算得到：
 +
 +$$\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 势中的范德瓦尔斯位点，因此该原子的参数为零。但它有静电场，因此带电荷。
 +</WRAP>
 +
 +从弹出的周期表中选择元素，然后从下拉菜单中选中标签。添加相应的非键合参数。例如，CT 型的输入信息如下所示。
 +
 +
 +{{ :atk:methanol_particle_type-20190929.png?600 |}}
  
  
 === 键合 === === 键合 ===
 +
 +键伸缩能 $E_{bond}$ 可以用简谐波弹簧方程近似表达为：
 +
 +$$E_{bond} = \sum_{ij} K_{ij}\left( r_{ij} - r_0 \right)^2$$
 +
 +甲醇中有三种不同的键类型，这些键的参数如下表所示：
 +
 +^  Type I  ^  **Type J**  ^  $K_{ij}$ / kcal mol<sup>-1</sup>Å<sup>-2</sup>  ^  **$r_0$ / Å**  ^
 +|  CT      |  OH          |  320.0                                           |  1.41           |
 +|  CT      |  HC          |  340.0                                           |  1.09           |
 +|  OH      |  HO          |  553.0                                           |  0.945          |
 +
 +
 +为添加键合势，在 //Potential Components// 框中点击 Add {{:atk:arrow.png?direct&5|}} Add potential component。从势类型的列表中选择 //HarmonicBondPotential//。选中一个势，会在右边显示它的相关参数。添加三个势并为每个键型输入参数。''CT-OH'' 键型的输入信息如下所示。
 +
 +
 +{{ :atk:methanol_bond_potential-20190929.png |}}
 +
 +
 +
  
  
 == 键角 == == 键角 ==
 +
 +键角弯折能 $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<sup>-2</sup>  |  $\theta_0$ / degree  |
 +|  CT      |  OH      |  HO      |  2.3861                          |  108.5                |
 +|  HC      |  CT      |  OH      |  1.5184                          |  109.5                |
 +|  HC      |  CT      |  HC      |  1.4317                          |  107.8                |
 +
 +
 +为添加键角势，在 //Potential Components// 框中点击 Add {{:atk:arrow.png?direct&5|}} Add potential component，然后从势类型列表中选择 //HarmonicAnglePotential//。选中一个势就会在右侧显示参数。添加三个势并输入每个键型的参数。''CT - OH - HO'' 键角类型的输入如下所示。
 +
 +{{ :atk:methanol_angle_potential-20190929.png |}}
 +
  
  
 == 扭转 == == 扭转 ==
 +
 +
 +在甲醇中，只有一种键扭转项。这将延伸到 ''HC - CT - OH - HO'' 类型。该项中的键角扭转能可以通过 Fourier 展开近似表达为：
 +
 +$$E_{torsion} = K (1+\cos(3\phi))$$
 +
 +OPLS-AA 势中的总扭转项实际上是更为复杂的 3 级傅立叶展开。这个势不能直接在界面上添加，但可以通过少量修改脚本而轻松设置。创建势的步骤，在 //Potential Components// 框中点击 Add {{:atk:arrow.png?direct&5|}} Add potential component，从势类型列表中选择 //CosineTorsionPotential//。稍后我们将在编辑最终脚本时设置参数。
 +
 +
 +
  
  
  
 == 范德瓦尔斯 == == 范德瓦尔斯 ==
 +
 +虽然已在粒子类型定义中为单个原子类型添加了参数，但还需要为类型组合添加项。原子类型 //HO// 不是范德瓦尔斯位点，因此不需要添加涉及该原子的配对项。其余 3 种类型可以给出 6 种唯一的组合。
 +
 +对于 //Potential Components// 框中的每个组合，点击 Add {{:atk:arrow.png?direct&5|}} Add potential component，从势类型列表中选择 //LennardJonesSplinePotential//。每一种组合的参数现在可以在右侧设置。
 +
 +  * 可以通过按下 //Apply Combination Rules// 从独立原子类型参数中填写 Sigma 和 epsilon。
 +  * 设置 $r_{cut}$ 为 ''10'' Å。这个参数规定了超过该距离就会忽略势。
 +  * 设置 $r_{i}$ 为 ''9'' Å。这个参数规定了超过该距离势会通过样条函数按比例缩小以使得它们在截断处为零。这是很有必要的，因此当原子在彼此的范围内时，力不会突然增加。
 +  * 设置 //Bonded mode scaling// 为 ''0.5''。这种缩放由 OPLS-AA 势指定以使原子连接 3 个键。
 +
 +一个 ''CT - OH'' 非键合势的例子如下所示。
 +
 +
 +{{ :atk:methanol_nonbond_potential-20190929.png |}}
 +
 +
 +
  
  
行 69: 行 209:
  
  
-=== 弛豫起始结构 ===+在粒子类型的定义中将原子电荷添加到模型中，且仍然需要添加静电项的求和法。
  
 +
 +在 //Potential Components// 框中点击 Add {{:atk:arrow.png?direct&5|}} Add potential component，从势类型列表中选择 //CoulombSPME//。这就选择了静电项的 Smooth Particle Mesh Ewald 求和。为了正确估算长程静电作用，本方法采用了 Fourier 转换。
 +
 +右侧的参数选择为：
 +
 +  * 设置 $r_{cut}$ 为 ''9'' Å。
 +  * 设置精确度为 ''0.001''。
 +
 +其他参数的默认值适用于本次计算。
 +
 +已经完成添加这些势之后，点击 **Potential Editor** 和 **ForceFieldCalculator** 的 //OK//。
 +
 +
 +
 +
 +
 +
 +=== 弛豫起始结构 ===
  
 == 优化结构 == == 优化结构 ==
 +
 +现在已经完成了力场的定义，可以将计算添加到脚本中。
 +
 +<WRAP center important 100%>
 +=== 注意 ===
 +本教程遵循运行分子动力学模拟的一般方法。有关方法和设置细节的更多信息，请参阅以下教程：
 +
 +  * [[https://docs.quantumatk.com/tutorials/md_basics/md_basics.html#basic-md-tutorial|Molecular Dynamics: Basics]]
 +  * [[https://docs.quantumatk.com/tutorials/diffusion_liquid_copper/diffusion_liquid_copper.html#diffusion-liquid-copper|Diffusion in Liquids from Molecular Dynamics Simulations]]
 +
 +</WRAP>
 +
 +第一步是几何优化以从初始构型中消除任何大的力。
 +
 +在 **Scripter** 中添加 Optimization {{:atk:arrow.png?direct&5|}} OptimizeGeometry 模块。此模块的默认参数就已足够满足计算需求。
 +
 +
 +
  
 == 获得合适的密度 == == 获得合适的密度 ==
 +
 +构建结构时，密度太低会使甲醇分子不能更轻易地添加到块体结构中。现在晶胞需要收缩到大致合适的密度，这个数据可以通过简单的分子动力学计算获得。
 +
 +添加一个采用 //NPT Berendsen// 恒压器的分子动力学模块。这个恒压器要比 MTK 恒压器更稳定，却没有那么准确，可用于从初始猜想结构中创建弛豫结构。100 ps 的运行通常足以使系统冷凝成液体密度。
 +
 +
  
 === 平衡结构 === === 平衡结构 ===
 +
 +添加一个采用 //NPT Martyna Tobias Klein// 恒压器的分子动力学模块。在本次计算中，时间设置为 100 ps 可以保证模拟达到平衡。平衡时间可能会有很大差异，取决于所研究的系统。模拟还应该以之前计算的速度开始，可以通过设置 //Initial Velocity// 类型为 ''Configuration Velocities'' 实现。
 +
 +
  
 === 设置生成计算 === === 设置生成计算 ===
 +
 +使用 //NPT Martyna Tobias Klein// 恒压器添加一个分子动力学块。将动力学模拟运行设置为 500 ps。生成模拟适当的长度在很大程度上取决于所研究的系统以及最终结果的期望精度。在这部分中应该要设置的两个特定项为：
 +
 +
 +  * 在 //Save Trajectory// 框设置 ''Production_Trajectory.hdf5''
 +  * 在 //Log interval// 框设置 1000
 +
 +这两个细节的设置都将用于本教程的分析部分。
 +
 +
 +
  
  
行 85: 行 282:
  
 == 添加扭转项 == == 添加扭转项 ==
 +
 +现在基本的计算已经创建完成，在 **Scripter** 设置输出文件为 ''Methanol_Viscosity.hdf5''。脚本生成后应该如下所示：
 +
 +{{ :atk:methanol_scriptor_starting-20190929.png |}}
 +
 +
 +按下 **Send To** 按钮 {{:atk:sendto.png?direct&20|}} 将脚本发送到 **Editor** {{:atk:sendto.png?direct&20|}}。现在可以设置扭转参数了，以下部分定义了扭转势：
 +
 +<code python>  
 +potentialSet.addPotential(_potential)
 +_potential = CosineTorsionPotential(
 +   particleType1 = ParticleIdentifier('H', ['HC']),
 +   particleType2 = ParticleIdentifier('C', ['CT']),
 +   particleType3 = ParticleIdentifier('O', ['OH']),
 +   particleType4 = ParticleIdentifier('H', ['HO']),
 +   k = [0.0*eV, ],
 +   n = [0, ],
 +   delta = [0.0*Radians, ],
 +)
 +</code>
 +
 +使用以下代码块替换上面的部分
 +
 +<code python>  
 +potentialSet.addPotential(_potential)
 +_potential = CosineTorsionPotential(
 +   particleType1 = ParticleIdentifier('H', ['HC']),
 +   particleType2 = ParticleIdentifier('C', ['CT']),
 +   particleType3 = ParticleIdentifier('O', ['OH']),
 +   particleType4 = ParticleIdentifier('H', ['HO']),
 +   k = [0.176*kiloCaloriePerMol, ],
 +   n = [3, ],
 +   delta = [0.0*Radians, ],
 +)
 +</code>
 +
 +保存编辑过的脚本。
  
  
 == 定义 hook 函数 == == 定义 hook 函数 ==
 +
 +为了计算粘度，需要得到模拟期间压力张量的值。由于模拟过程中压力变化迅速，因此数据的高精度是必要的。虽然将结构保存至轨迹时也保存了压力张量，但是使用短的时间步长保存轨迹文件会导致文件太大超出合理范围。我们所需要的是采用一种方法只保存每一步的压力信息。
 +
 +在 **QuantumATK** 中，上述问题可以通过 hook 函数轻松解决。这些函数会在每次分子动力学整合之前或之后运行。hook 函数被设计用于控制模拟，或是在这种情况下，记录和分析模拟器件的不同量值。
 +
 +<WRAP center important 100%>
 +=== 注意 ===
 +有关如何使用 hook 函数修改模拟的示例，请参阅教程：[[https://docs.quantumatk.com/tutorials/youngs_modulus/youngs_modulus.html#youngs-modulus|Young’s modulus of a CNT with a defect]]。
 +</WRAP>
 +
 +为使用 hook 函数，必须首先添加一个实现 ‘__call__’ 函数的 Python 对象。传递给这个函数的参数是：
 +
 +  * //step//：调用函数时模拟的步数编号
 +  * //time//：到目前为止模拟中已经累积的时间量
 +  * //configuration//：包含当前结构的 //BulkConfiguration// 对象
 +  * //force//：每个原子上的一系列力
 +  * //stress//：当前结构的应力张量
 +
 +利用这些数值可以创建许多自定义分析和数据记录方法。在当前的案例中，需要记录压力张量，这可以通过下面脚本中定义的 //PressureRecordingUtility// 类来完成。
 +
 +<code python>  
 +# Define the class that does the analysis for the MD simulation
 +class PressureRecordingUtility:
 +    ''' Calculate and record the volume and pressure tensor at each step '''
 +
 +    def __init__(self):
 +        ''' Initialse the arrays used to store the data '''
 +        self._time_list = []
 +        self._volume_list = []
 +        self._pressure_list = []
 +
 +    def __call__(self, step, time, configuration, forces, stress):
 +        ''' Function called at each step that calculates and records the pressure tensor '''
 +
 +        # 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,3,1) * velocity.reshape(-1,1,3)) * masses.reshape(-1,1,1)).sum(axis=0)
 +
 +        # Construct the pressure tensor
 +        pressure_tensor = (velocity_tensor/volume) - stress
 +
 +        # Add the pressure tensor to the list
 +        self._pressure_list.append(pressure_tensor)
 +
 +    def times(self):
 +        ''' Function that returns the list of times saved during the simulation '''
 +        return Units.PhysicalQuantity(self._time_list)
 +
 +    def volumes(self):
 +        ''' Function that returns the list of volumes saved during the simulation '''
 +        return Units.PhysicalQuantity(self._volume_list)
 +
 +    def pressure_tensors(self):
 +        ''' Function that returns the list of pressure tensors saved during the simulation '''
 +        return Units.PhysicalQuantity(self._pressure_list)
 +</code>
 +
 +将此代码粘贴到模拟脚本的顶部，然后保存编辑的脚本。
 +
 +
 +
 +
 +
  
 == 在分子动力学模拟中使用 hook 函数 == == 在分子动力学模拟中使用 hook 函数 ==
 +
 +现在已定义了适合的 hook 函数，接下来需要告诉 **QuantumATK** 如何在模拟过程中使用这个函数。
 +
 +找到脚本中执行分子动力学计算的最终生成阶段的部分。在该部分的顶端添加以下行：
 +
 +<code python>  
 +pressureRecorder = PressureRecordingUtility()
 +</code>
 +
 +该行代码声明了用于记录压力数据的类的一个实例，分配用于存储压力数据的存储空间。
 +
 +现在需要更改 **MolecularDynamics** 的调用以添加 hook 函数。编辑此函数调用以包含参数：
 +
 +<code python>  
 +post_step_hook=pressureRecorder
 +</code>
 +
 +这行代码告诉 **QuantumATK** 使用 //PressureRecordingUtility// 的 //pressureRecorder// 实例运行和存储每一步的数据。
 +
 +
 +<WRAP center important 100%>
 +== 注意 ==
 +**MolecularDynamics** 函数支持参数 //pre_step_hook// 和 //post_step_hook// 分别用于分子动力学整合步数之前或之后运行 hook 函数。
 +</WRAP>
 +
 +最后，需要将 //pressureRecorder// 数据块收集的数据写出到文件中供后续处理，可以通过将数据写入 hdf5 文件来完成。为此，请将以下代码添加到脚本的末尾。
 +
 +<code python>  
 +data_file = 'md_methanol_data.hdf5'
 +temperature = 300
 +
 +nlsave(data_file, pressureRecorder.times(), "Time")
 +nlsave(data_file, pressureRecorder.pressure_tensors(), "Pressure")
 +nlsave(data_file, pressureRecorder.volumes(), "Volume")
 +nlsave(data_file, pressureRecorder.volumes().mean(), "Volume_Average")
 +nlsave(data_file, temperature, "Temperature")
 +nlsave(data_file, md_trajectory.timeStep(), "Time_Step")
 +</code>
 +
 +保存编辑后的脚本。
 +
 +
 +
 +
 +
 +
  
 === 运行计算 === === 运行计算 ===
 +
 +脚本中的所有必要部分已准备就位。点击 {{:atk:sendto.png?direct&20|}} 按钮将脚本发送到 **Job Manager**。
 +
  
 ==== 结果分析 ==== ==== 结果分析 ====
 +
 === 观察压力变化 === === 观察压力变化 ===
 +
 +采用 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** 平台的一些强大功能和灵活性，不止运行，还可以分析模拟。
  
 ==== 扩展结果 ==== ==== 扩展结果 ====
行 104: 行 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/用分子动力学方法模拟液体的粘度.1569762368.txt.gz · 最后更改: 2019/09/29 21:06 由 xie.congwei

页面工具

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