版本：P-2019.03

在本教程中，您将以液态甲醇为例学习如何采用分子动力学（MD）模拟计算不同液体的粘度。理解粘度对许多工业化学过程的设计至关重要，因为粘度描述了液体的流动。本教程演示了您将如何利用 QuantumATK 工具模拟简单液体的粘性。该方法也可以应用于许多更为复杂和混合的液体。

为了模拟液体甲醇的行为，我们采用经典的分子动力学生成分子的运动轨迹。这些轨迹提供了统计计算液体综合性能所需的必要信息。本教程还演示了如何设置和使用为这类分子问题创建的特定键合力场。

粘度是液体中的分子相互通过时摩擦力的量度。有许多方法可以估算这种摩擦^[1]，但常用的一种方法是 Einstein 关系式，其中粘度 $\eta$ 可以表示为^[2]：

$$\begin{split}\eta = \lim_{t \to \infty} \frac{V}{2tk_BT}\left<\left( \int_0^t \tau_{\alpha\beta}(t') dt' \right)^2 \right>\end{split}$$

这里积分的量是压力张量的非对角线元素。标注的 $\alpha$ 和 $\beta$ 指代这些元素。此处其他的量代表的是时间 $t$、温度 $T$、体积 $V$ 和玻尔兹曼常数 $k_B$。

另外一种计算粘度的方法是通过 Green-Kubo 关系式^[3]：

$$\begin{split}\eta = \frac{V}{k_BT} \int_0^\infty \left< \tau_{\alpha\beta}(t_0) \tau_{\alpha\beta}(t) \right> dt\end{split}$$

这里尖括号中的数量是压力张量的自相关函数。由于该自相关函数在长时间尺度下会衰减为零，因此积分收敛可以得出流体粘度。

这两个方程都可以产生相似的结果。对于具有大量数据点的轨迹，Green-Kubo 表达式可能在计算上消耗过大。计算自相关函数尺度 $O(N^2)$ 中 $N$ 就是数据点的数量。Einstein 关系式的尺度仅仅只是 $O(N)$。由于计算粘度需要查看大量数据点，因此将把重点放在 Einstein 关系式上。

压力张量本身是由原子运动和原子与模拟中的边界条件相互作用两部分组成。压力张量可以表示为：

$$\tau_{\alpha\beta} = \sum_i m_i v_{i,\alpha} v_{i,\beta} - \frac{dE}{d\varepsilon_{\alpha\beta}}$$

在这里，$m_i$ 是原子 $i$ 的质量，$v_{i,\alpha}$ 是原子 $i$ 在笛卡尔坐标 $\alpha$ 方向的粘度，$E$ 是系统的能量，随应力 $\varepsilon_{\alpha\beta}$ 变化会有差异。

经典力场的原理是用被设计为近似总能量的简单函数代替计算上消耗过大的电子结构计算。这种方法可以确保有效地模拟更大的系统，远远超出了电子结构计算所能实现的。该方法的主要挑战之一是为替换函数拟合参数，以便它们能够准确地预测能量。

一种方法是将分子分成许多化学类型。然后使用这些类型找到不同类型交互作用的特定参数。例如，乙烷中的 sp³ 碳将与苯环中的 sp² 碳具有不同的关联势。由于将根据类型的定义排除一些组合，这就限制了需要找到参数的数量。因为参数适用于特定的原子类型和函数，这些势的类型通常会更准确。相对地，这些势往往更难以设置，因为每个原子的化学类型需要被结构指定。这种类型的势也不能模拟反应，因为需要原子连接来确定化学类型。

在本教程中，您将使用 OPLS-AA 力场模拟甲醇^[4]。OPLS-AA 力场为这些类型使用了具有特定拟合势的原子类型。本教程将介绍如何使用这些类型的力场设置计算。

打开 Builder ，选择 Add From Database。将数据库中甲醇的分子结构添加到 Stash，确保通过点击 Databases Molecules 选择分子数据库。

现在需要为每个原子指定化学类型。OPLS 甲醇有 4 种类型：

• 'CT' 为四面体碳
• 'HC' 为附着在碳上的氢
• 'OH' 为乙醇氧
• 'HO' 为乙醇氢

选择附着在碳上的三个氢原子，通过 Selection Tools Tags 将其类型设置为 'HC'，并输入标签 'HC'。您现在应该可以使用标记选择这些原子。对其他三种原子类型重复此过程。

接下来，选择 Builders Packmol 创建液态甲醇的块体模型。将甲醇分子从 stash 拖入分子类型栏，选择 250 个分子，晶胞长度设置为 28Å。点击 Create。

将鼠标悬停在新结构中的原子上，检查标签是否反映了每个原子的类型。

如果您的原子标签指定没有问题，点击 Send To 按钮 将块体结构发送到 Scripter 。

首先，在 Scripter 里添加一个 ForceFieldCalculator 模块。在此模块中需要定义一个力场，因此需要将新的力场定义添加到脚本中。打开 ForceFieldCalculator 并选择 New，在打开的 Potential Editor 对话框输入新的力场项。

给新的势命名为 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 种原子类型的非键合参数如下：

从弹出的周期表中选择元素，然后从下拉菜单中选中标签。添加相应的非键合参数。例如，CT 型的输入信息如下所示。

键伸缩能 $E_{bond}$ 可以用简谐波弹簧方程近似表达为：

$$E_{bond} = \sum_{ij} K_{ij}\left( r_{ij} - r_0 \right)^2$$

甲醇中有三种不同的键类型，这些键的参数如下表所示：

为添加键合势，在 Potential Components 框中点击 Add Add potential component。从势类型的列表中选择 HarmonicBondPotential。选中一个势，会在右边显示它的相关参数。添加三个势并为每个键型输入参数。CT-OH 键型的输入信息如下所示。

键角弯折能 $E_{angle}$ 可以用另一个简谐波弹簧方程近似表达为：

$$E_{angle} = \sum_{ijk} K_{ijk}\left( \theta_{ijk} - \theta_0 \right)^2$$

甲醇中有三种不同的键角类型，这三种键角参数如下表所示：

为添加键角势，在 Potential Components 框中点击 Add Add potential component，然后从势类型列表中选择 HarmonicAnglePotential。选中一个势就会在右侧显示参数。添加三个势并输入每个键型的参数。CT - OH - HO 键角类型的输入如下所示。

在甲醇中，只有一种键扭转项。这将延伸到 HC - CT - OH - HO 类型。该项中的键角扭转能可以通过 Fourier 展开近似表达为：

$$E_{torsion} = K (1+\cos(3\phi))$$

OPLS-AA 势中的总扭转项实际上是更为复杂的 3 级傅立叶展开。这个势不能直接在界面上添加，但可以通过少量修改脚本而轻松设置。创建势的步骤，在 Potential Components 框中点击 Add Add potential component，从势类型列表中选择 CosineTorsionPotential。稍后我们将在编辑最终脚本时设置参数。

虽然已在粒子类型定义中为单个原子类型添加了参数，但还需要为类型组合添加项。原子类型 HO 不是范德瓦尔斯位点，因此不需要添加涉及该原子的配对项。其余 3 种类型可以给出 6 种唯一的组合。

对于 Potential Components 框中的每个组合，点击 Add 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 非键合势的例子如下所示。

在粒子类型的定义中将原子电荷添加到模型中，且仍然需要添加静电项的求和法。

在 Potential Components 框中点击 Add Add potential component，从势类型列表中选择 CoulombSPME。这就选择了静电项的 Smooth Particle Mesh Ewald 求和。为了正确估算长程静电作用，本方法采用了 Fourier 转换。

右侧的参数选择为：

• 设置 $r_{cut}$ 为 9 Å。
• 设置精确度为 0.001。

其他参数的默认值适用于本次计算。

已经完成添加这些势之后，点击 Potential Editor 和 ForceFieldCalculator 的 OK。

现在已经完成了力场的定义，可以将计算添加到脚本中。

第一步是几何优化以从初始构型中消除任何大的力。

在 Scripter 中添加 Optimization 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

这两个细节的设置都将用于本教程的分析部分。

现在基本的计算已经创建完成，在 Scripter 设置输出文件为 Methanol_Viscosity.hdf5。脚本生成后应该如下所示：

按下 Send To 按钮 将脚本发送到 Editor 。现在可以设置扭转参数了，以下部分定义了扭转势：

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, ],
)

使用以下代码块替换上面的部分

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, ],
)

保存编辑过的脚本。

为了计算粘度，需要得到模拟期间压力张量的值。由于模拟过程中压力变化迅速，因此数据的高精度是必要的。虽然将结构保存至轨迹时也保存了压力张量，但是使用短的时间步长保存轨迹文件会导致文件太大超出合理范围。我们所需要的是采用一种方法只保存每一步的压力信息。

在 QuantumATK 中，上述问题可以通过 hook 函数轻松解决。这些函数会在每次分子动力学整合之前或之后运行。hook 函数被设计用于控制模拟，或是在这种情况下，记录和分析模拟器件的不同量值。

为使用 hook 函数，必须首先添加一个实现 ‘call’ 函数的 Python 对象。传递给这个函数的参数是：

• step：调用函数时模拟的步数编号
• time：到目前为止模拟中已经累积的时间量
• configuration：包含当前结构的 BulkConfiguration 对象
• force：每个原子上的一系列力
• stress：当前结构的应力张量

利用这些数值可以创建许多自定义分析和数据记录方法。在当前的案例中，需要记录压力张量，这可以通过下面脚本中定义的 PressureRecordingUtility 类来完成。

# 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)

将此代码粘贴到模拟脚本的顶部，然后保存编辑的脚本。

现在已定义了适合的 hook 函数，接下来需要告诉 QuantumATK 如何在模拟过程中使用这个函数。

找到脚本中执行分子动力学计算的最终生成阶段的部分。在该部分的顶端添加以下行：

pressureRecorder = PressureRecordingUtility()

该行代码声明了用于记录压力数据的类的一个实例，分配用于存储压力数据的存储空间。

现在需要更改 MolecularDynamics 的调用以添加 hook 函数。编辑此函数调用以包含参数：

post_step_hook=pressureRecorder

这行代码告诉 QuantumATK 使用 PressureRecordingUtility 的 pressureRecorder 实例运行和存储每一步的数据。

最后，需要将 pressureRecorder 数据块收集的数据写出到文件中供后续处理，可以通过将数据写入 hdf5 文件来完成。为此，请将以下代码添加到脚本的末尾。

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")

保存编辑后的脚本。

脚本中的所有必要部分已准备就位。点击 按钮将脚本发送到 Job Manager。

采用 NPT 系综时，检验模拟过程中压力的控制方式非常重要。理想情况下，压力应在模拟过程中最终达到目标值并趋于平衡。如果这种情况没有发生，则可能需要更长的平衡或模拟时间。

在 NPT 系综中，默认情况下会将压力写入 HDF5 文件。为了将压力可视化，可以使用脚本从 HDF5 文件中读取数据并绘图。为此，请点击 打开 Editor，粘贴以下脚本。

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()

也可以直接在此处下载脚本 ↓ pressure_analysis.py

将脚本发送到 Job Manager 运行。脚本应该显示与以下相似的图。

请注意，在动力学模拟期间，压力摆动约 0.2 GPa。这些大的压力变化是分子动力学模拟的典型特征。重要的是，模拟运行时平均压力收敛到目标 0.1 MPa。

在分子动力学模拟期间使用 hook 函数收集的数据，可用于计算和可视化粘度的估值。

以上可以通过使用 ATKPython 脚本来完成。在 Editor 中打开一个新脚本并粘贴以下脚本。

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()

也可以直接在此处下载脚本 ↓ 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^[3]。尽管该示例给出了类似的结果，但在计算的粘度中仍存在合理的噪声。理想地，该粘度应该收敛到一个不同的值。这表明需要更多的压力取样来提高模拟的准确性。

采用 Green-Kubo 关系可以获得类似的粘度结果。您可以在此处下载概述如何实现此目标的脚本 ↓ green_kubo_analysis.py。

我们鼓励 Python 的高级用户查看和运行脚本以了解其工作原理。该脚本还演示了 ATKPython 平台的一些强大功能和灵活性，不止运行，还可以分析模拟。

在上一节中，我们了解到为了提高估算粘度的准确性，需要进行更多的模拟。

在 QuantumATK 中，通过将它们放在 Python 循环中，可以轻松地重复模拟。在甲醇的例子中，一种策略是循环模拟的平衡和产出阶段。随着速度在平衡阶段开始时的重新随机化，将会产生一系列独立的轨迹。由此，采用标准统计方法，平均粘度及其误差都可以轻易估算。

在 Python 中循环计算，可以将要重复的代码放在 for 循环中。为了实现粘度计算循环 30 次以产生不同的轨迹，需要将模拟代码放置在 for 循环中，例如：

for loop in range(30):
    # Molecular dynamics simulation commands
 
# Futher commands outside the loop

创建多个轨迹的脚本中唯一必要的变化是修改输出的完成方式。在 HDF5 中，使用相同文本描述符写入数据的文件会覆盖以前的数据，因此需要唯一的描述符。这可以通过将循环变量（在上面的例子中变量循环）包含在名称中实现。更改文件写入代码为：

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))

为读取每条轨迹的数据进行分析，需要以相同的方式利用 nlread 命令更改名称。

为了演示如何使用多次运行提高模拟的准确性，模拟重复了 30 次。从这些模拟结构得到的平均值为 0.41 cP，标准误差为 0.05 cP。下图显示了平均粘度的变化和标准误差与模拟次数的关系。

• [1] B. Hess: Determining the shear viscosity of model liquids from molecular dynamics simulations. 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. 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. 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. J. Am. Chem. Soc. 118, 11225 (1996) 118