费米孔

Coulomb孔

由于Coulomb孔的存在，电子之间存在动态相关，如何描述原子分子体系中电子之间的动态相关是量子化学的主要任务

玻色子的波函数是对称的，没有反对称性的要求因此可以靠得很近

Hartree-Fock方程与自洽场计算

HF方程

HF方程是量子化学中计算多电子体系近似波函数和能量的基本方法，通过“平均场”近似把复杂的多电子问题拆解成单电子问题。

多电子波函数用单电子轨道（分子轨道）的Slater行列式表示 f̂(x_i)φ_k(x₁) = ε_kφ_k(x₁) Fock算符 $$ \hat{f}(x_1)=\hat{h}(x_1)+\sum_{j=1}^{N/2}(2\hat{J}-\hat{K_j}) $$

• 每个电子在其他电子构成的“平均场”中运动，忽略瞬时电子相关性（这是HF的局限性）。

类比理解

想象教室里的学生（电子）在考试：

• HF方法：每个学生只关心自己的试卷（轨道），并假设其他学生的行为是固定的“平均干扰”（库仑和交换势）。
• 现实：学生之间其实会互相偷看（电子相关作用），但HF忽略了这点。

自洽场计算（SCF）

一句话总结

SCF是求解HF方程的迭代计算过程，通过不断更新“平均场”直到结果自洽。

关键步骤

1. 猜初始轨道（如用原子轨道线性组合）。

2. 构建Fock算符 （依赖当前轨道）。

3. 解HF方程 ，得到新轨道。

4. 检查收敛

   新轨道与旧轨道是否几乎相同？

   • 是 → 计算结束，输出能量和波函数。
   • 否 → 用新轨道回到步骤2，继续迭代。

为什么叫“自洽”？

最终得到的轨道生成的“平均场”，与计算这些轨道时用的场是一致的（自圆其说）。

1
2
3
4
5
6
7

[多电子薛定谔方程]  
   ↓ 太难解 → 近似  
[HF方程]（单电子近似 + 平均场）  
   ↓ 需要迭代求解  
[SCF计算]  
   ↓ 收敛后  
[分子轨道和能量] → 用于计算化学性质

组态相互作用（Configuration Interaction, CI）计算

CI是一种后Hartree-Fock方法，通过将多个电子组态（Slater行列式）线性组合，考虑电子相关能，弥补HF方法的不足。

核心思想

• HF的缺陷： HF单行列式近似忽略了电子瞬时排斥（动态相关）和激发态贡献（静态相关），导致能量偏高。

• CI的改进： 将波函数表示为多个激发组态的叠加： Ψ_CI = c₀Φ_HF + ∑_ic_iΦ_i^single + ∑_jc_jΦ_j^double + ⋯

  • (Φ_i^single)：单激发组态（一个电子从占据轨道跃迁到虚轨道）。
  • ( Φ_j^double )：双激发组态（两个电子同时跃迁），贡献最大。

CI等级

关键问题：

• 计算量爆炸：FCI的组态数随电子数和轨道数阶乘增长，仅适用于极小体系。
• 大小一致性：CISD不满足大小一致性（如两个远离的H原子能量≠2×单个H能量），需更高激发或耦合簇（CC）修正。

1
2
3
4
5
6
7

[Gauss基函数]
   ↓ 基组展开
[Hartree-Fock方程] → 单行列式近似
   ↓ 电子相关修正
[组态相互作用(CI)] → 用二次量子化算符组合激发组态
   ↓ 更高精度
[耦合簇(CC)、多参考方法]

理论计算中的各种方法

简单总结了一下Sob老师的博文，个人的总结仅供参考，详细信息请阅读原文：http://sobereva.com/680

任务类型

• 优化极小点
• 优化过渡态
• 产生反应路径
• 振动分析
• 分子动力学
• 构象/构型搜索
• 波函数分析

计算方法

凡是基于量子理论思想提出的，在实际数值求解的过程中一般都要涉及分子轨道，绝大多数计算程序中都是把分子轨道展开成基函数的线性组合来描述的

最常见的基函数一类是原子中心基函数（如Gaussian等大多数量子化学程序以及CP2K等部分第一性原理程序用的高斯型基函数），其中心一般位于原子核，还有一类常见的基函数是平面波基函数，是大多数计算周期性体系为主的第一性原理程序用的，它的分布覆盖整个被计算的晶胞

基组（basis set）是对于原子中心基函数而言的，例如6-31G*、def2-TZVP、cc-pVDZ等都是很常见的基组，它定义了实际计算时对各种元素原子具体用多少、什么参数的基函数。做HF、DFT、后HF、MCSCF、多参考等方法计算时都需要告诉计算程序用什么基组，基组质量越好，也就是越接近于所谓的完备基组极限，这些方法自身的精度发挥得就越充分，但代价就是耗时越高。一个好方法配一个烂基组，以及一个烂方法配一个好基组，结果都不理想，必须好方法配好基组才能得到较准确结果。

Comment and share

pi-pi作用的定义和说法比较乱，Sob老师认为：pi-pi作用是相距较近的两个片段上彼此朝向相对的pi电子之间的独特的色散吸引作用。

1. pi-pi作用可以是分子间的也可以是分子内的，满足以上定义即可
2. 诸如苯分子里面的pi电子之间的作用不叫pi-pi作用，那属于共享电子作用
3. 两套pi电子的分布必须近乎彼此相对才可能算pi-pi作用。而诸如两套pi电子近乎肩并肩挨着就不能算pi-pi作用，像是不能说环丁二烯里面两套近乎定域的pi电子之间是分子内pi-pi作用
4. “相距较近”一般是在4.0-4.5埃以内。也不是说更远距离就完全没有pi-pi作用了，只不过由于色散作用随作用距离r呈$-\dfrac{1}{r^6}$快速衰减行为，因此距离稍微一远pi-pi作用就非常弱了，就不太值得一提了。但相互作用的片段间距离也不能太近，若显著小于相接触的原子的范德华半径之和，则显著的位阻互斥作用会远大于pi-pi吸引作用，使得pi-pi作用也相对不值得一提
5. pi-pi作用最常出现在碳原子间，因为碳最容易带显著的pi电子。碳的Bondi和CSD范德华半径分别为1.70和1.77埃。如果没有特殊的因素影响相互作用距离的话，在平衡结构（势能面极小点结构）下，C-C间的pi-pi作用出现在3.4-3.6埃左右是最常见的
6. pi-pi作用属于弱相互作用和非共价相互作用范畴。虽然名义上算作弱相互作用，但实际强度可大可小，直接取决于作用面积，详见后文
7. 有一个词叫pi-pi堆积（pi-pi stacking），这个词往往和pi-pi作用混用。但这个词更适合用来形容由于pi-pi作用的吸引效果使得相互作用的两个部分发生紧密结合从而产生的彼此堆积的结构特征

比如下面这个图是晕苯二聚体的局部极小点构型 （最小点为垂直构型）——平行错位构型，有一部分pi电子相对，算作pi-pi作用

pi-pi作用是一种非共价相互作用，内在本质是来自于相距较近的pi电子之间的色散作用。pi-pi堆积的体系有一个普遍特点是在极小点结构下，pi-pi作用区域的原子间通常不是正好对着，而是相互错位。下图是晕苯二聚体的俯视图，其中一个晕苯用红色显示，可以很清楚看到错位特征，即一层的碳对着另一层的六元环中心。一层层堆叠的石墨中每一层之间也同样是这样错位的。这可能是因为错位的结构下的静电相互作用能更负（静电吸引作用更强），也可能是因为错位的结构下比严格面对面的结构下的原子间的Pauli互斥更低。

下图为二聚体极小点结构下的单体的静电势着色的分子范德华表面的叠加图（18碳环），若上层静电势为蓝色则下层为红色，这样构型下的静电吸引作用比一个个碳原子彼此精确对着的时候更强。

虽然pi-pi相互作用与静电有些关系，但不可过分夸大静电的作用，pi-pi作用绝不是静电相互作用，比如下面这张描述pi-pi作用的图就是不合理的。

pi-pi相互作用能量的深入剖析可以借助Sob老师的提出的sobEDAw能量分解方法，本人比较菜，还没有阅读过相关内容。

pi-pi堆积的结构在现实环境中一般很容易发生分子间相对滑移（除非有额外的位阻效应等阻碍），这是由于==滑移导致能量变化很小。==

pi-pi相互作用的强度与涉及pi-pi作用的原子数目显著相关，若射击pi-pi作用的原子数目很多，相互作用能量甚至能够达到接近化学键的强度。

如果不知道pi-pi作用在哪， 不方便定义片段，则应当用使用IRI方法图形化考察化学体系中的化学键和弱相互作用介绍的IRI方法，可以把体系中所有相互作用区域全都展现出来，也包括化学键作用区域

考察体系的弱相互作用

对相互作用的认识

化学体系中的相互作用主要分为化学键作用和弱相互作用（区别于物理学中的术语）两大类。

从相互作用能上看，前者一般强度比较强，共价键和离子键都属于此类，也可以叫强相互作用；而后者一般强度比较弱，比前者通常弱一个数量级。弱相互作用既可以是分子间的，也可以是分子内的。弱相互包括范德华作用、pi-pi堆积作用、氢键、二氢键、卤键，以及后来很多人炒作概念而提出来的碳/硫/磷/金/银/铜…键等等，还有个词叫做非共价相互作用（noncovalent interaction），这个词的范畴相当于==弱相互作用和离子键的并集==。

弱相互作用形式多样，但主要本质只有两个：

• 静电相互作用。可以起到互斥作用也可以起到吸引作用，看具体情况（体系，以及相对位置）。这里说的静电相互作用也把极化作用包含进去了。
  • 色散作用。它起到==吸引作用==。必须从量子力学角度才能予以解释。而从量化理论角度来讲，它对应于==电子的长程的库仑相关作用==。

一般强度的氢/二氢/卤/硫/磷键等等是以==静电吸引主导，色散作用为辅构成的==，这点通过能量分解也可以体现出来。而一般说的范德华相互作用，以及pi-pi堆积作用的本质都是色散作用。

注：

1. 还有个作用叫交换互斥，它起到互斥作用，不管算什么类型相互作用都要考虑这个，==只有在较近距离（比如小于两个原子的范德华半径和）的时候才开始凸显出来，且原子间距离越近此作用越大==。正是这个作用使得弱相互作用的势能曲线总是有极小点，而不会令原子间距离因为静电吸引和色散作用而无限减小。
2. “范德华作用”这个词描述的是电中性的两个原子之间的非共价相互作用，故本质包含了交换互斥和色散吸引作用两个部分。这个词可以囊括比如两个氩原子间的相互作用的全部内涵。这个词也可以用来描述无极性分子间的全部相互作用（但实际上，没有绝对严格意义的无极性分子，比如就连一般说的无极性分子氮气也有四极矩，可以认为是局部极性，因此氮气之间也有静电相互作用，讨论见http://sobereva.com/209。但本文就不咬文嚼字了）。至于一些初级的教材里说的范德华作用包括“诱导力、色散力、取向力”，这种歪曲的说法应该彻底从量化研究者的脑海中删除
3. “范德华作用”这个词描述的是电中性的两个原子之间的非共价相互作用，故本质包含了交换互斥和色散吸引作用两个部分。这个词可以囊括比如两个氩原子间的相互作用的全部内涵。这个词也可以用来描述无极性分子间的全部相互作用（但实际上，没有绝对严格意义的无极性分子，比如就连一般说的无极性分子氮气也有四极矩，可以认为是局部极性，因此氮气之间也有静电相互作用，讨论见http://sobereva.com/209。但本文就不咬文嚼字了）。至于一些初级的教材里说的范德华作用包括“诱导力、色散力、取向力”，这种歪曲的说法应该彻底从量化研究者的脑海中删除
4. 《透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍》（http://sobereva.com/513）。

选择合适的方法分析说相互作用

• IGMH：如果你的目就是分析特定片段间或片段内的相互作用，很明确地知道片段该怎么定义，那么IGMH绝对是最理想的选择。从上一节的大量例子可看出，IGMH对各类体系表现得都很理想，图像很美观。

• IRI：如果你想将体系中所有相互作用（不分类型和强弱）在一张图里同时展现出来，肯定要用IRI，另外，IRI还有个重要优点是可以很好地展现出化学反应过程中化学键和弱相互作用的变化和过渡，==还可以做成生动的动画予以动态展现==，这在IRI介绍文章里都给出了，一定记得看。此外，IRI还有个变体IRI-pi，对于专门研究pi电子作用极有用处。

• IGM：有了图像效果好得多且物理上更严格的IGMH后，IGM基本就不用再考虑了。IGM仍有点用处的地方也就是这几个：(1)分析巨大体系时出于耗时考虑不得不用IGM (2)仅仅是想很粗略、快速地预览一下片段间弱相互作用出现的区域 (3)有坐标文件，但由于特殊情况不方便得到用于IGMH分析的波函数文件。

• NCI：有了IGMH和IRI后，NCI就基本没用了。IRI比NCI能展现更丰富的信息，尤其是能同时展现出强相互作用。而IGMH由于能自定义片段，比NCI方便灵活太多了，避免了无关区域的等值面妨碍分析，而且IGMH不需要特别精细的格点就能得到较平滑的等值面，从这个角度来说IGMH比NCI耗时还更低。若无特殊理由就没必要再用NCI了。

D3校正本质上校正的是体系的势能面，会令所有依赖于势能面的问题都受到影响，自然也包括优化出的结构（势能面极小点的位置）、振动频率（取决于势能面极小点处Hessian矩阵）、IRC（质权坐标下的势能面上的能量极小路径）等。因为D3校正仅仅是基于几何结构计算校正能，与电子态无关，而且只改变能量而不影响体系的波函数，因此对这些问题，加不加D3结果都一样：计算某个结构下的gap、轨道、偶极矩、极化率、NMR、原子电荷、键级、静电势分布等等。因为这些问题不涉及体系的能量或者能量对核坐标的导数。D3也不会影响这些问题的计算结果：垂直激发能、理论电子光谱（不考虑振动耦合时）、垂直电离能、垂直电子亲合能、垂直单-三重态能量差等等，因为计算不同电子态时用的几何结构相同，把D3校正能精确抵消了。但是，D3可以间接地对上述问题产生影响。比如一个二聚体，用B3LYP和B3LYP-D3优化出来的结构可能相差很大，那显然最终计算的偶极矩、吸收光谱等性质也会有明显不同。算垂直吸收/发射能和振子强度、计算激发态与基态的密度差、计算NTO等等，前面说了，D3丝毫不可能影响结果，因为结构没变，而且D3不影响波函数。而对于激发态几何优化，D3显然会影响结果，而且如果优化基态时加了D3则优化激发态也得加，得保持统一。

==不知道该不该加就加==

Non-covalent interaction(NCI)

NCI方法，2010年有Weitao Yang等人提出，核心思想是使用reduced density gradient(RDG)这一标量函数的三维等值面图来描述出现弱相互作用的区域，RDG函数定义为： $$ RDG(\bold r)=\dfrac{1}{2(3\pi^2)^{1/3}}\dfrac{\lvert \nabla \rho (\bold r) \rvert}{\rho(\bold r)^{4/3}} $$ 其中ρ是电子密度，可有量子化学程序计算出，也可由高分辨晶体衍射确定

在原子核附近，低能级上的电子高速运动，电子密度大，周围的电子密度变化迅速（即$\lvert \nabla \rho(\bold r) \rvert$比较大），总体上RDG($\bold r$)函数值中等。在化学键和弱相互作用区域，RDG相对较小，而化学键区域周围的ρ明显大于弱相互作用区域。因此，如果只考虑电子密度较低的区域，具有较小等值的RDG等值面将能够显示发生弱相互作用的区域。在NCI(非共价相互作用)分析中，通常使用等值为0.4-0.6的RDG等值面来直观展示弱相互作用，在ρ大于某一阈值的区域，将RDG设为任意较大值以抑制其等值面的出现。该阈值通常选择为0.05原子单位(a.u.)，因为在弱相互作用区域，ρ很少高于此值。然而，具有部分共价性质的弱相互作用，如带电氢键H-bond，其在相互作用区域的ρ可能明显更大，因此阈值需要根据实际情况进行适当调整。

如上图所示，(a)是RDG平面图，可以看出，在氢键区域（圆环中心），RDG值确实较小；与(b)相比，(c)的三维RDG等值面函数图屏蔽了$\rho (\bold r)$大于0.05 a.u.的区域（即屏蔽了原子核和键的区域，将这部分区域的RDG值设得较大）

此外，NCI(非共价相互作用)方法还提出通过不同颜色在RDG(电子密度梯度)等值面上映射$\operatorname{sign}\left( {\lambda }_{2}\right) \rho(\bold r)$函数，以直观区分相互作用的类型。$\operatorname{sign}\left( {\lambda }_{2}\right) \rho(\bold r)$表示$\rho (\bold r)$与λ₂符号的乘积，其中λ₂代表$\rho (\bold r)$的Hessian矩阵的第二大特征值。用于在RDG等值面上映射$\operatorname{sign}\left( {\lambda }_{2}\right) \rho (\bold r)$的常用色标如下图所示。由于化学体系中相互作用的复杂性和多样性，图中的标注仅适用于一般情况。

范德华(vdW)相互作用区域的ρ通常非常小，因此图2中对应的sign (λ₂)ρ颜色为绿色。此外，π − π堆叠具有与范德华吸引部分相同的物理性质，即色散效应(dispersion effect)。因此，π − π堆叠区域的ρ也相当小。对应明显的空间位阻效应或吸引性弱相互作用的区域(例如，中等强度的H-键和X-键)具有较大的ρ(通常为0.02-0.05原子单位)，因此在图2中分别对应红色和蓝色区域。虽然离子键的强度通常比弱相互作用大一个数量级，==但它们相互作用区域内ρ的大小是相当的==。例如，Na-Cl键的BCP(临界点)处的ρ仅约为0.03原子单位。因此，RDG等值面不仅可以揭示弱相互作用，还能显示离子键相互作用。相反，具有显著共价性的相互作用，尤其是共价键，==其在相互作用区域的ρ通常明显大于0.05原子单位==。因此，它们在采用密度截断的标准NCI(非共价相互作用)图中不可见。显然，NCI本质上是一种旨在揭示各种非共价相互作用的方法

在绘制NCI(非共价相互作用)图之前，波函数分析代码需要计算大量在矩形区域内均匀分布的网格点的RDG和sign (λ₂)ρ值，区域应合理设置以确保所有感兴趣的RDG等值面都位于其中。可以绘制RDG与sign (λ₂)ρ的==散点图==，以更好地理解NCI图的本质，这类图在文献中常用NCI分析进行讨论。下图(b)显示了对应于图(a)网格数据的散点图。从图中可以看到四个明显的峰，每个峰由多个点组成，延伸到图底。通过比较颜色，可以识别等值面图与散点图之间的对应关系，两个图中的圈内数字标示了这种对应关系。每个峰的最低点基本对应电子密度梯度为零的位置，实际上就是AIM理论中的==CP(临界点)==。图3(b)中的粉色虚线对应RDG为0.5，显然它只与上述四个峰相交，这也是(a)中<! − −swig0 − − >  = 0.5等值面能够揭示所有弱相互作用的原因。将RDG的等值线降低到低于0.5的值也是可以的，但如果设置得太低，NCI图中的RDG等值面会变得过小，难以观察。绘制NCI图的RDG等值线也不宜过大，否则可能出现不理想的等值面，严重干扰对感兴趣的弱相互作用的视觉分析。例如，如果RDG的等值线设置大于约0.7，散点图中的粉色虚线将与一些没有明显意义的网格点相交，从而导致出现意外的等值面。

如上图所示，分别描述了：

• (a)：N-I键的颜色为蓝绿色，故卤键的强度大于van der Waals相互作用
• (b)：苯甲酸钠粒子中，钠离子和阳离子之间等值面为蓝色表明二者的相互作用为离子键
• (c)：π氢键（苯环和氟化氢分子之间），强度较弱
• (d)：离子键+氢键
• (e)：一氧化碳与18碳环间的范德华相互作用
• (f)：$\rm C_{17}H{36}$的分子内显著的范德华相互作用（色散吸引作用），这也是该构象为该分子的最佳构象的原因（对于线性烷烃系统，当链长足够长时，发夹构象的能量低于直线构象的能量）
• (g)：由一个σ键连接的壬并苯二聚体中的显著的色散吸引作用(π − π堆叠)和空间位阻斥力

如果用于NCI分析的几何结构来自X射线晶体衍射实验且分辨率令人满意，并且打算研究晶体环境中的弱相互作用，则可以避免对非氢原子进行优化，但氢原子的位置仍需通过几何优化进行精细调整，尤其是在关注氢键的情况下，因为氢原子的坐标通常无法通过X射线衍射准确确定。否则，在进行NCI分析之前通常需要进行结构优化，否则可能出现下面的错误

NCI分析对波函数的质量要求不高，且几何优化任务本身比较基础，通常的做法是==结合中等规模的基组使用DFT进行计算==再将收敛后的波函数文件载入Multiwfn进行分析。尤其需要注意使用的DFT泛函需要能合理描述色散作用。

• 孤立体系：B3LYP-D3(BJ)和$\omega\rm B97XD$
• 周期性体系：PBE-D3(BJ)和TPSS-D3(BJ)

在基组的选择上，三重ζ(基组(如def-TZVP)且不含昂贵的高角动量极化函数是不错的选择。6 − 311G * * 提升基组质量通常不会带来明显的结果改善。使用包含极化函数的较便宜的双重ζ基组(如6 − 31G * *和def2-SVP也是完全可以接受的。

对于DFT超大体系，也可以使用GFN-xTB方法，它可以被是为DFT方法的半经验变体，GFN-xTB能够轻松在个人计算机上优化包含数百甚至上千个原子的系统。

NCI图的图形质量对RDG(电子密度梯度)和sign (λ₂)ρ的网格数据质量非常敏感，这主要体现在网格点之间的间距上。网格间距过大会导致等值面出现难看的锯齿边缘；但如果网格间距过小，则会带来过高的计算成本。请注意，如果空间范围固定，计算网格数据所花费的时间与网格间距的三次方成反比。根据经验，网格间距为0.2 Bohr通常是NCI图的最低可接受质量；如果对图形质量要求较高，明显应使用小于0.15 Bohr的网格间距。在开始计算之前，应正确设置网格数据的空间范围。一般指导原则是空间范围应充分覆盖所有可能出现感兴趣的微弱相互作用等值面(isosurfaces)区域，否则可能会遗漏或截断某些等值面。例如，如果微弱相互作用分布在整个系统中，并且希望在NCI(非共价相互作用)图中同时显示所有这些相互作用，则网格数据的范围应覆盖整个系统。相反，如果感兴趣的微弱相互作用仅发生在大型系统的局部区域，则可以合理定义空间范围，使得网格点仅分布在感兴趣的区域，这不仅大大减少计算时间，还能避免在无关区域出现等值面。用于绘制RDG等值面的等值值(isovalue)也可以根据实际情况适当调整，以改善图形质量。大多数情况下，选择0.5作为RDG值是合理的。通常可以将其略微提高到0.6，这样等值面会变得更宽，更易于观察。然而，如果设置得过大，则会导致无意义且杂乱的等值面。

IRI

IRI(Interaction Region Indicator)，由卢天老师提出 $$ \rm IRI(r)=\dfrac{\lvert \nabla \rho(r) \rvert}{[\rho(r)]^a } $$

在常规的量子化学计算之后可得波函数信息，再进一步可以计算出电子密度ρ及其梯度，即可计算IRI，一般a=1.1，这是最优的经验数值。

Hessian矩阵是电子密度函数的二阶导数矩阵，描述电子密度在空间中的==曲率变化==。数学形式为： $$ H_{ij}=\dfrac{\partial ^2 \rho}{\partial x_i \partial y_j}\quad(i,j=x,y,z) $$

• 物理意义
  • 正曲率（凸起）→ 电子密度在某个方向“向外膨胀”（如原子核附近）
  • 负曲率（凹陷）→ 电子密度在某个方向“向内收缩”（如化学键或分子间区域）

λ₂是其Hessian矩阵的第二本征值，在AIM理论中：

• 键临界点处sign(λ₂)=-1
• 环、笼临界点处sign(λ₂)=+1
• 接近临界点的区域其值与临界点处一般相同

可以将sign(λ₂)函数用不同色彩投影到RDG等值面上，用来表现某一个区域的相互作用类型，==这便是用RDG方法分析弱相互作用==，但不如==IRI方法==全面，IRI与RDG在展现弱相互作用方面效果是等同的，而IRI有个关键性的好处是可以把化学键作用区域一起直观地展现出来！因此可以一目了然地图形化考察体系当中存在的所有类型相互作用，更有价值。考察体系的弱相互作用时会绘制电等值面图，IRI方法绘制的等值面图横坐标是sign(λ₂)ρ(r)，纵坐标是IRI函数；而RDG方法绘制的等值面图横坐标只有ρ(r)，纵坐标为RDG函数

RDG

RDG(Reduced Density Gradient)，即约化密度梯度，是量子化学中==一种比较落后==的可视化弱相互作用的重要函数 $$ \rm RDG=\dfrac{\lvert \nabla \rho \rvert}{2(3\pi^2)^{\frac{1}{3}}\rho^{\frac{4}{3}}} $$ ρ(r)只能反映出强度，但类型需要由sign(λ_2)函数来反映，这个函数是电子密度Hessian矩阵的第二大的本征值λ_2的符号，在AIM理论中键临界点的sign(λ_2)=-1，环、笼临界点的sign(λ_2)=+1，在接近临界点的区域其值与临界点处一般相同。可以将sign(λ_2)函数用不同色彩投影到RDG等值面上，用来表现某一个区域的相互作用类型。

Comment and share

从固体材料的化学组成来区，主要有无机、有机和金属三大类材料。固体物质可以按照其中原子排列的有序成都分类为晶态固体和非晶态固体，后者类似于液体，只在及个原子间距的短程范围内具有源自有序的状态，比如玻璃和许多聚合物可以看作是过冷的液体，它们中间的原子排列是没有一定的格子的。

金属键和金属晶体

金属的自由电子模型

20世纪初，Drude和Lorentz提出了自由电子论。金属都是电离能第、颠覆性小的元素，这些元素容易失去外层价电子，形成正离子，正离子按照紧密堆积的方式形成晶体，价电子在整个金属晶体中自由运动，为所有原子所共有。自由电子处于类似于理想气体的状态，可以看作是自由电子气

在金属表面存在着一种把自由电子限制在金属范围内的是能查，但是在金属内部，势能是相对均匀的，好像自由电子是在一个均匀的正电场中运动着，势能相对为0.

当沿着与电流在金属晶体中运动方向垂直的方向上再施加一个外磁场时，电子将受到Lorentz力而偏向某一侧，这便是霍尔效应。大多数金属的霍尔系数小于零，即载流子是电子。然而有一些金属如Zn，Cd，Pb等，其霍尔系数大于零，载流子密度月底，霍尔系数值越高，但是由霍尔效应测得的载流子密度并不和价电子密度相同，这是自由电子论所不能解释的。

实际上，电子是在晶体中所有格点上的离子和其他所有电子所产生的势场中运动，金属中的电子的势能是位置的函数，并非为常数，并且电子在金属中的运动遵循量子力学而非经典力学。

金属的近自由电子论

电子的波长： $$ \lambda=\frac{h}{mv} $$ 电子处于一个很深的势箱中运动，晶体的边缘长度L就是势箱的边界。对于一维势箱，电子的动能为： $$ E=\frac{1}{2}mv^2 $$ 对于三维势箱： $$ E=\dfrac{h^2}{8mL^2}(n_x^2+n_y^2+n_z^2)=\dfrac{h^2}{8\pi^2m}(k_x^2+k_y^2+k_z^2) $$ 其中k_x, k_y, k_z分别是波矢k在三个方向上的分量

Comment and share

==为什么方形电池设计为带正电的？==，作为一个小白的无意间看到了这么一个问题，仔细学习了一下，参考许多资料。

动力电池的电芯封装模式：

• 硬壳
  • 圆柱：技术成熟成本较低、稳定耐用、单体一致性好，但能量密度的上升空间小、对电池管理系统的要求高
    • 铝壳
    • 钢壳
  • ==方形==：强度高、内阻小、寿命长、空间利用率高，应用最广，但厂商需求不一样，规格多且难统一
    • 铝壳
    • 钢壳
• 软包：能量密度极高、重量轻，但成本高，需要额外防护防止电池受损和热失控，更容易因为膨胀产生内应力、膨胀力，导致电池内部的结构都会发生改变
  • 铝塑膜
  • 钢塑膜

标题中的方形电池应该指的是铝壳封装的方形电池，其中的电芯通常与正极相连，出于以下五个方面考虑

1. 容易形成Al-Li合金。Al的立方晶格的八面体空隙大小与$\rm{Li^+}$半径匹配，在负极电位下，二者极易形成金属间隙化合物，随着嵌锂的深入，氧化锂、氢氧化锂等副产物会逐步形成，腐蚀反应会使得铝外壳逐渐粉化
2. 铝的嵌锂电位$\varphi \rm(vs~Li^+/Li)$较高，大致0.3V，高于石墨负极的嵌锂电位$\varphi (\rm vs~Li^+/Li)$的0.01-0.2V，如果石墨和铝同时作为负极那么铝将优先发生嵌锂反应。为了避免铝的嵌锂反应，需要将铝外壳的电位提升至嵌锂电位以上，

EIS在电池表征中的应用

电极材料的表征

• 电子电导性评估
• 离子扩散特性研究

界面行为分析

• SEI膜：中频区半圆直径对应电荷转移电阻$\rm R_{ct}$，循环后$\rm R_{ct}$突增可能预示着SEI膜破裂再生（如六氟磷酸锂分解产物积累）
• 复合电极层间接触：多层电极（如NCM/石墨全电池）中，中高频区额外半圆可揭示正负极界面接触电阻

电池老化与失效诊断

• 活性锂损失：全频阻抗缓升，$\rm R_{ct}$与$\rm R_W$同步增加（如负极析锂导致可逆锂损失）

• 结构衰退：低频Warburg斜率陡降（如NCM正极晶界裂纹阻碍离子传输）

• 极端失效识别：

  • 微短路：高频截距$\rm R_s$异常降低（如隔膜穿透引发电子泄露）
  • 析锂：低频区“扩散尾”畸变（如锂枝晶导致传输异质性）

锂电池中需要铜箔吗？

铜箔在锂电池中作为负极的集流体，主要用于收集锂电池在放电过程中的电子并有效地导出，确保电池的电流传输效率。此外，铜箔能够支撑负极活性材料（如石墨、硅等），确保电池在充放电过程中结构稳定。

机械性能

• 厚度及均匀性：会导致电池性能不均匀。厚度偏差需控制在±0.5μm以内
• 面密度一致性：单位面积质量极差≤1.5g/m²
• 抗拉强度：高抗拉强度可提高电池生产中的涂布碾压效率，避免断带现象。普抗(300-400MPa)、中抗（400-500MPa）、高抗(500-600MPa)、超高抗（>600MPa)，5μm产品要求≥430MPa
• 延伸率：高延伸率可增强极片的压实密度并降低电极片的厚度，更好地抑制电池循环过程中活性材料的变形，从而提高电池的耐久性和使用寿命6μm常规延伸>4%，高延>6%;8μm常规>5%，高延>8%，5μm铜箔延伸率要求≥6%

物理化学性能

• 表面粗糙度：表面粗糙度对负极活性物质的附着和涂布效果至关重要，以提升涂覆均匀性，减少活性材料脱落。
• 粗糙度Rz≤2μm
• 纯度：避免杂质在电池充放电过程中因杂质的存在而引发副反应，影响电池的容量、寿命和安全性。要求≥99.999%（六个九），高端产品达99.9999%(七个九)，杂质(如Fe、S)需<0.001%
• 电导率/电阻率：降低电池内阻，可以减少电池内部的能量损失，提高电池的充放电效率和功率性能；纯铜的电阻率约为1.694$\rm \mu\Omega \cdot cm$，锂电池用铜箔的电阻率需尽可能接近纯铜理论值

表面特性

• 清洁度：无压痕、皱褶、缺口、污物及氧化缺陷
• 亲水性：表面需良好润湿性，增强与负极材料及粘结剂的结合力
• 抗氧化性：180℃高温环境60分钟无明显氧化
• 良好的耐腐蚀性：防止其被电解液等物质腐蚀，从而保证电池的性能和寿命

热稳定性和微观结构

• 高温屈服强度：需在150C制造过程中保持原始机械性能，防止软化
• 微观结构：铜晶粒的取向和织构需优化避免晶界过多导致电子传输受阻，从而维持低电阻率

电解液浸润不良可能导致电池析锂

极化效应

• 欧姆极化

孔隙填充率下降，浸润不良区域，电解液渗透率<80%时，离子传输路径的曲折度增加40%-60%；浸润不良导致电极与电解液 的有效接触面积减小，界面电阻增大

• 浓差极化

扩散延迟，孔隙内电解液缺失使锂离子扩散系数降低

表面浓度降低，充电过程中电极表面Li浓度可达本体浓度的2-3倍，触发局部析锂的临界电流密度降低40%

• 电化学极化

活性位点减少，由于电解液进入不足导致电极与电解液活性位点减少，以及有效嵌锂位点减少，电荷转移电阻增大

动力学-热力学耦合效应

浸润不良打破了锂离子嵌入/沉积的动态平衡，沉积电势窗口变窄，局部区域有效电势低于锂沉积阔值的概率提 升2-3倍

SEI膜形成异常

浸润不良时，反应物无法到达特定区域，导致SEI膜无法连续覆盖，不均匀的SEI膜中薄弱区域易因体积膨胀或机械应力破裂，会优先成为锂金属沉积位点，诱导锂枝晶的形成

电池传输路径受阻

电解液浸润不良，导致电极可能会出现部分干涸区域，锂离子传输需要绕过浸润不良区域，增加了锂离子的迁移路径，导致局部锂离子扩散速率降低，传输受阻

电流密度

浸润不良的区域因电解液缺乏，阻抗升高，电流被迫集中于少数通道，高电流密度区域锂离子嵌入速度达到极限，超出负极动力学承受能力，引发局部析锂

关于SEI膜

• LiF化学稳定性极高，电子绝缘优异，可以有效提升界面抑制副反应；但其锂离子扩散能垒非常高（约0.5eV以上），当其在界面中过度堆积时，会显著阻碍Li+穿透SEI并进入负极，从而对倍率性能、低温性能产生不利影响

• Li₂O易形成多孔、富缺陷的结构，其离子电导率高于LiF，能够在一定程度上缓解界面阻抗，特别适合快充或低温场景。然而其稳定性略低，循环过程中容易参与副反应并诱导SEI逐渐增厚或成分演化，对长寿命形成挑战

• 理想的SEI应由多相协同构成，其中LiF负责界面钝化与电子阻断，Li₂O则提升离子通量与机械缓冲性能。通过精准调控电解液配方、FEC浓度、锂盐类型与化成条件，可以构建出具有“分散型LiF+结构调和型Li₂O”的纳米复合SEI结构，真正满足未来快充与高能电池对界面性能的双重要求。

电池中使用的各种表征手段

电池材料结构表征

电池材料成分与化学态分析

电池材料电子结构与能带分析

电池材料界面与动态过程表征

电池材料热安全与机械性能

电池材料电化学性能表征

安全与失效分析

一些基础知识

C rate

即倍率，描述了电池充电或放电的速度。假定使用的电池的额定==容量==Q= 40 ==mAh==，倍率为0.2 ==C==。这意味着理想情况下，根据以下等式电池能以8 mA的恒定电流放电5小时 Q = It

$$ I=\dfrac{Q}{t}=\dfrac{Q}{\frac{1}{C}}=\dfrac{40}{\frac{1}{0.2}}\rm mA=8mA $$

实验数据如下表所示：

可以看到，随着倍率C的增加，放电时间减小，IR压降更多，因此容量也逐渐减小。

等效串联电阻ESR随着倍率增大而降低，这是电池内部温度升高造成的。然而，容量和能源较低等缺点超过了这一优势。此外，较高的温度也可能会导致材料变质

电池循环测试

测试电池长期稳定性的典型实验是循环。为此，电池会被充电和放电数百甚至更多次并测量容量。

下图显示了电池的标准循环充放电（CCD，即恒电流充放电）实验。首先以1.0 C速率（40 mA电流） 将纽扣电池充电至4.2 V。然后，该电位保持至少72小时，或者如果电流达到1mA。然后以1.0 C的速率将电池放电至2.7 V。该序列重复 100 个周期。较暗的曲线显示容量，较浅的曲线显示相对于开始时的容量百分比。

电解质杂质或电极缺陷总是会导致容量损失。 本例中测试的电池显示出良好的循环行为。 纽扣电池的最大容量约为 28.7 mAh。 100 次循环后容量仅略有下降。 总容量损失总和约为 4.5%。极端温度、过度充电或过度放电会加速容量损失。一般情况下，当容量损失超过20%时应更换电池。

库伦效率被定义为： $$ \rm \eta_C=\dfrac{Q_{discharge}}{Q_{charge}} $$ 库伦效率总是小于100%，这也是测试中放电时的容量小于充电时容量的原因，本次测试纽扣电池表现出约98%的库仑效率。

漏电和自放电

理想情况下，当没有外部电流流动时，电池的电位是恒定的。然而，实际上，即使电池没有连接到外部负载，电位也会随着时间的推移而降低，这种效应称为自放电（self discharge），所有储能设备或多或少都会受到自放电 （SD） 的影响。

下图显示了使用新型纽扣电池的自放电实验图。首先将电池充电至 4.2 V，然后在此电位下保持恒电位三天。然后测量开路电位九天。

最初，电位降低超过 6 mV。之后，速率减慢至每天不到 1 mV。9 天后，电位总共降低了 15.6 mV。这相当于相对于初始值仅下降 0.37%，具体数据如下表所示：

自放电是由内部电流引起的，称为漏电流（$\rm I_{leakage}$）。自放电率主要受电池的使用年限和使用情况、初始电位以及温度影响。

下图显示了两个纽扣电池的漏电流测量值。一块电池是新的，另一块电池短时间加热到 100°C。两块电池最初都充电至 4.2 V。然后保持电位恒定并测量电流。

许多制造商将$\rm I_{leakage}$ 指定为 72 小时后的测量值，但请注意，即使在72小时候测量电流依旧不断下降，不是恒定的。

在此情况下，新电池的漏电流约为4.7 μA，而老化的纽扣电池显示出两倍大的值，约10 μA。

Gamry公司建议不要使用恒电位测试来测量漏电流。

电化学阻抗谱EIS

关于EIS的基础知识，之前也是写过一篇博文

下图显示了不同电位下的四种不同的Nyquist图。纽扣电池首先分别充电至3.9 V、4.1 V、4.3 V和4.5 V。然后保持电位，直到电流降至1 mA以下。此步骤可确保电位在EIS实验期间保持恒定。恒电流EIS实验在100 KHz和10 MHz 之间进行。直流电流为零，交流电流设置为10 mArms

奈奎斯特图的形状在很大程度上取决于电池的电位。在较低电位下，即 3.9 V 和 4.1 V，两条曲线几乎重叠。电池的阻抗在更高的电位下增加。4.3 V 和 4.5 V 处的奈奎斯特曲线分别向右移动，半圆更大。下图显示了锂离子电池的典型模型

$\rm R_{ESR}$代表电池的ESR电阻。它是高频下的极限阻抗，可以很容易地估计为奈奎斯特曲线和x轴（$\rm Z_{real}$） 之间的交点。

此外，假设每个电极/电解质界面具有双层电容和电荷转移电阻$\rm R_{ct}$。这些元件的每个并联电路代表奈奎斯特图中的一个半圆。 为了解决两个电极的孔隙率和不均匀性问题，双层电容被恒相元件（CPE） 取代，它总结了双层在非理想电极/电解质界面处的极化效应。==理想情况下，CPE 可以被视为电容器。==

==所有奈奎斯特曲线在低频下都显示一条-45°角的对角线==。该区域可以通过Warburg阻抗$\rm Z_W$进行建模，它描述了扩散层无限厚度的线性扩散。为了简化问题，仅考虑一个电极上的扩散。

其中，参数Y及其无量纲指数a定义了CPE，Y的单位为$\rm S\cdot s^a$。

CPE定义为： $$ Z_{CPE}=\dfrac{1}{Y_0(j\omega)^a} $$ Z_CPE的单位为Ω，ω^a的单位为$\rm s^{-a}$

• 当a = 1时：

$$ Z_{CPE}=\dfrac{1}{Y_0(j\omega)^1}=\dfrac{1}{j\omega Y_0} $$

与电容器的阻抗形式$Z_C=\dfrac{1}{j\omega C}$相同，此时CPE元件等价于理想电容器，Y₀的单位$\rm S\cdot s^1$等价于法拉第(F)。CPE元件代表一个完美的、没有能量损耗的电荷存储元件。在电化学中，这可能对应于一个完全平整、均匀且不发生任何化学反应的理想电极表面

• 当a = 0时：

$$ Z_{CPE}=\dfrac{1}{Y_0} $$

这是一个不随频率变化的常数，也就是一个纯电阻。此时，Y₀是这个电阻的电导G，而$\frac{1}{Y_0}$就是电阻值R，CPE元件代表一个纯粹的能量耗散元件，通常模拟电解液的电阻或电路中的其他欧姆电阻。

• 当a = 0.5时：CPE模拟的是一种特殊的电化学现象——==扩散过程==，这被称为瓦氏阻抗。这代表了电化学反应中，反应物从溶液主体传输到电极表面，或者反应产物从电极表面扩散离去的过程所遇到的阻碍。这种阻碍与频率的平方根成反比，是==受扩散控制==的电化学反应的典型特征。在Nyquist图上，它表现为一条与实轴成45°角的直线。

对称电池是啥？

具有相同电极的电池或对称电池用于==研究电池电极的行为、研究插入反应以及确定电池中正极和负极的阻抗==

对称电池由两个发生氧化还原反应的相同电极组成，电极间通过液态或固态离子导体隔离，电荷通过离子传输实现迁移。使用对称电池的主要优势在于可以研究两个相似的界面，而全电池甚至半电池中涉及的是两种不同界面。在电池系统中，电极是由多种化学物质组成的复合体，包括活性材料、粘结剂、电子导体等成分……对称电池使我们能够了解在特定电解质或特定电极条件下，这些化合物在氧化和还原过程中的相互作用及电化学稳定性。

让我们考虑一个由两个锂电极组成的对称电池，两电极间填充含有锂离子的电解质，该电池两端的电位差为零。若施加正向电流（ I>0），那么其中一个电极将成为阳极，其界面会发生如下氧化反应： $$ \rm Li \longrightarrow Li^+ +e^- $$ 另一电极则成为阴极，其界面会发生如下还原反应：

$$ \rm Li^+ +e^- \longrightarrow Li $$ 下图展示了对称锂电池在正负电流阶跃作用下的电位响应。τ₁被称为过渡时间常数，超过该值后”（在还原反应情况下）流向电极表面的О通量[…]不足以满足施加电流需求，电位将偏移至更极端的数值，==此时可能发生其他电极过程（比如下面所说的集流体的氧化）==

就I>0的正向电流而言，t = τ₁处电位的快速上升可能对应着锂电极的完全消耗，由下图所示：

上图中，阳极厚度在转变时刻(t = τ₁)为0，进一步的氧化过程可能是集流体的氧化。初始条件下为E₁ = E₂，随后变为 E₁ > E₂且U = E₁ − E₂

然而，在正向电流的情况下，电位的快速上升也可能对应于阴极电解质 | 锂界面处锂离子的完全耗尽，如下图所示：

阴极界面处的锂离子浓度降至0，==进一步的还原过程可能是溶剂的还原==。初始条件下为E₁ = E₂，随后变为 E₁ > E₂且U = E₁ − E₂。

对称电池是一种简单且经验性的方法，==用于测试电极在恒电流循环中的稳定性==。对于锂金属电极而言，氧化可能导致表面降解，而金属沉积（锂离子被还原）则可能导致枝晶形成。下图a)展示了稳定电池和稳定电极的电压响应，其电压幅值不随时间变化；图b)则显示了由不稳定电极组成的不稳定电池的电压响应，其电位幅值会随时间发生变化：

对过渡时间的讨论

我们知道了在t = τ时电位会发生阶跃，对应这两种可能的情况：阳极锂金属耗尽或者阴极与电解质的界面处的锂离子耗尽，下面分别讨论一下这两种情况：

假设阳极氧化过程完全遵循方程$\rm Li \longrightarrow Li^+ +e^-$，则过渡时间τ仅取决于电极质量（该质量决定了完全氧化阳极所需的电荷量）以及施加电流。若我们定义单位面积的电荷量与电流密度分别为Q和i，则可得到： $$ \tau=\dfrac{Q}{i} $$ 将单位面积电子数量Q替换为单位面积锂离子数量，可得：

1
2
3
4
5

{% raw %}
$$
\tau=\dfrac{m_{\rm Li}}{{\rm M_{Li}}i}
$$
{% endraw %}

假设阴极还原过程完全遵循$\rm Li^+ +e^- \longrightarrow Li$，则过渡时间τ现可通过 Sand 方程计算得出： $$ \sqrt{\tau}=\dfrac{FC^*_{\rm Li^+}\sqrt{\pi D_{\rm Li^+}}}{2i} $$ 其中$C^*_{\rm Li^+}$和$D_{\rm Li^+}$分别表示电解液中锂离子的本体浓度和扩散系数，了解电流密度如何影响过渡时间，有助于我们更好地理解是哪种机制决定了电位升降所对应的时间过渡特性。

简化版的sand方程： $$ j(\tau)=\dfrac{z_eFc\sqrt{\pi D_{\rm Li^+}}}{2\sqrt{\tau}} $$ 当电流密度超过电解质的$\rm Li^+$传输能力时，电极表面$\rm Li^+$浓度随时间下降，直至耗尽(τ时刻），引发电压骤升（阻塞型极化）。方程表明，τ与j²成反比，与D_Li⁺和c成正比

下图展示了在锂对称电池（Li | Li）上，以三种不同数值、正负相反的恒电流密度阶跃测试时，电压随时间变化的演变过程。

可以看出，当电流密度增大一倍时，两个电流方向的==过渡时间均减半==。这表明上述方程可用于描述过渡时间与时间的依赖关系，且==控制机制为锂的完全氧化过程==。

对称电池可有效用于研究和测试锂电极与固体电解质接触时的稳定性，以及这些电解质防止枝晶形成的能力。通过使用两个相同的电极，研究者无需参比电极即可==轻松研究单一电极 | 电解质界面随时间和循环次数的阳极与阴极行为==。然而，若要==深入研究降解机制，建议采用三电极配置==。

一些基本的术语

锂离子电池的文献中有很多关于电池性能的测试，要怎么判断什么是好的锂电池呢？循环次数

参考资料

美国Gamry公司技术资料，https://www.gamry.com/application-notes/battery-research/testing-lithium-ion-batteries/

Biologic公司，https://www.biologic.net/topics/introducing-symmetric-cells/

研之成理微信公众号，https://mp.weixin.qq.com/s?__biz=MzUxMDMzODg2Ng==&mid=2247522824&idx=1&sn=2bc2f226c7348a1f89b1a11ea91fdf86&chksm=f906a855ce7121434d626bfd0c840c22786744c056b41a557149832e7cd5360af8a15dd12e9c&scene=21#wechat_redirect

Comment and share

IRI等值面图

启动Multiwfn，载入输入文件，我使用的输入文件是.fchk文件，别的如.wfn文件等也可以

使用功能20 Visual study of weak interaction，进行弱相互作用分析，选择4 IRI: Interaction region indicator analysis (Chemistry-Methods, 1, 231)进行IRI分析，选择精度，这里我选的是高精度3 High quality grid, covering whole system, about 1728000 points in total（精度是根据格点数目确定的，对于巨大体系可能高精度也不够用需要手动指定格点数目）

至此，格点计算已经完成，然后我们导出数据3 Output cube files to func1.cub and func2.cub in current folder，借助VMD 分子可视化程序画出更好看的图，将载入文件所在文件夹下导出的func1.cub和func2.cub(这两个文件分别记录的是sign(λ₂)ρ和IRI格点数据)复制到VMD程序的根目录并在根目录下新建一个脚本文件，名为IRIfill.vmd，内容为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

mol new func1.cub
mol addfile func2.cub
mol delrep 0 top
mol representation CPK 1.0 0.3 18.0 16.0
mol addrep top
mol representation Isosurface 1.0 1 0 0 1 1
mol color Volume 0
mol addrep top
mol scaleminmax top 1 -0.04 0.02
mol modstyle 0 top CPK 0.700000 0.300000 18.000000 16.000000
color scale midpoint 0.666
color scale method BGR
color Display Background white
axes location Off
display depthcue off
display rendermode GLSL
light 3 on
color Element N iceblue
mol modcolor 0 top Element

然后启动VMD，在文本窗口中输入source IRIfill.vmd，

Comment and share

we

we $$ WSS=\sum_i w_i(I_i^{exp}-I_i^{calc})^2, w_i=\dfrac{1}{I_i^{exp}} $$

I_i^calc = S_F∑_kL_k|F_k|²S(2θ_i − 2θ_k)P_kA + bkg_i

其中，P_k是优选取向函数，S(2θ_i − 2θ_k)是剖面形状函数，PV=η(HWHM)，HWHM²=Utan²θ + Vtan θ + W $$ P_k=(r^2\cos^2\alpha+\dfrac{\sin^2\alpha}{r})^{-\frac{3}{2}} $$

Comment and share

• Multiwfn版本：3.8(dev)，Last update: 2025-May-5
• VMD版本：1.9.3

VMD版本请尽可能使用1.9.3而不要使用bug较多的1.9.4

本文参考了Sob老师的博文

目标分子的结构优化

本次的目标分子为二茂铁分子

首先建立二茂铁分子的模型，并且手动设置全重叠的D5d点群，这样可以极大的减少优化用时，关键词如下：

1

#p tpssh/def2tzvp opt freq

体系是过渡金属配合物，因此选择tpssh泛函，由于初猜结构比较合适，算得还是挺快的

1
2
3
4
5

                                   -- CHARLES MCCARRY IN "THE GREAT SOUTHWEST"
Job cpu time:       0 days  1 hours 57 minutes 15.2 seconds.
Elapsed time:       0 days  0 hours  5 minutes 58.7 seconds.
File lengths (MBytes):  RWF=    147 Int=      0 D2E=      0 Chk=     16 Scr=      1
Normal termination of Gaussian 16 at Wed May 28 21:43:45 2025.

如下图，优化完成，没有虚频

绘图前的准备工作

如Sob老师所说，绘图之前用户要做的事包括：

1. 把examples.txt和.bat文件拷到Multiwfn可执行文件所在目录
2. 把.bat文件里的VMD路径改成实际路径
3. 把examples.vmd拷到VMD路径下
4. 在vmd.rc里加入上述proc语句

画图

有两种风格的分子表面静电势图(ESP)，分别是表面顶点着色方式和电子密度等值面着色方式的ESP图，我个人更喜欢后者，以后者为例画图吧。

首先运行ESPiso.bat产生电子密度和静电势的cube文件：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451

PS C:\Multiwfn_3.8_dev_bin_Win64> ls

    Directory: C:\Multiwfn_3.8_dev_bin_Win64

Mode                 LastWriteTime         Length Name
----                 -------------         ------ ----
d----           2025/5/15    12:20                examples
-a---           2025/5/23     9:20        8114156 1.fch
-a---           2024/9/21    17:46           3057 使用Multiwfn发表文章必须在正文里进行引用（包括代算）.txt
-a---           2025/5/28    17:51            105 ESPext.bat
-a---           2018/10/7     5:24             32 ESPext.txt
-a---           2025/5/28    17:51            460 ESPiso.bat
-a---          2021/11/18    10:43             26 ESPiso.txt
-a---           2025/5/28    17:51            344 ESPpt.bat
-a---           2018/9/30    13:53             29 ESPpt.txt
-a---            2025/2/9     9:28         200205 How to cite Multiwfn.pdf
-a---            2023/3/2     2:26        2047000 libiomp5md.dll
-a---           2024/8/27    20:00           1254 LICENSE.txt
-a---           2025/4/13     9:17         416416 Multiwfn quick start.pdf
-a---            2025/5/5     2:50       35666432 Multiwfn.exe
-a---           2025/5/16     1:26          20490 settings.ini

PS C:\Multiwfn_3.8_dev_bin_Win64> .\ESPiso.bat

C:\Multiwfn_3.8_dev_bin_Win64>Multiwfn 1.fch -ESPrhoiso 0.001  0<ESPiso.txt
 Multiwfn -- A Multifunctional Wavefunction Analyzer
 Version 3.8(dev), update date: 2025-May-5
 Developer: Tian Lu (Beijing Kein Research Center for Natural Sciences)
 Multiwfn official website: http://sobereva.com/multiwfn
 Multiwfn English forum: http://sobereva.com/wfnbbs
 Multiwfn Chinese forum: http://bbs.keinsci.com/wfn
 ( Number of parallel threads:   4  Current date: 2025-05-30  Time: 20:33:12 )

 Both following papers ***MUST BE CITED IN MAIN TEXT*** if Multiwfn is used:
  Tian Lu, Feiwu Chen, J. Comput. Chem., 33, 580 (2012) DOI: 10.1002/jcc.22885
  Tian Lu, J. Chem. Phys., 161, 082503 (2024) DOI: 10.1063/5.0216272
 See "How to cite Multiwfn.pdf" in Multiwfn binary package for more information

 Please wait...
 Loading various information of the wavefunction
 The highest angular moment basis functions is F
 Loading basis set definition...
 Loading orbitals...
 Converting basis function information to GTF information...
 Back converting basis function information from Cartesian to spherical type...
 Generating density matrix based on SCF orbitals...
 Generating overlap matrix...

 Total/Alpha/Beta electrons:     98.0000     49.0000     49.0000
 Net charge:     0.00000      Expected multiplicity:    1
 Atoms:     23,  Basis functions:    427,  GTFs:    708
 Total energy:   -1651.955643142885 Hartree,   Virial ratio:  2.00246837
 This is a restricted single-determinant wavefunction
 Orbitals from 1 to    49 are occupied
 Title line of this file: Generated by autoGau.

 Loaded 1.fch successfully!

 Formula: H12 C10 Fe1      Total atoms:      23
 Molecule weight:       188.04779 Da
 Point group: C1

 "q": Exit program gracefully          "r": Load a new file
                    ************ Main function menu ************
 0 Show molecular structure and view orbitals
 1 Output all properties at a point       2 Topology analysis
 3 Output and plot specific property in a line
 4 Output and plot specific property in a plane
 5 Output and plot specific property within a spatial region (calc. grid data)
 6 Check & modify wavefunction
 7 Population analysis and calculation of atomic charges
 8 Orbital composition analysis           9 Bond order analysis
 10 Plot total DOS, PDOS, OPDOS, local DOS, COHP and photoelectron spectrum
 11 Plot IR/Raman/UV-Vis/ECD/VCD/ROA/NMR spectrum
 12 Quantitative analysis of molecular surface
 13 Process grid data (No grid data is presented currently)
 14 Adaptive natural density partitioning (AdNDP) analysis
 15 Fuzzy atomic space analysis
 16 Charge decomposition analysis (CDA) and plot orbital interaction diagram
 17 Basin analysis                       18 Electron excitation analysis
 19 Orbital localization analysis        20 Visual study of weak interaction
 21 Energy decomposition analysis        22 Conceptual DFT (CDFT) analysis
 23 ETS-NOCV analysis                    24 (Hyper)polarizability analysis
 25 Electron delocalization and aromaticity analyses
 26 Structure and geometry related analyses
 100 Other functions (Part 1)            200 Other functions (Part 2)
 300 Other functions (Part 3)
 -10 Return to main menu
 -2 Obtain deformation property
 -1 Obtain promolecule property
 0 Set custom operation
             ----------- Available real space functions -----------
 1 Electron density (rho)     2 Gradient norm of rho     3 Laplacian of rho
 4 Value of orbital wavefunction         44 Orbital probability density
 5 Electron spin density
 6 Hamiltonian kinetic energy density K(r)
 7 Lagrangian kinetic energy density G(r)
 8 Electrostatic potential from nuclear charges
 9 Electron localization function (ELF)
 10 Localized orbital locator (LOL)
 11 Local information entropy
 12 Total electrostatic potential (ESP)
 13 Reduced density gradient (RDG)       14 RDG with promolecular approximation
 15 Sign(lambda2)*rho      16 Sign(lambda2)*rho with promolecular approximation
 17 Correlation hole for alpha, ref. point:   0.00000   0.00000   0.00000
 18 Average local ionization energy (ALIE)
 19 Source function, mode: 1, ref. point:   0.00000   0.00000   0.00000
 20 Electron delocal. range func. EDR(r;d)  21 Orbital overlap dist. func. D(r)
 22 Delta-g (promolecular approximation)    23 Delta-g (Hirshfeld partition)
 24 Interaction region indicator (IRI)    25 van der Waals potential (probe=C )
 100 User-defined function (iuserfunc=   -1), see Section 2.7 of manual

 Please select a method to set up grid
 -10 Set extension distance of grid range for mode 1~4, current:  6.000 Bohr
 1 Low quality grid,    covering whole system, about 125000 points in total
 2 Medium quality grid, covering whole system, about 512000 points in total
 3 High quality grid,   covering whole system, about 1728000 points in total
 4 Input the number of points or grid spacing in X,Y,Z, covering whole system
 5 Input original point, grid spacings, and the number of points
 6 Input center coordinate, number of points and extension distance
 7 The same as 6, but input two atoms, the midpoint will be defined as center
 8 Use grid setting of another cube file
 10 Set box of grid data visually using a GUI window
 11 Select a set of atoms, set extension distance around them and grid spacing
 100 Load a set of points from external file
 Coordinate of origin in X,Y,Z is     -12.142516  -10.564851  -10.153107 Bohr
 Coordinate of end point in X,Y,Z is   11.971265   10.285411   10.334542 Bohr
 Grid spacing in X,Y,Z is    0.181307    0.181307    0.181307 Bohr
 Number of points in X,Y,Z is  134  116  114   Total:     1772016
 Note: Virtual orbitals higher than LUMO+9 have been temporarily discarded for saving computational time

 Note: All exponential functions exp(x) with x< -40.000 will be ignored
 Unique GTFs have been constructed. Number of unique GTFs:     708
 Progress: [##################################################]   100.0 %     \
 Calculation of grid data took up wall clock time         4 s
 Note: Previous orbital information has been restored

 Electric dipole moment estimated by integrating electron density
 X component:          0.281807 a.u.      0.716283 Debye
 Y component:         -0.255113 a.u.     -0.648432 Debye
 Z component:          0.140547 a.u.      0.357235 Debye
 Total magnitude:      0.405279 a.u.      1.030117 Debye

 The minimum is  0.60334894E-17 at -12.14252 -10.56485  10.33454 Bohr
 The maximum is  0.95394143E+03 at   0.00503   0.13224   0.00006 Bohr
 Summing up all value and multiply differential element:
   100.725159050150
 Summing up positive value and multiply differential element:
   100.725159050150
 Summing up negative value and multiply differential element:
  0.000000000000000E+000

                   ---------- Post-processing menu ----------
 -1 Show isosurface graph
 0 Return to main menu
 1 Save graph of isosurface to file in current folder
 2 Export data to a Gaussian-type cube file in current folder
 3 Export data to output.txt in current folder
 4 Set the value of isosurface to be shown, current:   0.05000
 5 Multiply all grid data by a factor
 6 Divide all grid data by a factor
 7 Add a value to all grid data
 8 Substract a value from all grid data
 9 Multiply all grid data by Hirshfeld weights of a fragment (can be used to only make isosurface around interested fragment visible)
 Exporting cube file, please wait...
 Done! Grid data has been exported to density.cub in current folder
 Hint: If you want to add input file name as prefix of the outputted cube file, you can set "iaddprefix" in settings.ini to 1

                   ---------- Post-processing menu ----------
 -1 Show isosurface graph
 0 Return to main menu
 1 Save graph of isosurface to file in current folder
 2 Export data to a Gaussian-type cube file in current folder
 3 Export data to output.txt in current folder
 4 Set the value of isosurface to be shown, current:   0.05000
 5 Multiply all grid data by a factor
 6 Divide all grid data by a factor
 7 Add a value to all grid data
 8 Substract a value from all grid data
 9 Multiply all grid data by Hirshfeld weights of a fragment (can be used to only make isosurface around interested fragment visible)

 Note: A set of grid data is presented in memory
 "q": Exit program gracefully          "r": Load a new file
                    ************ Main function menu ************
 0 Show molecular structure and view orbitals
 1 Output all properties at a point       2 Topology analysis
 3 Output and plot specific property in a line
 4 Output and plot specific property in a plane
 5 Output and plot specific property within a spatial region (calc. grid data)
 6 Check & modify wavefunction
 7 Population analysis and calculation of atomic charges
 8 Orbital composition analysis           9 Bond order analysis
 10 Plot total DOS, PDOS, OPDOS, local DOS, COHP and photoelectron spectrum
 11 Plot IR/Raman/UV-Vis/ECD/VCD/ROA/NMR spectrum
 12 Quantitative analysis of molecular surface
 13 Process grid data
 14 Adaptive natural density partitioning (AdNDP) analysis
 15 Fuzzy atomic space analysis
 16 Charge decomposition analysis (CDA) and plot orbital interaction diagram
 17 Basin analysis                       18 Electron excitation analysis
 19 Orbital localization analysis        20 Visual study of weak interaction
 21 Energy decomposition analysis        22 Conceptual DFT (CDFT) analysis
 23 ETS-NOCV analysis                    24 (Hyper)polarizability analysis
 25 Electron delocalization and aromaticity analyses
 26 Structure and geometry related analyses
 100 Other functions (Part 1)            200 Other functions (Part 2)
 300 Other functions (Part 3)
 -10 Return to main menu
 -2 Obtain deformation property
 -1 Obtain promolecule property
 0 Set custom operation
             ----------- Available real space functions -----------
 1 Electron density (rho)     2 Gradient norm of rho     3 Laplacian of rho
 4 Value of orbital wavefunction         44 Orbital probability density
 5 Electron spin density
 6 Hamiltonian kinetic energy density K(r)
 7 Lagrangian kinetic energy density G(r)
 8 Electrostatic potential from nuclear charges
 9 Electron localization function (ELF)
 10 Localized orbital locator (LOL)
 11 Local information entropy
 12 Total electrostatic potential (ESP)
 13 Reduced density gradient (RDG)       14 RDG with promolecular approximation
 15 Sign(lambda2)*rho      16 Sign(lambda2)*rho with promolecular approximation
 17 Correlation hole for alpha, ref. point:   0.00000   0.00000   0.00000
 18 Average local ionization energy (ALIE)
 19 Source function, mode: 1, ref. point:   0.00000   0.00000   0.00000
 20 Electron delocal. range func. EDR(r;d)  21 Orbital overlap dist. func. D(r)
 22 Delta-g (promolecular approximation)    23 Delta-g (Hirshfeld partition)
 24 Interaction region indicator (IRI)    25 van der Waals potential (probe=C )
 100 User-defined function (iuserfunc=   -1), see Section 2.7 of manual

 Please select a method to set up grid
 -10 Set extension distance of grid range for mode 1~4, current:  6.000 Bohr
 1 Low quality grid,    covering whole system, about 125000 points in total
 2 Medium quality grid, covering whole system, about 512000 points in total
 3 High quality grid,   covering whole system, about 1728000 points in total
 4 Input the number of points or grid spacing in X,Y,Z, covering whole system
 5 Input original point, grid spacings, and the number of points
 6 Input center coordinate, number of points and extension distance
 7 The same as 6, but input two atoms, the midpoint will be defined as center
 8 Use grid setting of another cube file
 10 Set box of grid data visually using a GUI window
 11 Select a set of atoms, set extension distance around them and grid spacing
 100 Load a set of points from external file
 Coordinate of origin in X,Y,Z is     -12.142516  -10.564851  -10.153107 Bohr
 Coordinate of end point in X,Y,Z is   11.789958   10.321672   10.298280 Bohr
 Grid spacing in X,Y,Z is    0.435136    0.435136    0.435136 Bohr
 Number of points in X,Y,Z is   56   49   48   Total:      131712
 Note: Virtual orbitals higher than LUMO+9 have been temporarily discarded for saving computational time

 Unique GTFs have been constructed. Number of unique GTFs:     708
 Note: ESP will be calculated only for the grids around isosurface of electron density of   0.001000 a.u.
 Detecting the grids for calculating ESP...
 Number of grids to calculate ESP:       10900

 Initializing LIBRETA library (fast version) for ESP evaluation ...
 LIBRETA library has been successfully initialized!

 NOTE: The ESP evaluation code based on LIBRETA library is being used. Please cite Multiwfn original papers (J. Comput. Chem., 33, 580-592 (2012) and J. Chem. Phys., 161, 082503 (2024)) and the paper describing the efficient ESP evaluation algorithm adopted by Multiwfn (Phys. Chem. Chem. Phys., 23, 20323 (2021))
 Progress: [##################################################]   100.0 %     \
 Setting ESP of the grids neighbouring to boundary grids...
 Calculation of grid data took up wall clock time        13 s
 Note: Previous orbital information has been restored

 The minimum is -0.32612104E-01 at  -0.39385  -0.99186   4.64151 Bohr
 The maximum is  0.52214027E-01 at   2.21697  -2.29727  -5.36661 Bohr
 Summing up all value and multiply differential element:
   3.07509251072690
 Summing up positive value and multiply differential element:
   8.67387639395045
 Summing up negative value and multiply differential element:
  -5.59878388322356

                   ---------- Post-processing menu ----------
 -1 Show isosurface graph
 0 Return to main menu
 1 Save graph of isosurface to file in current folder
 2 Export data to a Gaussian-type cube file in current folder
 3 Export data to output.txt in current folder
 4 Set the value of isosurface to be shown, current:   0.05000
 5 Multiply all grid data by a factor
 6 Divide all grid data by a factor
 7 Add a value to all grid data
 8 Substract a value from all grid data
 9 Multiply all grid data by Hirshfeld weights of a fragment (can be used to only make isosurface around interested fragment visible)
 Exporting cube file, please wait...
 Done! Grid data has been exported to totesp.cub in current folder
 Hint: If you want to add input file name as prefix of the outputted cube file, you can set "iaddprefix" in settings.ini to 1

                   ---------- Post-processing menu ----------
 -1 Show isosurface graph
 0 Return to main menu
 1 Save graph of isosurface to file in current folder
 2 Export data to a Gaussian-type cube file in current folder
 3 Export data to output.txt in current folder
 4 Set the value of isosurface to be shown, current:   0.05000
 5 Multiply all grid data by a factor
 6 Divide all grid data by a factor
 7 Add a value to all grid data
 8 Substract a value from all grid data
 9 Multiply all grid data by Hirshfeld weights of a fragment (can be used to only make isosurface around interested fragment visible)
forrtl: severe (24): end-of-file during read, unit -4, file CONIN$
Image              PC                Routine            Line        Source
Multiwfn.exe       00007FF60135CEB4  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF600F5B65C  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF60127BDD6  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601CCFB7B  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601FE53F0  Unknown               Unknown  Unknown
KERNEL32.DLL       00007FF9259B7374  Unknown               Unknown  Unknown
ntdll.dll          00007FF92667CC91  Unknown               Unknown  Unknown

C:\Multiwfn_3.8_dev_bin_Win64>move /Y density.cub density1.cub
移动了         1 个文件。

C:\Multiwfn_3.8_dev_bin_Win64>move /Y totesp.cub ESP1.cub
移动了         1 个文件。

C:\Multiwfn_3.8_dev_bin_Win64>Multiwfn 2.fch -ESPrhoiso 0.001  0<ESPiso.txt
 Error: Unable to find the input file you specified in argument

 Multiwfn -- A Multifunctional Wavefunction Analyzer
 Version 3.8(dev), update date: 2025-May-5
 Developer: Tian Lu (Beijing Kein Research Center for Natural Sciences)
 Multiwfn official website: http://sobereva.com/multiwfn
 Multiwfn English forum: http://sobereva.com/wfnbbs
 Multiwfn Chinese forum: http://bbs.keinsci.com/wfn
 ( Number of parallel threads:   4  Current date: 2025-05-30  Time: 20:33:30 )

 Both following papers ***MUST BE CITED IN MAIN TEXT*** if Multiwfn is used:
  Tian Lu, Feiwu Chen, J. Comput. Chem., 33, 580 (2012) DOI: 10.1002/jcc.22885
  Tian Lu, J. Chem. Phys., 161, 082503 (2024) DOI: 10.1063/5.0216272
 See "How to cite Multiwfn.pdf" in Multiwfn binary package for more information

 Now input file path, for example, E:\Planetarian\Yumemi_Hoshino.mwfn
 (.wfn/wfn/wfx/fch/molden/pdb/xyz/mol2/cif/cub... see Section 2.5 of manual)
 Hint: Pressing ENTER button directly can select a file in a GUI window. To reload the past file, inputting "o". Input such as ?miku.fch can open the miku.fch in the same folder as the past file
"5" cannot be found, input again
"1" cannot be found, input again
"3" cannot be found, input again
"2" cannot be found, input again
"0" cannot be found, input again
"5" cannot be found, input again
"12" cannot be found, input again
"1" cannot be found, input again
"2" cannot be found, input again
forrtl: severe (24): end-of-file during read, unit -4, file CONIN$
Image              PC                Routine            Line        Source
Multiwfn.exe       00007FF601370B01  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF60127CBB3  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601CCFB7B  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601FE53F0  Unknown               Unknown  Unknown
KERNEL32.DLL       00007FF9259B7374  Unknown               Unknown  Unknown
ntdll.dll          00007FF92667CC91  Unknown               Unknown  Unknown

C:\Multiwfn_3.8_dev_bin_Win64>move /Y density.cub density2.cub
系统找不到指定的文件。

C:\Multiwfn_3.8_dev_bin_Win64>move /Y totesp.cub ESP2.cub
系统找不到指定的文件。

C:\Multiwfn_3.8_dev_bin_Win64>Multiwfn 3.fch -ESPrhoiso 0.001  0<ESPiso.txt
 Error: Unable to find the input file you specified in argument

 Multiwfn -- A Multifunctional Wavefunction Analyzer
 Version 3.8(dev), update date: 2025-May-5
 Developer: Tian Lu (Beijing Kein Research Center for Natural Sciences)
 Multiwfn official website: http://sobereva.com/multiwfn
 Multiwfn English forum: http://sobereva.com/wfnbbs
 Multiwfn Chinese forum: http://bbs.keinsci.com/wfn
 ( Number of parallel threads:   4  Current date: 2025-05-30  Time: 20:33:30 )

 Both following papers ***MUST BE CITED IN MAIN TEXT*** if Multiwfn is used:
  Tian Lu, Feiwu Chen, J. Comput. Chem., 33, 580 (2012) DOI: 10.1002/jcc.22885
  Tian Lu, J. Chem. Phys., 161, 082503 (2024) DOI: 10.1063/5.0216272
 See "How to cite Multiwfn.pdf" in Multiwfn binary package for more information

 Now input file path, for example, E:\Planetarian\Yumemi_Hoshino.mwfn
 (.wfn/wfn/wfx/fch/molden/pdb/xyz/mol2/cif/cub... see Section 2.5 of manual)
 Hint: Pressing ENTER button directly can select a file in a GUI window. To reload the past file, inputting "o". Input such as ?miku.fch can open the miku.fch in the same folder as the past file
"5" cannot be found, input again
"1" cannot be found, input again
"3" cannot be found, input again
"2" cannot be found, input again
"0" cannot be found, input again
"5" cannot be found, input again
"12" cannot be found, input again
"1" cannot be found, input again
"2" cannot be found, input again
forrtl: severe (24): end-of-file during read, unit -4, file CONIN$
Image              PC                Routine            Line        Source
Multiwfn.exe       00007FF601370B01  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF60127CBB3  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601CCFB7B  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601FE53F0  Unknown               Unknown  Unknown
KERNEL32.DLL       00007FF9259B7374  Unknown               Unknown  Unknown
ntdll.dll          00007FF92667CC91  Unknown               Unknown  Unknown

C:\Multiwfn_3.8_dev_bin_Win64>move /Y density.cub density3.cub
系统找不到指定的文件。

C:\Multiwfn_3.8_dev_bin_Win64>move /Y totesp.cub ESP3.cub
系统找不到指定的文件。

C:\Multiwfn_3.8_dev_bin_Win64>Multiwfn 4.fch -ESPrhoiso 0.001  0<ESPiso.txt
 Error: Unable to find the input file you specified in argument

 Multiwfn -- A Multifunctional Wavefunction Analyzer
 Version 3.8(dev), update date: 2025-May-5
 Developer: Tian Lu (Beijing Kein Research Center for Natural Sciences)
 Multiwfn official website: http://sobereva.com/multiwfn
 Multiwfn English forum: http://sobereva.com/wfnbbs
 Multiwfn Chinese forum: http://bbs.keinsci.com/wfn
 ( Number of parallel threads:   4  Current date: 2025-05-30  Time: 20:33:30 )

 Both following papers ***MUST BE CITED IN MAIN TEXT*** if Multiwfn is used:
  Tian Lu, Feiwu Chen, J. Comput. Chem., 33, 580 (2012) DOI: 10.1002/jcc.22885
  Tian Lu, J. Chem. Phys., 161, 082503 (2024) DOI: 10.1063/5.0216272
 See "How to cite Multiwfn.pdf" in Multiwfn binary package for more information

 Now input file path, for example, E:\Planetarian\Yumemi_Hoshino.mwfn
 (.wfn/wfn/wfx/fch/molden/pdb/xyz/mol2/cif/cub... see Section 2.5 of manual)
 Hint: Pressing ENTER button directly can select a file in a GUI window. To reload the past file, inputting "o". Input such as ?miku.fch can open the miku.fch in the same folder as the past file
"5" cannot be found, input again
"1" cannot be found, input again
"3" cannot be found, input again
"2" cannot be found, input again
"0" cannot be found, input again
"5" cannot be found, input again
"12" cannot be found, input again
"1" cannot be found, input again
"2" cannot be found, input again
forrtl: severe (24): end-of-file during read, unit -4, file CONIN$
Image              PC                Routine            Line        Source
Multiwfn.exe       00007FF601370B01  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF60127CBB3  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601CCFB7B  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601FE53F0  Unknown               Unknown  Unknown
KERNEL32.DLL       00007FF9259B7374  Unknown               Unknown  Unknown
ntdll.dll          00007FF92667CC91  Unknown               Unknown  Unknown

C:\Multiwfn_3.8_dev_bin_Win64>move /Y density.cub density4.cub
系统找不到指定的文件。

C:\Multiwfn_3.8_dev_bin_Win64>move /Y totesp.cub ESP4.cub
系统找不到指定的文件。

C:\Multiwfn_3.8_dev_bin_Win64>move /Y *.cub "C:\VMD"
C:\Multiwfn_3.8_dev_bin_Win64\density1.cub
C:\Multiwfn_3.8_dev_bin_Win64\ESP1.cub
移动了         2 个文件。

然后打开VMD软件，输入iso即可画图

再输入ext即可显示极值点

上图中的材质已被我改为transparent，因此能看到背面的极值点，修改方法可以通过Grapgics-Materials进行修改，然而答主修改时没有任何变化，于是在ESPiso.vmd中直接将EdgyGlass改成了Transparent

怎么添加刻度呢？操作如下图

接下来是自定义图片的导出，File-Render，渲染器选择Tachyon，开始渲染，然后在VMD的根目录下得到了文件vmdscene.dat，接下来调用渲染器进行渲染，渲染的.bat脚本是：

1

tachyon_WIN32.exe vmdscene.dat -aasamples 24 -mediumshade -trans_vmd -res 1024 768 -format BMP -o vmdscene.bmp

其中的1024是长，768是宽，可以改成想要的分辨率

Comment and share

什么是电化学阻抗谱？

有一种叫做黑箱的未知特性系统。可采用一系列测量方法，通过施加输入信号并检测输出响应来探查黑箱。例如，假设将黑箱置于暗室中，然后施加特定波长的光照。若观察到响应信号（如电流），则可判定箱内物质具有光活性。输入与输出之间的关系称为“传递函数”，阻抗谱是传递函数的一种特例。

阻抗：电路中电阻、电感、电容对交流电个阻碍作用个统称。阻抗是一个复数，实部称为电阻，虚部称为电抗。电化学过程主要涉及电荷在电极和电解质界面上的传输和存储。例如，电极表面的双电层效应会形成双电层电容，电解质中的离子扩散也会导致电容行为。这些机制使电容成为电化学系统中主要的电抗元件。 相比之下，电感通常与磁场的变化相关，比如电流通过线圈时产生的磁场。然而，在常规的电化学系统中，磁场变化通常不显著，因此电感的影响非常小，通常电化学系统中考虑电阻、容抗而不考虑感抗。

电化学阻抗谱（EIS）通过施加交流电压（或电流）并测量产生的电流（或电压），能够研究电极和电化学池的电化学特性，利用等效电路进行分析，可以确定与电极电化学反应相关的重要参数，如电荷转移电阻、双电层电容和Warburg阻抗

阻抗的复数表达形式

交流电（AC）电路中的电流、电压以及阻抗通常使用复数来表示，主要原因有两个：

1. 交流信号（以及许多其他正弦波现象）具有一个幅值和一个相位，它们分别与复数的模和幅角非常相似
2. 复数的基本运算，如加法、减法、乘法和除法，更容易执行和编程

注意：

• 因为符号i 在交流电路中用于表示电流，这里我们用j作为虚数单位，即j² = −1
• 符号 ℜe代表复数的实部

预备知识

复数的标准形式 Z = a + bj 复数的模 $$ \lvert Z \rvert =\sqrt {a^2+b^2} $$

欧拉公式 e^ix = cos x + isin x 复数的指数形式 Z = re^jθ = rcos θ + rjsin θ 取实数部分： ℜe(re^jθ) = rcos θ

单变量复函数f(t)的导数的实部等于f(t)的实部的导数，证明如下： $$ $$ **复函数$f(t)$的不定积分的实部等于$f(t)$的实部的不定积分，证明如下** $$ $$

复数形式的阻抗

• 电压源

v(t) = V₀cos (ωt) = ℜe(V₀cos (ωt) + V₀jsin (ωt)) = ℜe(V₀e^jωt)

• 电阻

对于如下图的交变电路

定义Z_R = R，那么电流$i=\Re e(\dfrac{V_R}{R})$或$I=\dfrac{V_R}{Z_R}$

• 电容

电容上电流的定义为：$I=C(\dfrac{d v(t)}{dt})$，则有如下推导： $$ \begin{align} v(t)=&\dfrac{1}{C}\int i dt\\ v(t)=&V_0\cos (\omega t)=\Re e(V_0e^{j\omega t})\\ \Re e(V_0e^{j\omega t})=&\dfrac{1}{C}\int idt\\ \dfrac{d}{dt}(\Re e(V_0e^{j\omega t}))=&\dfrac{d}{dt}(\dfrac{1}{C}\int idt)\\ \Re e(j\omega V_0e^{j\omega t})=&\dfrac{1}{C}i\\ i=&\Re e(j\omega CV_0e^{j\omega t}) \end{align} $$ 电容器C两端的电压为V_C = V₀e^jωt，那么电容对电流的阻碍作用$Z_C=\dfrac{1}{j\omega C}=-\dfrac{j}{\omega C}$

与电阻的阻抗（电阻）Z_R不同，电容的阻抗（容抗）Z_C是一个纯虚数

• 电感

在电化学系统中，电感并不常见

电感L上的电流i和电压v(t)之间的关系为：$v(t)=L\dfrac{di}{dt}$，则有： $$ \begin{align} \Re e(V_0e^{j\omega t})=&L\dfrac{di}{dt}\\ \int \Re e(V_0e^{j\omega t})=&\int L\dfrac{di}{dt}dt\\ \Re e(\dfrac{1}{j\omega}V_0e^{j\omega t})=&Li\\ i=&\Re e(\dfrac{1}{j\omega L}V_0e^{j\omega t}) \end{align} $$ 电感的对电流的阻碍作用为：Z_L = jωL，类似的，电感的阻抗（感抗）也是一个纯虚数

交变电路中的电阻和电容

电压源为：E(t) = E₀sin (ωt)

电容表现出与电阻相似的行为，即电压与电流的幅值比恒定，但电流相位相对于外加电压存在-90°偏移:$I=E_0\omega C\sin(\omega t-\dfrac{\pi}{2})$

电阻和电容的电流响应可用两个参数表征：

• 阻抗模量|Z|
  • 电阻：|Z_R|恒定
  • 电容：$\lvert Z_C \rvert=\dfrac{1}\omega C{}$
• 电压与电流之间的相位差θ
  • 电阻：θ = 0^∘
  • 电容：θ = −90^∘

使用阻抗表达式Z = |Z|(cos θ + sin θ)，在交流测量中可以统一处理这两种电路元件

电阻和电容的串联电路比较简单，而并联的情况稍复杂，讨论一下：

总阻抗： $$ \begin{align} \dfrac{1}{Z}=&\dfrac{1}{Z'}+\dfrac{1}{Z''}\\ =&\dfrac{1}{R}+j\omega C\\ Z=&\dfrac{1}{\dfrac{1}{R}+j\omega C}\\ =&\dfrac{R}{1+(\omega CR)^2}-\dfrac{\omega CR^2 j}{1+(\omega CR)^2}\\ \lvert Z \rvert=&\sqrt{Z'^2+Z''^2}\\ =&\dfrac{R}{\sqrt{1+(\omega CR)^2}} \end{align} $$ 总结一下串联和并联 RC 电路的复阻抗特性：

阻抗的图解

如上所述，阻抗Z既可用阻抗模值|Z(ω)|和相位角θ(ω)表示，也可用阻抗的实部Z^′(ω)和虚部Z^″(ω)表示，分别对应两种图：Bode图和Nyquist图

串联RC电路的阻抗是Z_R与Z_C之和。因此在高频区，由于Z_R ≫ Z_C ，电路阻抗等于电阻R；而在低频区，由于Z_R ≪ Z_C ，阻抗等于Z_C 。如Bode图所示，串联RC电路在高频区表现为电阻特性（|Z|与频率无关，相位角θ为 0°），在低频区表现为电容特性（|Z|与频率成反比，相位角θ为−90°）。Nyquist图中，串联 RC 电路的阻抗轨迹为一条从实轴垂直延伸的直线，与实轴交点的值等于 R。该电路的时间常数为RC，意味着当频率$ω=\dfrac{1}{RC}$时，|Z_R| = |Z_C|且相位角θ = −45°

并联RC电路的阻抗Z通过公式$\frac{1}{Z}=\frac{1}{Z_R}+\frac{1}{Z_C}$计算。因此，并联 RC 电路，阻抗较小的电路元件特性会显现。在高频区，并联RC电路的Z 等于Z_C ；而在低频区则等于Z_R，这与串联RC电路的情况相反。在并联 RC 电路的Nyquist图中，阻抗轨迹 Z 呈现为半圆形，圆弧直径等于 R。该电路的时间常数τ = RC，与串联电路相同。当频率$\omega =\frac{1}{RC}$时，|Z_R| = |Z_C|且相位角θ = −45°，此时对应Nyquist图中圆弧的顶点。

与Bode图相比，Nyquist图更适合通过轨迹形状来辨识电路结构。但需指出的是，当两个并联 RC 电路的时间常数相近时，即便其电阻和电容值不同，它们的半圆也难以区分。如下图所示，当两个并联 RC 电路的时间常数仅相差10倍时，可观察到由两段圆弧重叠形成的畸变半圆。只有当时间常数差异达到 100 倍以上时，两段圆弧才能被清晰区分。因此在构建等效电路时，仅依Nyquist图并不可取，必须验证该等效电路能否合理解释目标电化学现象中的每个基本过程

下面列举一些常见的电路模型以及相应的Nyquist、Bode图吧

• 电感的Nyquist图在第四象限
• “电阻等价于点”
• “电容等价于垂直线”

• 电容和电阻的并联反映到Nyquist图上即为一个半圆弧
• 有几个电阻-电容并联电路就有可能出现几个半圆弧
• 电阻$\rm{R}_0$+电阻-电容并联电路就等价于将半圆弧向右平移$\rm{R_0}$个单位

定义法拉第阻抗Z_f = R_ct + Z_W

• 上图左侧为理想高频部分，此时电荷转移电阻R_ct很大，阻抗由电荷转移主导
• 上图右侧为理想低频部分，此时“电容断路（电容特性：通高频阻低频）”，阻抗由物质扩散主导（溶液离子扩散、固体电解质中离子扩散等等）
• 如果是一个频率非常广泛的、既包含高频又包含低频的电化学阻抗测试，那么Nyquist图中可能同时出现高频部分的半圆弧和低频部分的直线，如下图所示

电极的等效电路

在电化学阻抗谱（EIS）中，电极的电化学行为通常由包含电阻和电容的等效电路表示。其中，Randles 电路常被用作电极的等效电路模型

• 定义法拉第阻抗Z_f = R_ct + Z_W
• 溶液电阻R_s：即使使用参比电极，IR降也不可忽略。电解质传导过程中涉及的电阻称为溶液电阻（R_s ）， 可通过Nyquist奎斯特图中高频极限处与实轴的交点获得
• 电荷转移电阻R_ct：当过电位较大时，电极发生极化，过电位与电流之间呈现指数依赖，而当电极电位接近平衡电位时，过电位与电流的指数关系可简化为线性关系。因此，在极低过电位区域（<10 mV）内，电极表现为电阻特性，其电阻值称为电荷转移电阻，对应于Nyquist奎斯特图中半圆的直径。由于电荷转移电阻的概念仅适用于获得线性关系的极小过电位区域，因此电化学阻抗谱中的电压振幅通常设置为≤10 mV。
• 双电层电容C_dl：电极与电解质界面处的电双层类似于电容器，虽然无法直接从Nyquist图中确定C_dl 的数值，但如果确定了R_ct ，则可通过以下关系式利用圆弧顶点频率对应并联电路的时间常数计算出C_dl ：τ = R_ct × C_dl
• Wurburg阻抗：源于电解液中活性物质的扩散过程。在半无限扩散条件下，Z_W可通过以下方程表示：

$$ Z_W=\dfrac{RT}{n^2F^2Ac\sqrt{D}}\cdot\dfrac{1}{\sqrt{\omega}}\cdot\dfrac{1-j}{\sqrt{2}} $$

其中，n为反应中转移的电子数，A 为电极的几何面积，c为活性物质的浓度，D为扩散系数。Wurburg阻抗可用一个电容器与电阻器无限级联的电路表示。在Wurburg阻抗中，Z’和-Z’‘相等，因此在Nyquist图中可观察到一条斜率为-45°的直线。通过绘制 Z’和 Z’’随ω变化的关系曲线，可从两条平行直线的斜率确定扩散系数。值得注意的是，该式源自能斯特方程中浓度电势依赖项($\dfrac{\partial E}{\partial c}$)。对于不符合能斯特方程的电极（如发生固态扩散的电极），必须考虑($\dfrac{\partial E}{\partial c}$)项才能计算扩散系数

如下图所示，当扩散系数 D 增大时，活性物质能快速响应电极表面浓度变化而迅速补充至电极界面。因此较大的 D 值会导致Z_W值减小。反之，Z_W随D值减小而增大，从而使直线向低频方向移动。最终在Bode图中，Z_W的线性区会与并联RC电路引起的|Z|变化区域重叠。这一现象在下下图（😛）所示的Nyquist图中更为明显：当 D 值较大时，-45°斜率的直线与圆弧可清晰区分；但随着 D 值减小，圆弧与直线开始重叠。这种重叠现象表明活性物质的扩散速度低于电荷转移反应的时间常数，本质上意味着扩散过程成为速率控制步骤。

一个性能良好的固态电解质材料，其离子电导率应当较大，扩散系数D应当较大，那么Nyquist图应当接近于a或b的形貌

从理想电极到真实电极——常相位元件CPE

CPE被定义为： $$ Z_{CPE}=\dfrac{1}{Y_O(j\omega)^n} $$

θ = 90^∘(1 − n)

其中θ是与理想情况（φ=90°）的相位偏差，n 是从 0 到 1 的常数，定义了与理想行为的偏差

两个极端情况：

• n=1，CPE表现为理想电容器
• n=0，CPE 表现为电阻器
• n=0.5，$Z(\omega)=\dfrac{1}{Y_O}\sqrt{j\omega}$，即等价于 Warburg 阻抗

用传输线模型来直观理解扩散阻抗

Warburg阻抗是一个复杂的元素，它表示氧化还原物质向电极表面的质量传递，并被描绘为奈奎斯特图的低频范围上的-45°线。这种行为是指化学物质的时间相关（非稳态）半无限扩散，如下图A所示，在电极/电解质界面处存在单个边界，距离 x=0。朝向本体溶液，并且在静止条件下，扩散层扩展到由电化学电池的尺寸限定的无限长度（x→∞），同时浓度梯度随时间下降。

将扩散层视为由无数个 “电阻（R_d）- 电容（C_d）” 单元串联，其中：

• R_d：单位长度的扩散电阻（反映扩散阻力）

• C_d：单位长度的 “化学电容”（反映浓度存储能力）

• 半无限扩散——无限长传输线，阻抗为Warburg直线

• 有限扩散——有限长传输线，末端接电阻（投射）或电容（反射）

  • 当有限扩散区对扩散物质可渗透时，随着时间和电流流过电化学电池，建立了稳态浓度梯度（dC/dx=常数）。这是透射或可渗透边界的情况。
  • 当有限扩散区在一段时间后对扩散物质不渗透时，电荷转移是不可能的，并且扩散物质的浓度梯度变为零（dC/dx=0）。这是反射或不渗透边界的情况

==具有透射边界的有限长度线性扩散的典型例子包括固体氧化物或聚合物膜燃料电池，和旋转盘电极实验中的扩散层==

固态电解质的阻抗谱分析

以陶瓷固态电解质为例，假设其仅存在单离子传输而不具有电子导电性

在极高频率区域（>1 MHz）的阻抗数据对于通过 EIS 测量分析固体电解质至关重要，这与其它材料的 EIS 分析存在显著差异。若某固体电解质的体相离子电导率为10⁻⁴S ⋅ cm⁻¹ ，且离子迁移相关电容为10⁻¹¹F ⋅ cm⁻² ，则对应1/(2π × τ)（τ：时间常数）的特征频率约为∼10⁷Hz（10 MHz）。若未对该频率区域采取特殊防护措施进行EIS测量，高频区的阻抗数据将会缺失，仅能观测到部分圆弧。在此情况下，EIS数据中包含的来自电缆和电池的电子电阻分量无法消除。此外，由晶粒内部电阻和晶界电阻构成的离子电阻也无法分离，从而导致固体电解质重要信息的丢失。

由于大多数陶瓷固体电解质是由不同尺寸和取向的晶粒聚集而成的多晶材料，除了晶粒内部电阻外，离子在晶粒间界面处的传输电阻对离子电导率具有重要影响。晶粒内部电阻$R_{\rm{gi}}$ （又称体电阻）是材料固有的离子传输阻力。因此，若测试样品不含孔隙，通过晶粒内部电阻倒数计算得到的体离子电导率等于该固体电解质单晶的离子电导率。晶界电阻$R_{\rm{gb}}$ 表示离子穿越晶粒间界面（晶界）及/或位于晶界的化学组成不同微相时所受的传输阻力。下图展示了陶瓷固体电解质的等效电路及Nyquist图，该电路由两个串联的并联 RC 回路构成：一个电阻对应$R_{\rm{gi}}$ ，另一个对应$R_{\rm{gb}}$ 。由于两个并联RC电路具有不同的时间常数（τ = R × C），电化学阻抗谱能够分别测定这两个电阻值。

由于开发和改进陶瓷固体电解质最有效的方法是通过降低$\rm{R_{gb}}$ 使离子电导率接近单晶水平，因此分离$\rm{R_{gi}}$ 和$\rm{R_{gb}}$ 并分别评估它们至关重要。要通过高频区半圆直径确定$\rm{R_{gi}}$ ，应通过增加电解质厚度、减小电极面积或降低测量温度来降低特征频率。然而这些方法存在局限性，例如当电极面积小于样品几何面积时会干扰样品中的电流分布，因此需要测量直至高频区的阻抗数据。本文介绍了我们测量固体电解质高达 100 MHz 高频阻抗的方法。

在固体电解质的电化学阻抗谱（EIS）测量中，通常采用电极/电解质/电极的三明治结构组件。这种配置的优势在于能够实现简单的两端测量（注：设备与电极之间需连接两根电压线和两根电流线，共计四根电缆）。采用电极/电解质/电极结构进行EIS测量时，电极材料的选择至关重要。电极材料需满足以下必要条件：能够构建稳定界面且不与电解质发生化学反应，同时具备足够的电子导电性。金或铂等贵金属常被用作阻塞电极，某些情况下也会采用能与移动离子发生可逆反应的可逆电极。圆盘状固体电解质样品通常通过粉末压片成型，必要时还需进行煅烧处理。 样品的相对密度最好大于 90%，理想情况下应超过 95%。

在电极涂覆工艺中，有两个关键因素会影响固体电解质的阻抗测量结果。其一是电极厚度：过薄的电极可能导致阻抗数据偏离原始值，进而造成对固体电解质离子传输行为的误判。其二是固体电解质的表面状态，该因素会显著影响阻抗测量。本节将详细阐述这两个因素如何影响阻抗数据。

下图展示了不同电极厚度及固体电解质表面是否抛光处理的氟离子导体Nyquist图。所有情况下均采用溅射法沉积直径为7mm 的铂电极。阻抗测试使用 Keysight E4990A 阻抗分析仪在室温空气环境下完成。所有Nyquist图均呈现两个特征：由晶粒内部与晶界电阻（具有相近时间常数）形成的畸变弧线（两个半圆弧融合），以及低频区对应阻塞电极与固体电解质界面电容特性的直线段。不同电极厚度的Nyquist图存在显著差异：对于 200 nm 厚电极（图 2a），离子传输电阻（晶粒内部与晶界电阻）形成的阻抗弧起始于坐标原点；而 50 nm 薄电极（图 2b）的阻抗弧则发生水平偏移，且由阻抗弧外推得到的高频极限阻抗值偏离原点。 这种因电极厚度不同而产生的差异可归因于样品中的电流分布；对于较薄的电极，其横向电阻不容忽视

如上图，当样品经 2000 目砂纸抛光并沉积 50 纳米铂电极后，所得阻抗弧严重偏离原点。这表明抛光过程引入了除离子传输过程之外的额外电阻。此外，抛光样品的离子传导过程对应阻抗弧略有扩大，直线斜率趋于平缓。这种斜率变化反映了抛光前后样品阻塞行为的差异：抛光样品在<1 MHz 频率区间出现与阻塞行为重叠的寄生电阻。固态电解质抛光前后的阻抗行为变化归因于抛光过程中（为去除煅烧样品表面偏析层而进行的）样品表面劣化。

为消除电缆电感对高频测量的干扰，设备与电解池（样品）之间的连接电缆应尽可能缩短。建议采用特性阻抗为 50Ω的同轴电缆，并将其紧贴样品安装。此外，必须合理设置电流回流路径。对于高达 100MHz 的高频测量，需确保电缆不会引起相位差畸变。在进行高频测量前，应对系统进行开路/短路/负载校准

离子电导率的计算

$$ \sigma'(\omega)=\dfrac{L}{S}\times\dfrac{Z'(\omega)}{\vert Z'(\omega) \rvert ^2} $$

对于使用固态电解质模具进行电化学阻抗谱测试的情况，离子电导率的计算公式为： $$ \sigma =\dfrac{L}{S\times R} $$ 其中L为固态电解质pellet圆片厚度，S为圆片pellet的面积，R为电阻（Nyquist图拐点处的横坐标）

实战演练——以$\rm{Li_{5.5}PS_{4.5}Cl_{1.5}}$的电化学阻抗谱图为例

这是之前测试的图，简单画了一下，2Tons（249MPa）压力下的测试结果。可以看到，这与常见的阻抗谱差别非常大，并没有出现超高频的半圆弧，只出现了低频的扩散部分和中高频的拐点。这是收阻抗仪的最高频率所限，使用的阻抗仪最高频只能达到约1MHz左右，对于10mS/cm的固态电解质，想要完整的测出半圆弧可能需要100MHz甚至更高的频率

• 由于设备频率上限较低，只能得到粒子抵抗和粒界阻抗的总值
• 活化能分析中无法判断随温度变化特性值发生变化的是粒子阻抗还是粒界阻抗
• 在通过变温方式尝试分离粒子阻抗和粒界阻抗时，可控温度范围有限，且无法对变温引起的参数条件变化进行校正

==此外，交变电流源的最高频率不等于最高有效频率==，有效的最高频率常常受限于测试端子类型和连线长度，在进行变温控制的时候，==夹具及导线材料的介电常数因为温度变化而变化，如果不及时矫正阻抗会使结果产生较大偏离==

下图展示了超高有效频率下的LLTO型固态电解质的阻抗谱

参考资料

1. Lazanas, A. Ch.; Prodromidis, M. I. Electrochemical Impedance Spectroscopy─a Tutorial. ACS Meas. Sci. Au 2023, 3 (3), 162–193. https://doi.org/10.1021/acsmeasuresciau.2c00070.
2. Lazanas, A. Ch.; Prodromidis, M. I. Correction to “Electrochemical Impedance Spectroscopy─a Tutorial.” ACS Meas. Sci. Au 2025, 5 (1), 156–156. https://doi.org/10.1021/acsmeasuresciau.5c00007.
3. Ariyoshi, K.; Mineshige, A.; Takeno, M.; Fukutsuka, T.; Abe, T.; Uchida, S.; Siroma, Z. Electrochemical Impedance Spectroscopy Part 2: Applications. Electrochemistry 2022, 90 (10), 102008–102008. https://doi.org/10.5796/electrochemistry.22-66080.
4. Ariyoshi, K.; Siroma, Z.; Mineshige, A.; Takeno, M.; Fukutsuka, T.; Abe, T.; Uchida, S. Electrochemical Impedance Spectroscopy Part 1: Fundamentals. Electrochemistry 2022, 90 (10), 102007–102007. https://doi.org/10.5796/electrochemistry.22-66071.
5. Hu, W.; Peng, Y.; Wei, Y.; Yang, Y. Application of Electrochemical Impedance Spectroscopy to Degradation and Aging Research of Lithium-Ion Batteries. J. Phys. Chem. C 2023, 127 (9), 4465–4495. https://doi.org/10.1021/acs.jpcc.3c00033.

Comment and share

实验部分

球磨

• 球磨的原理是什么？
• 球磨的作用是什么？
• 200rpm合不出来，300rpm就能合出来？

X-ray diffraction

• X射线是什么？
• X射线衍射的原理和用途是什么？
• XRD谱图的横坐标和纵坐标分别是什么意义？
• XRD为什么会产生峰？峰是什么意义？
• 横坐标2theta角度为什么通常选择20^∘~80^∘
• 单晶衍射和粉末衍射分别是什么？
• 关于消光？

近年来锂离子电池大规模应用的快速发展要求系统地了解锂离子电池的退化机制，这在很大程度上取决于先进表征方法的利用和发展。基于 X 射线的技术，如 X 射线吸收光谱（XAS）、X 射线衍射（XRD）和 X 射线光电子能谱（XPS），可以确定化学成分和元素的价态和晶体结构的变化。 电子和扫描探针显微技术，（3）如扫描电子显微镜（SEM）、透射电子显微镜（TEM）和低温透射电镜，可以观察电极的形貌变化和固体电解质界面（SEI）的厚度。基于滴定的技术，如滴定气相色谱（TGC）（4）和质谱滴定（MST），（5,6）可以定量死 Li 金属和 SEI 组分。此外，固态核磁共振（NMR）（7）和电子顺磁共振（EPR）（8）能够实现对 Li 镀层的时间分辨和定量检测。然而，这些方法大多是侵入性技术，需要额外昂贵以及具有特殊电池设计的精密仪器进行测试，在实际应用中是不可行的。 除了这些物理特性之外，电化学方法，如差分电压分析（DVA）、（9,10）差分充电电压（DCV）、（11）增量容量分析（ICA）、（12）和库仑效率（CE）、（13,14）也在电池老化研究中发挥了重要作用。然而，它们是相当复杂的确定方法，需要很长的监测时间。

电化学阻抗谱测试

• 为什么测出来的Nyquist图没有半圆部分呢？

理论知识

记得看Sob老师的博客，多逛逛计算化学公社呀

Comment and share

前言

1917 年，赫尔首次描述了简单粉末衍射仪的构造，这是在 1895 年威廉·康拉德·伦琴发现 X 射线后不久。衍射仪测量 X 射线反射的角度，从而获得其包含的结构信息。如今，这项技术的分辨率有了显著提高，它被广泛用作分析相信息和解决固态材料晶体结构的工具。X 射线衍射技术用于确定样品的组成或晶体结构。对于大分子和无机化合物等较大的晶体，它可用于确定样品中原子的结构。如果晶体尺寸太小，它可以确定样品的组成、结晶度和相纯度。

什么是X射线？

X射线是一种电磁波，由内壳层电子跃迁所得：

对于原子序数为Z、仅含一个电子的类氢原子，内壳层电子跃迁所辐射的能量可由Rydberg方程计算： $$ \bar{v}=(\dfrac{1}{n_i^2}-\dfrac{1}{n_f^2})RZ^2 $$ 可以明显看出，与电子跃迁相关的能量差随着原子序数的增加而大幅增大，并且在这种跃迁过程中发射的辐射波长随着原子序数的增加从 10⁻⁷ m (1000 Å) 范围移至 10⁻¹⁰ m (1 Å)范围（此时辐射的电磁波被定义为X射线，此辐射类型被称为特征辐射）、

如下图，要引发这种内壳层跃迁，需要产生一个电子空位，必须从K层中移除一个电子。由高压电场加速过的阴极电子书轰击靶材，将其部分能量传递给靶材料的电子，从而导致电子激发。如果入射电子的能量足够高，一些电子可能会击出靶中 K 壳层的一个电子，从而产生一个空位

不同能级跃迁回K层对应着不同的谱线，如$\rm{K_{\alpha}}$线和$\rm{K_{\beta}}$线。对于铜的X射线光谱，在低能分辨率下只有2条特征线被观察到，但是，在更高的分辨率下，很容易看出Kα线是双峰，它被标记为$\rm{K_{α1}}$和$\rm{K_{α2}}$。因为铜的2p轨道，能级L分裂为$\rm{L_{II}}$和$\rm{L_{III}}$，间距非常小（0.020kev），波长Kα1 (= 1.54056 Å) 和Kα2 (= 1.54439 Å) 也非常接近。因此，通常说的$\rm{K_{\alpha}}$其实是$\rm{K_{α1}}$和$\rm{K_{α2}}$的加权平均值

以铜靶为例 $$ E=\dfrac{hc}{\lambda}=\dfrac{6.626\times10^{-34}\times2.998\times 10^8}{1.542\times 10^{-10}\times1.6\times10^{-19}}=8.05~\rm{kV} $$ 8.05 kV为铜的$\rm{K_{\alpha 1}}$线的最低激发电压，实际操作时通常选择40kV。激发电压选择过低会造成谱线不明显、电流过大热损耗大等问题，选择40kV是长期实验优化的结果。

由于种种原因，如量子力学的选择定则：电子跃迁需要满足角动量变化Δl=±1、散热和成本等问题，通常X射线光谱中不会出现$\rm{K_{\gamma}}$甚至更高的线，即使能够发生N→K的跃迁（比如高Z元素W)，也会因为强度极弱而湮灭在噪声中，因此X射线光谱中通常见到的是$\rm{K_{\alpha}}$和$\rm{K_{\beta}}$

此外，在高速电子轰击靶材的过程中除了由内壳层电子跃迁所辐射的特征X射线，还存在着轫致辐射：高速电子被靶材原子核电场偏转，动能转化为连续电磁辐射。

如下图，当加速电压（单位：kV）较小、高速电子无法引起内壳层跃迁时（加速电压小于临界电压），那么此时原子的辐射均为轫致辐射

$$ V_{临界}=\dfrac{hc}{e\lambda_{K-edge}} $$

在所得的一系列X射线中，只有满足布拉格条件（波长与晶面间距d匹配）的射线才可用于样品检测，因此需要通过滤光片等手段只保留窄带波长，去除了除 K_α X 射线之外的所有电磁波

XRD谱图上的峰是什么？

XRD——X-ray diffraction（X射线衍射），其本质是波的干涉与叠加。

当波（如光波、声波、X射线）遇到障碍物或穿过周期性结构时，波前会发生弯曲和重新分布，这种现象称为衍射。障碍物或结构的尺寸与波长相当，这是发生衍射的关键条件，正因如此，与原子间距匹配的X射线适合用作发生衍射的电磁波。X射线是电磁波，照射到晶体时，每个原子中的电子受X射线的作用，成为次级波源，向四周散射X射线（这种散射也被称为相干散射），这些不同路径的X射线相遇并发生干涉（衍射发生的本质是弹性散射波，（即此时的次级X射线）的干涉）

干涉分为两类：

• 相长干涉（Constructive Interference）：满足布拉格条件时，所有晶面的散射波相位一致，两列波波峰对齐，叠加后振幅增强 → 形成峰。
• 相消干涉（Destructive Interference）：波峰与波谷对齐，叠加后振幅抵消 → 形成暗区。

布拉格定律： nλ = 2dsin θ 其中λ是外加波长，θ是衍射角，d 是原子平面之间的距离。在大多数情况下，我们假定n=1，即一级衍射，因为总可以将 n=2,3,…的衍射峰解释为来自(nh nk nl) 晶面的衍射——即来自晶面间距为$\rm{d_{hkl}}$的$\dfrac{1}{n}$倍的晶面$\rm{d_{nhnknl}}$。然后，原子平面之间的距离可用于确定成分或晶体结构。

当路径 ABC 与 A’B’C’之间的距离相差整数个波长（λ）时，衍射的 X 射线表现出相长干涉

X 射线衍射的结果绘制出了在各自的 2θ位置处不同衍射角的信号强度。2θ位置对应于样品中晶体或原子之间的特定间距，该间距由入射到样品中的 X 射线束的衍射角决定。峰的强度与该相或具有该间距的分子数量有关。峰的强度越大，具有该特定间距的晶体或分子数量就越多。

消光是什么？为什么只有特定的晶面处才会出现峰？

衍射峰强度取决于晶胞内所有原子散射辐射之间的相位关系，但经常出现这种情况：虽然布拉格定律预测某个峰应该存在，但其实际强度却为零（这是因为布拉格定律不涉及原子位置，仅与晶胞的大小和形状有关）。

例如，对于具有体心立方BCC晶胞的(100)晶面衍射峰强度：相位关系表明，在晶胞顶面和底面（即(100)晶面）散射的 X 射线虽然发生相长干涉，但与晶胞中心原子散射的X射线存在180^∘的相位差，因此最终强度为零。下表给出了不同立方晶格中特定衍射峰出现与否的选择定则。

为什么谱图的横坐标是2θ？

实验中，样品和探测器的转动是围绕同一轴进行的，因此 总偏转角的计算需要结合样品和探测器的相对运动，如下图所示

峰宽度的意义是什么？

峰的宽度与晶体尺寸成反比。较窄的峰对应较大的晶体。较宽的峰意味着可能存在较小的晶体、晶体结构缺陷，或者样品在本质上可能是非晶态的，即一种缺乏完美结晶度的固体。对于较小的样品，将由X 射线衍射（XRD）所得的图谱与数据库中标准样品的谱图对比可确定样品的组成。此结论可由谢乐公式（Debye-Scherrer）直接给出 $$ D=\dfrac{K\lambda}{B \cos\theta} $$

• K为Scherrer常数，也称形状因子，若B为衍射峰的半峰高宽，则K=0.89；

• D为晶粒垂直于晶面方向的平均厚度（Å）；

• B为实测样品衍射峰半高宽度

• θ为布拉格衍射角，单位为角度

• λ为X射线]波长，对于Cu kα来说，一般为1.54056 Å

为什么晶相能相长干涉形成峰，而非晶相不能？

晶体是长程有序的周期性结构，满足布拉格条件，即所有晶面(hkl)的的间距d恒定，当入射角θ满足布拉格公式时，所有晶面的散射波相位一致 → 相长干涉 → 尖锐峰。

非晶体可能短程有序，而长程无序，d值不固定， 散射波来自不同局部区域，相位杂乱，无法形成尖锐的峰

举个例子

例子仅供参考

晶体硅的晶胞参数a=0.543nm，X射线波长为0.154nm，那么我们可以预测某个晶面会在何处出现峰

对于立方晶系： $$ d_{hkl}=\dfrac{a}{\sqrt{h^2+k^2+l^2}} $$ 以(111)晶面为例 $$ d_{111}=\dfrac{0.543}{\sqrt{3}}=0.314nm $$ 根据布拉格公式： $$ \theta =\arcsin (\dfrac{\lambda}{2d_{111}})=14.3^\circ $$

如何确定2theta角度？

虽然理论上2θ范围是0°~180°，但实际操作中极少覆盖全部范围，原因包括：

• 物理限制：大多数晶面间距d有一个最小值，那么$\theta_{\rm{max}} \le \arcsin(\dfrac{\lambda}{2d_{\rm{min}}})$
• 信号强度：高角度区（如2θ>120°）的衍射强度极弱，信噪比差。
• 时间成本：全范围扫描耗时过长，而关键峰通常集中在20°~80°

一般选择20~80°，一般情况

常规晶体材料的晶面间距d大多在1~10Å之间

• Cu靶X射线波长λ≈1.54Å，代入公式：
  • 当d=10Å时，( $\sinθ = λ/(2d) ≈0.077$ → θ≈4.4° → 2θ≈8.8°)
  • 当d=1Å时，( $\sinθ = 0.77$ → θ≈50° → 2θ≈100°)

参考资料

1. Sadoway, D. Introduction to Solid State Chemistry.
2. Reig, A. CHEM322: Inorganic Chemistry.

Comment and share