粗粒化晶格动力学计算与独立刚度近似方法的结合：多态分子晶体的比较研究

《Physical Chemistry Chemical Physics》：Coarse-grained lattice dynamics calculations combined with independent stiffness approximation: a comparative study on polymorphic molecular crystals

　　我们应用了一种先前报道的基于分子振动分析的晶格动力学计算方案，该方案结合了粗粒化（CG）方法和独立刚度近似（ISA）方法，以p-二氯苯晶体的多晶型作为测试案例。为了将这一理论框架的适用性扩展到复杂的晶体结构，我们利用了晶格内的分子对称性，从而使ISA能够有效地将计算成本降低到可管理的水平。我们仔细研究了在CG空间截断维度下，声子模式之间的耦合对它们频率准确性的影响。获得了三种多晶型的声子能带结构后，我们计算了它们的热力学函数，这些函数显示了由晶体结构引起的微妙但有意义的差异。

分子晶体的晶格动力学在最近的晶体工程研究中变得越来越重要。随着密度泛函理论（DFT）的广泛采用，集体晶格振动（声子）的理论计算变得更加容易。目前，主流的计算方法包括有限位移超胞方法（冻结声子方法）和密度泛函微扰理论（DFPT）。它们之间的一个根本区别在于初始计算过程中获得的力常数的形式不同：实空间中的力常数矩阵K与倒易空间中特定波矢下的动态矩阵D。这些矩阵的大小对应于原子振动的自由度，赋予了这些方法高精度和广泛的适用性。然而，对于单元晶胞内包含许多原子的分子晶体，原子级计算显著增加了计算成本。此外，这些方法输出的结果是单个原子的运动，需要事后分析将原子位移与分子作为一个整体刚体的振动行为联系起来。因此，为了将这种方法应用于分子晶体，已经做出了多种尝试来降低计算成本。

为了解决这些问题，我们开发了一种减少正常模式振动原子位移向量维度的方法，称为粗粒化（CG）方案。在这种方法中，通过使用分子组装的正常模式坐标变换特征值矩阵，可以很容易地获得描述与分子刚体运动和振动相关的刚度矩阵Φ，而无需显式推导力常数矩阵K。与有限位移方法类似，我们需要覆盖所需范围的分子簇，这样就可以使用DFPT方法中使用的解析微分协议来加速计算。此外，我们最近提出了一种称为独立刚度近似（ISA）的方法，它将整个问题划分为多个部分，从而将多体计算问题转化为多组两体计算问题。使用CG-ISA方案，我们之前已经完成了对简单分子晶体（p-二氯苯的β相）的声子计算。通过系统的理论验证和其他计算结果的比较，CG-ISA方案对于组成分子可以通过平移对称性关联的简单系统是成功的。然而，当每个单元晶胞包含两个或更多分子（Z）时，CG-ISA方法还面临着分子振动在更复杂对称关系中的挑战以及自由度的增加。这些挑战直接反映在Φ的结构上，而Φ与分子晶体的对称性密切相关。以下部分通过示例计算展示了CG-ISA方案对于Z > 1系统的新的理论扩展。在本文中，我们再次选择p-二氯苯作为案例研究，原因如下：（1）它是一个简单的小分子，但具有三种不同的多晶型，并且在相同温度下的晶体结构数据是可用的；（2）它基本上是范德华晶体，没有复杂的氢键网络；（3）有实验性的声子能带结构和拉曼数据。这些特点有助于我们理解分子晶体材料中的结构-性质关系，作为我们研究的初始阶段。

粗粒化方法的概述

晶格中粒子的振动可以通过求解质量加权Hessian矩阵的特征值方程来在谐波近似下理解，其中关键参数对应于惯性矩阵M和刚度矩阵K。在基于原子的分析中，这些矩阵的元素对应于原子振动的自由度。对于包含N个原子的分子组装，矩阵的维度是3N × 3N。通过引入一个粗粒化矩阵B，将原子位移向量投影到一个由表示刚体（三个平移和三个旋转）运动（Tx, Ty, Tz, Rx, Ry, Rz）和α个分子内振动模式（V1,…,Vα）的基组张成的缩减空间上，刚度矩阵被转换为对应于分子的振动模式，每个分子的大小为6 + α（方程（1）–（3）和图1）。我们将6 + α称为CG基的维度。

Γ?1 = BTMB

Φ = BTKB

B = (Tx, Ty, Tz, Rx, Ry, Rz, V1, …, Vα)

图1

粗粒化方案中刚度矩阵的变换。K的下标表示原子，而Φ的下标表示晶格节点。等价于运动方程的特征值方程可以使用矩阵U和Ω2来表示，它们分别收集特征向量u和ω2，如方程（4）所示：

Γ1/2ΦΓ1/2U = UΩ2

特征向量的元素直接对应于分子的振动模式，这与使用原子作为振动单位的传统方法不同。这种方法显著简化了分子振动模式的分析。通过从CG基中排除分子内振动模式，可以近似低频振动模式，从而有助于降低大型分子系统的计算成本。当α设置为0时，系统采用最小基，将分子近似为刚体，只考虑分子间的平移和旋转模式。当α设置为3N ? 6时，系统采用完整基，完全描述了对应于3N个自由度的振动模式。对于任意矩阵A，我们使用方程（5）定义其傅里叶变换的表示法[A]，其中{A}σ,σ′表示标记为σ和σ′的构成单元之间的非对角块。

在方程（5）中，aj表示作为分子σ的任意晶格节点Rσ的n1a1 + n2a2 + n3a3线性展开基础的j-th晶格向量。此外，P表示单个周期内包含的分子数量，k假设为当P → ∞时的连续倒易向量。使用这个表达式，我们在周期性边界条件（PBC）下制定了动态方程（方程（6）。

独立刚度近似

在基于CG方案的实际声子计算中，需要与中心晶格节点（编号为0）相关的刚度矩阵Φ的所有块的信息。通过块的傅里叶变换引入PBC以获得[Φ]（方程（7）。

为了完全计算一个周期内的分子数量，我们定义σ = 9n1 + 3n2 + n3，这是一个“n1n2n3”的十进制表示。在这个编号系统中，中心的1个晶格节点和周围的26个晶格节点（n1, n2, n3 = {?1, 0, 1}）用从?1到13的整数标记。

为了以降低的计算负担获得每个块{Φ}0,σ，我们引入了ISA的概念，其中作用在分子上的总分子间力可以计算为成对相互作用的总和，假设任意两个分子之间的刚度与其他所有分子无关。当每个节点中包含的分子数量（Z）为1时，可以通过选择通过平移操作相关的两个分子来计算非对角块{Φ}0,σ。当Z > 1时，单个节点包含通过反演、镜像反射、旋转、滑移和螺旋操作相关的多个分子。在6 + α维度的CG空间中，每个晶格节点包括Z × (6 + α)个振动模式，这意味着Φ的大小相应增加。以Z = 2为例，{Φ}σ,τ的组成如图2所示。在这里，单个节点内的不同分子被指定为A和B，它们可以通过给定空间群内允许的某种对称操作进行转换。通过组合字母A、B、C等来表示不同的操作和节点索引σ，可以唯一地识别每个节点内的不同分子。为了简单起见，我们引入了如{φ}0A,σB这样的表示法，表示{Φ}0,σ块的（A, B）子块（方程（8）。

{φ}0X,σY ≡ {{Φ}0,σ}X,Y?(X,Y = {A,B,C, …})

图2

在不同振动单位定义下刚度矩阵的结构：从包含两个分子的晶格节点（左）变换到包含一个分子的晶格节点（右）。对于Z > 1，有两种方法可以填充{Φ}0,σ。一种方法是通过将ISA簇中的振动单元范围扩展为单个单元来处理所有Z个分子（图2，左侧）。这种处理导致一个包含2Z个分子的“移动对”，简化了后续的理论问题，使其等同于Z = 1的系统。然而，随着Z的增加，正常模式计算的计算成本大致按（2Z）3 ～ (2Z)?的比例增加，而分子内振动模式的维度（α）至少需要6Z ? 6来包括Z个分子的相互位移。另一种方法是在本研究中采用的，即通过“镶嵌”每个节点内每两个分子对应的较小的{φ}0X,σY子块来构建{Φ}0,σ（图2，右侧）。随着Z的增加，子块的数量按Z2的比例增加，它们的排列变得越来越复杂。因为单个ISA簇计算只产生一个{φ}0X,σY子块，所以计算例程的总数与Z2成比例增加。在将CG-ISA方案应用于Z > 1的系统时，正确构建{Φ}0,σ的同时最小化计算例程的数量是核心理论挑战。在CG-ISA方案中，每个离散分子被分配给一个振动单元σX，{Φ}0,σ块中的数值关系通过设置适当的坐标系完全反映了给定分子晶体的对称性。如之前对Z = 1系统的研究所述，{Φ}0,σ的转置与{Φ}0,?σ相同：这种关系将非对角块的数量减少了一半。对于近似球形的分子，一个分子通常有12–18个相邻分子对其结合能有显著贡献，这意味着只需要计算6–9组分子对。类似的关系（方程（9）也适用于对应于任意分子对的子块。

对于处于平移关系（X = Y）的分子，方程（9）中的关系可以将计算所需的子块数量减少一半。对于处于非平移关系（X ≠ Y）的分子，适当的{φ}0X,σY子块的坐标变换可以提供另一个对应于共享等效几何关系的分子对的子块。这种变换关系已经通过理论制定并验证，其细节将在未来的研究中呈现。这种关系将独立子群的数量减少到1/h，其中h是空间群的阶数。为了充分利用方程（9）来降低计算成本，需要指定处于与给定分子对（0X, σY）等效的几何关系中的分子对（0X′, σ′Y′）。例如，对于p-二氯苯的γ相的计算，图3中显示的六个ISA簇是必不可少的，其中分子0A和0B可以通过21螺旋操作（沿b轴）进行转换。对于分子对（0A, 3A），即σ = 3（a = 0, b = 1, c = 0）且X = Y = A（相同操作），C2操作仅反转a轴和c轴的方向，导致σ′值为3（因为a和c都为零）。因此，我们发现分子对（0B, 3B）的几何关系与（0A, 3A）等效，而{φ}0B,3B是通过使用P[C2]（一个绕b轴旋转180度的矩阵）对{φ}0A,3A进行坐标变换得到的（公式（10）。

{φ}0B,3B = P[C2]{φ}0A,3AP[C2]T

（10）

图3

ISA簇用于分子对，由0X, σY组合表示。通过一系列这样的变换关系，对单个分子施加的力场就足以对晶格节点中所有Z分子所需的子块进行对齐。因此，Z值的增加并不会增加计算所需的ISA簇的数量。如表1所示，通过计算六个ISA簇，最终获得了30个{φ}0X,σY子块的矩阵元素。在这个表中，移动对（0X, σY）对应于图3中显示的ISA簇。通过计算粗粒度刚度矩阵，获得了“等效”分子对的子块{φ}0X,σY。这些子块的转置对应于“转置”列中的分子对。通过P[C2]的变换得到“C2旋转”列中的对，它们的转置对应于最后一列。随着系统中对称操作数量的增加，坐标变换步骤变得更加复杂。更多细节超出了本研究的范围。因此，CG-ISA框架内的变换将在未来的出版物中讨论。表1

列出p-二氯苯γ形式的等效分子对（0X, σY），其子块{φ}0X, σY是从包含指定移动对的ISA簇中计算得出的。有关其他列中的分子对，请参阅正文。

移动对（对称操作）
等效
转置
C2旋转
C2旋转 + 转置

0A, 0A（相同）
0A, 0A
—
0B, 0B
—
0A, 3A（b平移）
0A, 3A
0A, ?3A
0B, 3B
0B, ?3B
0A, 7A（a ? b + c平移）
0A, 7A
0A, ?7A
0B, ?13B
0B, 13B
0A, 10A（a + c平移）
0A, 10A
0A, ?10A
0B, ?10B
0B, 10B
0A, 0B（21螺旋）
0A, 0B
0B, 0A
0B, 3A
0A, ?3B
0B, 1A
0A, ?1B
0A, ?4B
0B, 4A
0A, 9B（21螺旋 + a平移）
0A, 9B
0B, ?9A
0B, ?6A
0A, 6B
0B, 10A
0A, ?10B
0A, ?13B
0B, 13A

如公式（8）或（9）所示，我们的方法在Z数值上没有任何限制，这意味着该方案适用于具有任意空间群的各种晶体。然而，在实际计算这种晶体时，我们必须提供考虑的空间群所需的完整P矩阵列表。同时，应该注意的是，对于具有两个或更多独立单元的晶体系统来说，情况并不简单。在这种情况下，我们必须分别处理这些单元，这意味着可能的ISA簇数量会迅速增加，簇的选择也变得更加复杂。目前，我们还没有找到一种有效的方法来选择对声子带计算至关重要的移动对。这更多是一个数学问题，而不是化学问题。

整个过程的方案在图S1（补充信息）中进行了说明。我们选择了p-二氯苯的α相和γ相作为计算目标（图4）。α相，CCDC参考代码DCLBEN07，单斜晶系，P21/a，Z = 2，a = 14.664 ?，b = 5.74 ?，c = 3.925 ?，β = 111.77°。γ相，CCDC参考代码DCLBEN03，单斜晶系，P21/c，Z = 2，a = 8.624 ?，b = 6.021 ?，c = 7.414 ?，β = 127.51°。对于这两种情况，对称中心都位于分子的中心；因此，只有一半的分子是独立的。为了计算方便，我们将晶体系统重新定义为P21空间群，使得一个晶格节点包含两个分子。通过这种处理，只有“平移”和“21螺旋”作为晶体系统的对称操作。否则，“反演”操作会使一个分子与自身重叠，从而失去分子方向的唯一性，进而失去CG矩阵的唯一性。尽管对称性有所降低，但仍然可以利用原始的对称性来减少独立分子对的数量。例如，由于分子的中心对称性，（0A, 0B）的几何关系与（0B, 1A）是相同的，这在P21空间群中需要区分（图S2）。我们确认{φ}0A,0B和{φ}0B,1A在可忽略的误差范围内是等价的，这可能归因于截断的ISA簇中的非严格对称性。表1列出了“等效”列中的分子对。

p-二氯苯多晶型的晶胞结构：（a）α相（参考代码：DCLBEN07）和（b）γ相（DCLBEN03）。在构建ISA簇时，应根据分子的重心距离（dMole）、原子间接近度（dAtom）以及{φ}子块的迹（作为分子间相互作用的度量）来选择必要的移动对。对于β-DCB，忽略了dAtom > 5 ?的分子，这种选择对声子带轮廓几乎没有影响。在当前两个系统中，由于存在相似的分子间相互作用，我们使用相同的标准来确定ISA簇中的移动对：γ-DCB的0A, 3A, 7A, 10A, 0B, 9B（表1）和α-DCB的0A, 1A, 3A, 4A, 0B, 1B。构建的ISA簇的结构经过迭代优化，直到在给定的分子环境中达到显著的最小势能。计算使用Gaussian16程序在B3LYP-D3/6-31G(d)理论水平上进行，这是一种结合了经验色散校正的混合DFT。同样的理论水平也应用于ISA簇的简正模式计算，在该计算中，移动对被设置为自由移动，其周围的分子被冻结。粗粒度刚度常数是通过使用正常坐标数据的适当过程获得的，如我们之前的研究中所描述的。系统的Z值增加并不一定意味着ISA簇的大小必须扩大。在这项研究中，我们使用的簇大小与我们之前对DCLBEN06的研究中使用的簇大小相当（Z = 1）。由于簇大小的截断而产生的误差通过声学求和规则（ASR）校正得到了补偿，重新定义了对角块中与平移相关的部分，以保持空间的各向同性和均匀性。图5显示了DCLBEN03（γ相）和DCLBEN07（α相）的第一布里渊区（FBZ），以及连接高对称点的k路径。从倒易空间的特性可以看出，即使两个晶体系统具有相同的空间群，FBZ的形状也可能不同。FBZ形状中k路径的选择决定了水平轴。因此，它对声子带的外观有主要影响，但不会改变态密度（DOS）和后续计算。

CG空间维度的影响

刚度矩阵的{φ}0X,σY子块是通过在B3LYP-D3/6-31G(d)理论水平上进行的简正模式分析计算得出的，我们的先前研究已经证明这种方法能够为β相晶体提供足够可靠的结果。然后，根据晶体对称性构建{Φ}0,σ块，并在假设PBC的情况下使用傅里叶卷积将其转换为倒易空间表示[Φ]（见理论部分）。图S3显示了在CG-ISA框架内获得的声子带低频区域（黑线）。图6将这个结果与非CG（全基组；α = 3N ? 6 = 30）计算（蓝线）进行了比较。全基组结果的全范围表示在补充信息中提供（图S4）。两组结果之间的相似性证明了将CG-ISA方法应用于Z = 2和Z = 1系统的计算的可行性。即，{φ}0X,σY的镶嵌和ASR校正对于单个晶格节点包含两个分子的情况是有效的。请注意，由于晶格内分子数量的增加，声子分支的数量增加到了12到140 cm?1，与β-DCB相比。图6

CG空间最小基组（α = 0，黑色）和全基组（α = 3N ? 6 = 30，蓝色）之间的声子带和DOS比较。（a）γ相和（b）α相。仔细检查γ-DCB的结果（图6(a)）发现，在#11和#12分支（约120–140 cm?1）附近，最小基组和全基组结果之间存在显著差异，而α-DCB（图6(b)）的频率分布几乎相同。为了评估CG空间对于单个ISA簇中移动对的质量，我们使用了保真度指数（FI），这是一个衡量原始和质量加权位移矩阵的Frobenius内积的指标。根据经验，当FI为0.8或更高时，CG空间的质量是令人满意的。对于α-DCB和γ-DCB，最小基组条件下的平均FI值分别为0.93和0.98，都达到了足够的水平。观察到的差异表明除了Φ的可重复性之外，还需要进一步考虑。当与ISA结合使用时，CG计算必然涉及具有不同FI值的多个簇，这意味着不能仅基于这个参数来确定近似质量。换句话说，适用于一个ISA簇的CG空间（选择有限数量的分子内振动模式）可能并不总是适用于所有其他簇。之前对β-DCB晶体的研究显示，最小基组和全基组结果在低频声子带和DOS方面有非常好的一致性。这个晶体系统可能是一个理想的情况，其中组成分子足够刚性，可以将刚体运动和分子内振动模式之间的频率区域分开。为了确定影响CG-ISA方案近似质量的因素，我们检查了低频分子内模式的影响。图7显示了使用全基组条件计算的#1–#16模式的声子带。对于γ-DCB，模式#11–14共享相同的频率区域，而对于α-DCB，这些模式被很好地分为{#11, #12}和{#13, #14}两组。在量子力学理论中，具有相同对称性的微观状态会相互作用，并由于避免交叉（AC）规则而改变彼此的频率，这伴随着不同状态的混合。如果CG空间中选定的基组缺少一些与AC效应相关的模式，模式混合将不会影响声子带，从而导致与全基组结果的显著偏差。通过使用α = 1计算声子带和DOS（图8）验证了这一假设，我们可以看到与全基组条件下的结果非常吻合。这个结果表明，在保持足够高的FI值的同时，包括共享相同频率区域的分子内振动模式是必要的。

全基组下模式#1至#16的声子带：（a）γ相和（b）α相。模式#11–14分别用紫色、黄色、红色和灰色表示。图8

CG基组α = 1（绿色）和全基组（蓝色）下γ相的声子带和DOS比较。Γ点频率的比较

全基组CG结果是“细粒度”极限，因此不会降低计算精度。因此，比较降维和全维CG空间可以评估CG方案的近似质量。与此同时，构建ISA簇的任意性影响了计算精度。在ISA方法中，错误主要源于簇之间的分子环境不平等，这是由于分子放置范围的截断以及结构优化导致的对称性微妙扰动。在我们当前的计算协议中，引入ASR校正可以补偿上述原因导致的大部分错误。然而，对于需要更多不同类型ISA簇的具有更复杂空间群的系统，可能会出现其他意外的错误。这些问题将在其他地方详细讨论。目前，我们仅将中等大小簇的计算结果与其他方法和实验值进行了比较。鉴于文献中缺乏可靠的实验或理论声子能带图，我们通过比较Γ点的频率来验证我们的计算结果（表2）。为了与实验结果进行比较，我们使用了一个与计算水平相对应的缩放因子来校正频率。总体而言，我们观察到计算值和实验值之间有很好的一致性，但在某些模式中观察到了显著差异。具体来说，对于α-DCB，CG近似表现得很好，计算结果与实验结果之间的差异很小（α = 1或30）。其中，模式#6–9的频率与参考计算得到的模式#7–10的频率更为接近。这种“不对齐”可能是由于模式排序定义的不同，需要使用基于群论的不可约表示来标记模式。在γ-DCB中，CG-ISA方案表现出更好的一致性，尤其是对于模式#8，其频率低于100 cm?1，与参考计算相比有了显著改进。然而，对于模式#11和#12，计算值与实验数据不符。相反，模式#13和#14的计算频率分别为134和141 cm?1，与实验值非常吻合。表2

在Γ点的频率（cm?1）与实验数据以及使用另一种方法计算的值的比较

模式#
α形式
γ形式

CG-ISAa
Calcd.b
Obsd.c
CG-ISAa
Calcd.b
Obsd.c

α = 1
α = 30
α = 1
α = 30

1–3
0
0
—
—
0
0
—
—

4
27.3
25.0
35.8
27
48.5
48.8
44.6
—

5
37.0
37.2
29.5
33
51.0
50.3
44.8
—

6
52.3
53.1
45.8
46
52.3
52.2
52.5
52

7
60.9
60.3
53.2
—
64.1
63.1
65
67

8
66.2
65.8
57
59
73.8
73.5
55.2
76

9
69.4
68.1
67
66
83.3
83.1
81.2
86

10
70.9
71.3
66
67
90.8
91.3
102.4
—

11
102.6
101.1
108.5
109
112.9
112.7
114.1
133

a
模式按照Γ点频率的升序编号。
b
来自参考文献39和40的计算数据。
c
来自参考文献41的拉曼数据（T = 1.2K）。晶格结构的变化主要影响低频声子，即分子间振动模式。与低频带中显示的适度变化相比，高频声子分支在三种多晶型中几乎相同，表明p-二氯苯晶体中的分子内振动模式受晶格结构变化的影响很小（图9）。这种不敏感性归因于这些晶体中缺乏强分子间相互作用，如氢键，以及分子的小尺寸和高刚性。然而，仔细检查差异后，发现分子周围环境的一些共同效应。图9显示了多晶型的本征波数减去孤立p-二氯苯分子的本征波数的结果。差异范围从-20到+40 cm?1，且波数依赖的位移趋势在三种多晶型中大致相似。值得注意的是，在低频区域（高达300 cm?1），α和γ形式的图谱相似，而在中频区域（800–1000 cm?1），α和β形式的图谱相似。图9

在高频区域（表现为分子内模式）的Γ点频率比较。鉴于β相的FBZ大小是其他相的一半，每个频率出现两次。关于DOS和热化学的比较研究

目前，多晶型的比较研究在晶体工程中具有核心重要性。在晶体结构预测的背景下，已经做出了大量努力，以高精度计算晶体多晶型的相对稳定性，同时考虑了熵和焓的贡献。尽管如此，按照正确的稳定性顺序列出多晶型仍然很困难。接下来，我们评估我们的方法解释分子晶体多晶型相对稳定性的适用性。比较了三种相（α、β和γ-DCB）在低频区域的声子DOS函数（图10）。在35 cm?1以下，三个相的频率分布非常相似，这一范围主要由声学模式主导。鉴于缺乏强分子间相互作用，如氢键，这一范围内的振动模式几乎不受分子晶体结构变化的影响。图10

p-二氯苯三种多晶型在低频区域的声子DOS（表现为分子间模式）。在35到100 cm?1的范围内，振动模式主要是低频光学模式。对于Z = 1的β相，这一频率区域主要由分子摆动主导，而其他两个相则表现出分子内振动模式的混合，如平面外的C–Cl弯曲运动。对于α和γ-DCB，不同的DOS峰出现在声学和光学分支的边界附近。在这些频率下，具有相当低群速度的声子模式倾向于积累。相比之下，β相没有明显的DOS峰，表明分子排列更为松散，从而导致声子频率的分布更广。在100 cm?1以上，振动模式从晶格运动转变为分子内振动，其中α和γ-DCB的分子内分支显示出不同的DOS峰，这取决于它们各自的晶体结构。实验上，热稳定性顺序为β > α > γ：γ相在271 K时转变为α相，转变焓（ΔtrH）为1.256 kJ mol?1，然后在304 K时转变为β相，ΔtrH为0.214 kJ mol?1，最后在326 K时熔化。相应的，转变熵（ΔtrS）计算为γ到α的转变是4.6 J K?1 mol?1，α到β的转变是0.70 J K?1 mol?1。这项研究的主要目标之一是从热化学角度解释这种稳定性顺序。如上所述，只要使用了足够的维度，CG方案对频率分布的影响可以忽略不计，表明其对其他声子衍生属性（如热力学属性）的影响也可以忽略不计。图11展示了p-二氯苯三种相的热化学函数。正如分子系统的统计动力学所建立的，振动对恒容热容量有贡献，因此也对熵和亥姆霍兹/吉布斯自由能有贡献。尽管这三个相之间的这些函数只有轻微的差异，但我们发现低频声子对总函数的贡献产生了不可忽视的差异（图11(a)）。

使用全基组ISA条件计算的三种多晶型的热化学函数。(a) 恒容热容量，以及(b) 熵，其中插图用虚线显示了低频贡献。(c) 相对于γ形式的吉布斯自由能，其中α和β形式的焓分别假设为+0.376和+0.820 kJ mol?1。在约100 K时，熵值的大小顺序为β > α > γ（图11(b)），表明随着温度的升高，?TS项有助于稳定β相。然而，差异如此之小，我们无法使用实验获得的焓值（Δγ→αH = 1.256 kJ mol?1，Δα→βH = 0.214 kJ mol?1）绘制出合理的G–T曲线：ΔH的0.1 kJ mol?1误差会导致曲线截距处的温度偏差40–60 K。在未来的研究中，我们计划计算不同相之间的焓变化，例如，考虑组成分子的电子能量和分子间的相互作用能量，尽管已经发现计算需要0.01 kJ mol?1的精度来合理预测转变温度。相反，我们尝试获得理想的焓值，以得到实验确定的转变温度（图11(c)）。结果，Δγ→αH和Δα→βH分别估计为0.376和0.820 kJ mol?1，与观察到的值相当接近。两者的总和为1.196 kJ mol?1，与实验确定的值1.470 kJ mol?1相当接近。为了更严格的量化，我们应该考虑晶格的温度、压力或应力依赖的各向异性应变，这会引起声子能带图的变化，这些变化并不容易预测。严格来说，从单声子能带图计算出的热化学函数并不能完全再现它们的温度依赖性；然而，这个问题并不是CG-ISA方法所特有的。鉴于这些限制，当前的热化学结果将分析限制在定性范围内。此外，虽然声子之间的相互作用在各种属性（如热导率）中起着关键作用，但谐波近似提供了“静态”声子，这些声子缺乏声子间的相互作用，因此无法准确再现“动态”热现象。传统的解决方案需要将非谐性效应纳入DFT计算中。目前的CG-ISA方案是基于谐波近似建立的，修改理论以纳入非谐性效应仍面临许多挑战。

在这项研究中，我们改进了用于Z > 1的分子晶体系统晶格动力学的CG-ISA方案。通过分析刚度矩阵的子块分量{φ}0X,σY与晶体对称性之间的关系，提出了一个全面的理论来覆盖ISA簇的构建和Φ的镶嵌。目前的进展使得能够计算Z = 2的p-二氯苯的α和γ相晶体。我们通过比较CG空间不同尺寸的计算来评估CG方案的近似质量。通过选择足够数量的分子内振动模式并考虑避免交叉规则，成功再现了截断基组的结果。在Γ点计算的频率与实验值非常吻合，而基于群论的声子分支的分配对于与实验值的严格匹配是必要的。对三种多晶型的比较显示，声学和高频分子内模式受晶格结构的影响最小，而低频光学模式对晶格结构表现出明显的敏感性。不同晶体相中的相似振动模式由于结构变化而表现出频率偏移。热化学计算表明，三种多晶型仅表现出微妙的差异，但在定性上是有意义的，这鼓励我们通过考虑额外的效应来进行更准确的评估。作者贡献

所有作者都参与了研究的构思和设计。Yue Wang主要进行了计算、数据处理和数据分析。Hirohiko Houjou组织了整个研究项目。两人共同研究了研究的理论内容。Y.W.准备了初稿，H.H.审阅和编辑了手稿。所有作者都批准了手稿的最终版本。利益冲突

作者声明没有竞争性财务利益。数据可用性

本文的数据，包括高斯输入文件、输出文件和粗粒化分析文件，可在OSF上获取：https://osf.io/583ux/overview?view_only=fc04744fdbd74e62af4f60c2a8fab9dd。补充信息（SI）也可获取。详见DOI: https://doi.org/10.1039/d6cp00159a。致谢

这项工作部分得到了JST SPRING（项目编号JPMJSP2108）和东京大学的世界领先创新研究生材料研究、信息和技术计划（MERIT-WINGS）的支持。注释和参考文献