使用VASP计算离子电导率，需要进行从头算分子动力学(AIMD)模拟，对象为$\rm Li_6PS_5Cl$
Li-Argyrodite型硫化物固态电解质。
收敛性测试
在进行昂贵的AIMD模拟之前需要进行收敛性测试，测试最优的ENCUT、k点等参数。如果设置太小则有失精度，设置太大又增加成本，收敛性测试的脚本网上应该能找到，我之前用Claude写了一个，发到了公社论坛上面：http://bbs.keinsci.com/thread-55151-1-1.html
变胞优化
一个能量较低、结构合理的起始三维构型是进行MD模拟的第一步，本次模拟使用变胞优化后的惯用晶胞进行。
AIMD
分子动力学模拟的基本步骤如下图：
这次的AIMD跑的很短，总共30ps(30000fs)，15000步，即每步2fs(POTIM = 2)。在30ps长的AIMD的模拟中，最初的一段时间体系并未达到平衡，因而需要舍弃，真正的AIMD有效数据此时间段向后截取。
以下是本次AIMD的INCAR:
1 | Global Parameters |
BTW，这里的SMASS参数可以设置0。我先进行了2000fs左右的预模拟，如下：
1 | [ctan@baifq-hpc141 SMASS]$ grep SMASS OUTCAR |
OUTCAR中给出的SMASS值是0.12，设置为0也可以。
简单介绍一下其他参数：
MDALGO设置为2，可以控制使用Nose-Hoover热浴TEBEG初始温度TEEND终止温度IBRION设置为-1启用分子动力学模拟
接着就开始漫长的模拟吧，大概模拟了八天左右：
1 | [ctan@baifq-hpc141 pre-pro]$ tail OUTCAR |
AIMD后处理
本次模拟的目的主要是为了得到离子电导率，模拟结束之后可以使用vaspkit方便地计算扩散系数，否则可能需要使用pymatgen自己编写代码。
首先使用下面的代码查看体系在合适趋于平衡，来自github上的一个老师，仓库地址是https://github.com/tamaswells/XDATCAR_toolkit
1 | #!/bin/python |
运行python程序：
1 | (base) [storm@X16 vaspkit]$ python 1.py |
查看温度和能量曲线，如下图
大概在5800fs时候体系达到平衡，于是可以使用5800fs的数据进行后处理，使用vaspkit的722功能
1 | ------------>> |
722功能可以得到MSD，均方位移，接着可以使用723功能得到扩散系数：
1 | 723 |
MSD数据是5800fs-15000fs的，所以723功能计算扩散系数的时候默认读取的这个区间，因此Please input the time steps to do linear fitting选择全部范围就行，最后可以看到锂离子的扩散系数是$\rm 3.555\times 10^{-6} cm^2/s$
根据Nernst-Einstein方程计算离子电导率： $$ \begin{align} \sigma&=\dfrac{ne^2Z^2}{k_BT}D_s\\ \end{align} $$ 其中的n为数密度，单位体积中的锂离子数量
kB为玻尔兹曼常数，Z为离子所带电荷，Ds为扩散系数，离子电导率可计算为 $$ \begin{align} \sigma&=\dfrac{\dfrac{24}{1076\times(10^{-8})^3}\times (1.602\times10^{-19})^2\times 1^2}{1.380\times10^{-23} \times 300}\times 3.555\times10^{-6}\\ &=0.49 \rm mS/cm \end{align} $$
计算机网络
第一章：概论
什么是Internet？
网络
计算机网络
互联网
节点
- 主机节点
- 画方框
- 数据交换节点：路由器（工作在网络层）、交换机（工作在链路层）、中继器、负载交换设备
- 既不是源，也不是目标，用作中转节点
- 画圆
链路：边，把各个节点连在一起
接入网链路(access)：主机连接到互联网的链路，出现“方框”
主干链路(bakcbone)：路由器（交换机）之间的，圆和圆连在一起
协议
- 物理层
- 网络层
从具体构成角度讲，Internet包括：
- 计算设备如主机
- 通信链路
- 光纤、同轴电缆无线电
- 传输速率：带宽(bps)
- 分组交换设备：用于转发分组
- 路由器
- 交换机
SEI膜的结构
锂离子电池工作电位范围为2～4.3V。其中，石墨类负极工作电位范围在0～1.0V（vs.Li+/Li），正极工作电位范围一般在2.5～ 4.3V（vs.Li+/Li）。而目前商用电解液不发生氧化还原反应的电化学窗口一般为1.2～3.7V（vs.Li+/Li）。因此，在负极侧，当负极点位 <1.2 V时会生成SEI固态电解质中间相，类似的在正极侧会生成CEI结构。
从机理来看，当负极的Fermi能级高于电解质的LUMO时，电子容易从Fermi能级转移到LUMO轨道，即发生氧化还原反应形成SEI界面
如果电极材料的充放电电位范围较窄，例如负极的嵌锂电位高于1.2V（ vs.Li+/Li ） ， 正极的脱锂电位低于 3.5V（vs.Li+/Li），则正负极表面可以不发生电解质的氧化还原反应，不会形成SEI膜。
SEI界相的厚度可能需要>2 nm以防止电子隧穿效应，<50 nm来确保离子能够顺利传输。它的结构如下图所示：
SEI膜的成分和结构通常认为是靠近电极材料的为无机物层，对于锂离子电池，主要包含$\rm Li_2CO_3、LiF、Li_2O$等；中间层为有机物层，包括$\rm ROCO_2Li、ROLi、RCOO_2Li$（R为有机基团）等；最外面为聚合物层，例如PEO-Li等
SEI的成分与微观结构基本不存在普适性的规律，针对特定的电极与电解质体系需要具体问题具体分析，且不 易表征。Novák等对SEI膜表征方法做了总结，包括：俄歇电子能谱（AES），飞行时间二次离子质谱（ToF-SIMS），扫描探针显微 镜（SPM），扫描电子显微镜（SEM），透射电子显微镜（TEM），红外吸收光谱（IRAS），拉曼光谱（RS），X射线衍射（XRD），电子能量损失谱（ELLS），X射线近边吸收谱（XANES），电化学阻抗谱（EIS），差分扫描量热仪（DSC），程序升温脱附仪（TPD），核磁共振仪（NMR），原子吸收光谱（AAS），电化学石英晶体微天平（EQCM），离子色谱（IC），二次电子聚焦离子束与元素线性扫描分析仪（FIB-ELSA），傅立叶变换红外光谱（FTIR），X射线光电子能谱（XPS）等
TEM适合探测纳米材料的微观结构，常用来研究电极材料的SEI膜形貌，但是SEI膜较敏感，在被电子束照射后会收缩甚至消失，在常规透射电镜下难以保持原有的化学状态，无法实现纳米尺度的原位观测。崔屹等人首次拍摄出了具有原子级分辨率的SEI膜透射电镜照片（冷冻电镜）。研究发现，具有马赛克结构SEI的金属锂脱锂不均匀，从而形成大量失 去与电极接触的金属锂，俗称“死锂”导致电池循环效率降低。与之相比，具有层状结构SEI的金属锂脱锂均匀，残留的“死锂”较少，所以循环效率也较高。
SEI膜的化学成分可以使用表面增强拉曼光谱（SERS）研究
| 项目 | 3C设备 | 动力电池 | 储能 |
|---|---|---|---|
| 质量能量密度 /(W·h·kg⁻¹) | 260~295 | 240~300 | 140~200 |
| 体积能量密度 /(W·h·L⁻¹) | 650~730 | 500~600 | 320~450 |
| 循环寿命 /周 | 1000 | 1500~3000 | 5000~15000 |
| 倍率 | 3C | 1~3C | 0.2~0.5C |
| 工作温度 /℃ | ‑20~55 | ‑30~55 | ‑30~55 |
| 成本 /(元·(W·h电芯)⁻¹) | 1.2~2.0 | 0.5~1.2 | 0.5~0.8 |
| 电池容量 /A·h | 3~20 | 3.2~200 | 50~400 |
| 工程能力指数 (Cpk) | 1.33~1.66 | — | — |
| 安全性（欧洲汽车危险等级） | 4 | — | — |
上表是三种常见的不同应用场合锂离子电池的性能参数
电化学池的环境可近似为恒温、恒压且非体积功只有电功，那么反应中吉布斯自由能的改变即为电功： ΔrGs = −nFEs markdown不方便打出化学标准状态，就用standard上标代替了。
n为每摩尔电极材料在氧化还原反应中转移电子的量，F为法拉第常数（F=96485C/mol），E为电化学池的电动势。
电池的能量密度可以用两种方式表示：质量能量密度（W·h/kg）和体积能量密度（W·h/L），根据不同的体系选择合适的标度，分别定义为： $$ \varepsilon_M=\dfrac{\Delta_rG^s}{\sum M} $$
$$ \varepsilon_V=\dfrac{\Delta_rG^s}{\sum V_M} $$
使用国际单位制，电量定义为： $$ Q=I\times t=\rm 1A\times 3600s=3600C $$ 所以： $$ \rm 1mAh=3.6C $$ 所以比容量定义为： $$ \text{Capacity}=\dfrac{Q}{m}=\dfrac{nF}{m}=\dfrac{F}{M}(\text{C/g})=\dfrac{F}{3.6M}(\text{g/mol}) $$ 常见的非理想情况：
- 反应物气体来自于外界，若在计算理论能量密度时不考虑气体的质量，计算出的理论能量密度会显著高于考虑气体质量的计算方法
- 计算采用的吉布斯生成能一般为不含缺陷的体材料（perfect bulk material）的测量数据，实际材料由于存 在缺陷和尺寸效应，导致生成能会偏离理想材料的生成能，因此需要考虑各类缺陷能的贡献：
ΔfGs(real material) = ΔfGs(perfect material) − ∑ΔfGis(defect)
在实际电池电芯中，存在多种非活性物质，如集流体、导电添加剂、黏结剂、隔膜、电解质溶液、引线、封装材料等。
终于也是把文献整理得差不多了，博客也翻新了一下，继续分享一些之前的计算。
我的计算体系是Argyrodite硫化物固态电解质体系，计算以$\rm Li_6PS_5Cl$为例吧。能带计算最好是使用primitive cell，即原胞。这是惯用晶胞或者超胞计算得到的能带是折叠的，固体物理的相关书籍中有提到，（当然计算成本也是不同的），参考公社论坛上的两篇讨论帖子：
- http://bbs.keinsci.com/thread-11373-1-1.html
- http://bbs.keinsci.com/thread-19361-1-1.html
原胞的变胞优化
计算使用我稍微熟一点的VASP来，首先从Materials
Project网站下载晶体结构文件，下载primitive
cell，如果已经有了conventional
cell的cif文件，也可以使用vaspkit的602功能生成PRIMCELL.vasp原胞文件。
然后进行原胞的变胞优化，除了POSCAR的剩下三个输入文件可以使用vaspkit方便地生成，使用101功能，选择lattice relaxation进行变胞优化，102功能生成自洽计算的K点，选择Gamma撒点方案，0.03间隙，原胞计算量比较小，间隙填小一点也无妨。
笔者在Intel
5218R的集群上进行计算，分配了40个物理核心，计算时在INCAR中设置KPAR=2和NCORE=20来最大化效率，一般而言NPAR参数不需要设置，KPAR和NCORE的参数设置众说纷纭，我参考的公社某一篇帖子里的经验，前者设置为2，后者设置为物理核心数的一般，当然这是对于CPU版本的VASP而言的。
原胞的变胞优化用时很短：
1 | reached required accuracy - stopping structural energy minimisation |
在确认收敛并确认结构正常后进行能带的计算
PBE泛函计算能带
能带计算的KPOINTS与普通计算的KPOINTS不一样，通常需要第一布里渊区内的一条或几条高对称点路径来计算能带性质
VASP的计算中PBE泛函是被使用最广泛的，需要修改K点和输入文件。参考vaspkit官网的教程：http://vaspkit.cn/index.php/27.html
PBE能带的计算可以分两步进行。
首先进行自洽计算，此时也即（最常用的）静态计算，将变胞优化的CONTCAR更名为POSCAR，使用vaspkit的101功能，选择Static-Calculation，k点选择可以不变。
接着进行能带计算，需要使用vaspkit的303功能生成高对称K点，将生成的KPATH.in文件更名为KPOINTS
1 | (base) [storm@X16 temp]$ ls |
可以看到，303功能也生成了原胞结构，当然前面我们已经提到过。
在进行能带计算的目录中，将上一步自洽计算的CHGCAR文件拷贝过来并在INCAR中设置电荷密度不更新，即ICHARG=11，这是必须的，为了加速计算，建议读取上一步的WAVECAR文件，虽然ISTART设置为1或者2程序都可以读取WAVECAR，但是为了保持整个计算任务中（第一步和第二步）的平面波基组一致，建议使用ISTART=2。
此种方法计算能带分两步的原因是使用高密度 k
点网格进行计算时耗时较长，因此用不耗时的普通静态计算多次迭代得到的CHGCAR和WAVECAR来加速高密度k点网格的计算。
其实对于我的小体系，直接使用高密度k点网格进行静态计算也是完全可接受的，并且SCF 和能带完全自洽，会更精确一些，用时并不长：
1 | [ctan@baifq-hpc141 primitive_cell-band_PBE]$ ls |
HSE06杂化泛函计算能带
PBE最常被诟病的方面即能带计算，它常常低估带隙。By the way，Materials Project上的能带数据刚好能和我用PBE算出来的能带数据对上，明显小于HSE06泛函计算的数据，估计是使用PBE算的。
HSE (Heyd–Scuseria–Ernzerhof)杂化泛函在2003年首次被提出（HSE03），在2006年HSE06被提出，它能更加精确地计算能带。使用HSE06杂化泛函做计算是十分昂贵的，用时是PBE泛函的数十倍不止。
受vasp算法所限制，使用HSE06杂化泛函计算能带时只能使用自洽计算（此时为静态计算），不能进行非自洽计算。中文互联网上有一些教程是错误的，PBE可以自洽+非自洽两步，也可以自洽一步，但是HSE06杂化泛函只能一步自洽。不过可以在能带计算之前也先做一下静态计算（自洽计算），然后将产生的WAVECAR文件用于能带计算以加速。
使用vaspkit的251功能生成所需要的KPOINTS文件，这是一个“混合”的文件：
1 | ------------>> |
首先需要生成自洽（静态）计算所需要的k点网格，然后生成能带计算的高密度k点网格。INCAR文件可以沿用普通静态计算的然后改一下泛函，也可以使用101功能生成HSE06
Calculation的INCAR。
HSE06杂化泛函极为昂贵，大概在第五个电子步时，用时超过8h仍未计算完，遂放弃。之前注册某家超算平台时送了一点机时，开通了GPU版本的vasp计算，发现快了不少：
1 | General timing and accounting informations for this job: |
数据后处理与Origin绘图
能带计算结束后，将目录拷贝到可以使用vaspkit的机器上。
对于PBE计算的能带，使用211功能处理数据，得到带隙、CBM、VBM、Fermi Energy等信息：
1 | ------------>> |
对于HSE06计算的能带，使用252功能处理数据，得到带隙、CBM、VBM、Fermi Energy等信息：
1 | ------------>> |
绘图时我习惯把价带顶归零，即纵坐标为$\rm E-E_{VBM}$，操作如下。
把FERMI_ENERGY文件就拷贝一份，改名为FERMI_ENERGY.in。将文件FERMI_ENERGY.in中的费米能级数据改为价带顶数据，然后再重新运行252功能，此时得到的REFORMATTED_BAND.dat文件中的数据是归零的。
接下来就是origin绘图了
将REFORMATTED_BAND.dat文件拖入工作表内，确认只有一个X轴，剩下全是Y轴。新建一个sheet。打开KLABELS文本文件，手动复制数据，如下图
接着选中sheet1表格中的所有数据，绘制折线图：
设置一下横坐标，改成高对称点的样子，看图：
这样横坐标就被设置好了：
接下来就是美化了，放一张最后画好的图吧：
费米孔
Coulomb孔
由于Coulomb孔的存在，电子之间存在动态相关，如何描述原子分子体系中电子之间的动态相关是量子化学的主要任务
玻色子的波函数是对称的，没有反对称性的要求因此可以靠得很近
Hartree-Fock方程与自洽场计算
HF方程
HF方程是量子化学中计算多电子体系近似波函数和能量的基本方法，通过“平均场”近似把复杂的多电子问题拆解成单电子问题。
多电子波函数用单电子轨道（分子轨道）的Slater行列式表示 f̂(xi)φk(x1) = εkφk(x1) 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方程的迭代计算过程，通过不断更新“平均场”直到结果自洽。
关键步骤
猜初始轨道（如用原子轨道线性组合）。
构建Fock算符 （依赖当前轨道）。
解HF方程 ，得到新轨道。
检查收敛
新轨道与旧轨道是否几乎相同？
- 是 → 计算结束，输出能量和波函数。
- 否 → 用新轨道回到步骤2，继续迭代。
为什么叫“自洽”？
最终得到的轨道生成的“平均场”，与计算这些轨道时用的场是一致的（自圆其说）。
1 | [多电子薛定谔方程] |
组态相互作用（Configuration Interaction, CI）计算
CI是一种后Hartree-Fock方法，通过将多个电子组态（Slater行列式）线性组合，考虑电子相关能，弥补HF方法的不足。
核心思想
HF的缺陷： HF单行列式近似忽略了电子瞬时排斥（动态相关）和激发态贡献（静态相关），导致能量偏高。
CI的改进： 将波函数表示为多个激发组态的叠加： ΨCI = c0ΦHF + ∑iciΦisingle + ∑jcjΦjdouble + ⋯
- (Φisingle)：单激发组态（一个电子从占据轨道跃迁到虚轨道）。
- ( Φjdouble )：双激发组态（两个电子同时跃迁），贡献最大。
CI等级
| 类型 | 包含的激发组态 | 计算成本 | 典型应用 |
|---|---|---|---|
| CIS | 仅单激发 | 低 | 激发态初步探索 |
| CISD | 单+双激发 | 中等 | 小分子基态相关能修正 |
| CISDT | 单+双+三激发 | 高 | 高精度小体系 |
| FCI | 全组态（所有可能激发） | 天文数字 | 极限精度（如H₂O极小基） |
关键问题：
- 计算量爆炸：FCI的组态数随电子数和轨道数阶乘增长，仅适用于极小体系。
- 大小一致性：CISD不满足大小一致性（如两个远离的H原子能量≠2×单个H能量），需更高激发或耦合簇（CC）修正。
1 | [Gauss基函数] |
理论计算中的各种方法
简单总结了一下Sob老师的博文，个人的总结仅供参考，详细信息请阅读原文：http://sobereva.com/680
任务类型
- 优化极小点
- 优化过渡态
- 产生反应路径
- 振动分析
- 分子动力学
- 构象/构型搜索
- 波函数分析
计算方法
| 计算方法 | VASP | |
|---|---|---|
| 分子力场 | 计算耗时极低 | 无法描述化学反应 |
| 机器学习势 | 通过机器学习的思想构造依赖于原子坐标的分子描述符与能量之间的抽象关系 | |
| 密度泛函理论（QM） | 性价比非常高 | 耗时高于分子力场N个数量级 |
| Hartree-Fock（QM） | 完全过时 | |
| 后Hatree-Fock类（QM） | 额外把动态相关在一定程度上考虑了进来 | 耗时高了很多 |
| MCSCF如CASSCF（QM） | 弥补HF缺乏对静态相关的考虑 | 几乎没有或很少考虑动态相关 |
| 多参考方法如CASPT2、NEVPT2、MRCI（QM） | 在MCSCF基础上进一步把动态相关考虑进来，精度整体很好，普适性很强 | 昂贵 |
| 半经验类方法如AM1、PM3、PM6（QM） | 对HF的简化以巨幅降低耗时，耗时只是HF的微小零头 | 精度低，只能用于有限的元素 |
| GFN-xTB（QM） | 相当于纳入了一定DFT思想的半经验级别的方法，整体精度和可靠性>=主流的半经验方法 | |
凡是基于量子理论思想提出的，在实际数值求解的过程中一般都要涉及分子轨道，绝大多数计算程序中都是把分子轨道展开成基函数的线性组合来描述的
最常见的基函数一类是原子中心基函数（如Gaussian等大多数量子化学程序以及CP2K等部分第一性原理程序用的高斯型基函数），其中心一般位于原子核，还有一类常见的基函数是平面波基函数，是大多数计算周期性体系为主的第一性原理程序用的，它的分布覆盖整个被计算的晶胞
基组（basis set）是对于原子中心基函数而言的，例如6-31G*、def2-TZVP、cc-pVDZ等都是很常见的基组，它定义了实际计算时对各种元素原子具体用多少、什么参数的基函数。做HF、DFT、后HF、MCSCF、多参考等方法计算时都需要告诉计算程序用什么基组，基组质量越好，也就是越接近于所谓的完备基组极限，这些方法自身的精度发挥得就越充分，但代价就是耗时越高。一个好方法配一个烂基组，以及一个烂方法配一个好基组，结果都不理想，必须好方法配好基组才能得到较准确结果。
仅供学习、交流。
好久没来博客了，分享一下近期折腾的东西以及踩的坑吧。
之前尝试过使用红米AC2100路由器实现校园网认证，但是没有记录，并且其处理MT 7621性能孱弱，干不了太多事，这次换成了稍微热门一点的——捷稀 JCG Q30 Pro，处理器是MT 7981B,架构为Arm Cortex-A53 (1.3 GHz, dual-core)，可以干的事更多了，不过个人还是主要用于dogcom认证JLU的校园网以及搭建openclash实现透明代理。
安装ImmortalWRT
首先是替换路由器官方的固件，捷稀JCG Q30 Pro似乎是中国移动的固件，可以为其刷入OpenWRT，也可刷入衍生版比如Immortalwrt、QWRT等等，说来话长，就不提供教程了，具体教程在恩山论坛。
SDK编译
路由器固件用的是Immortalwrt，爱来自恩山大佬😀，原帖链接：https://www.right.com.cn/forum/forum.php?mod=viewthread&tid=8398454&highlight=%E6%8D%B7%E7%A8%80%2BJCG%2BQ30%2BPro
泥吉校园网使用的是哆点客户端，即Dr.COM认证(包括DHCP、PPPoE、802.1x三种认证方式，即D版、P版、X版)，桌面端(windows、Linux、MacOS)下载客户端认证，然而Dr.COM备受吐槽、嫌弃，Arm路由器也无法使用Dr.COM。各位民间大佬分析了DrCOM的通信协议，开发了drcom-generic。当时的开发者是通过python2脚本来代替Dr.COM实现校园网认证的，下图是项目仓库。
然而，现在python3才是主流，OpenWRT及其衍生版本的官方仓库中早已抛弃了python2，在离线环境安装python2还是比较痛苦的。因此，后续大佬使用C语言重新复现drcom-generic的功能，将新项目命名为dogcom，以及后续有大佬开发了泥吉专属的C语言版客户端。
本次使用前者，开发者大佬开源了源代码，需要将源代码编译为MT 7981可用的二进制文件，由于MT 7981性能远不如AMD64处理器，因此编译在Linux(WSL)上完成。笔者使用的机子为2*2676V3，固件为Immortalwrt，首先搭建交叉编译环境。
交叉编译环境需要固件的SDK，然而Immortalwrt并没有现成的关于MT 7981的SDK，还得先手动编译SDK😅。进入适用于处理器的Immortalwrt项目仓库，如下图。
一定一定一定听劝，使用推荐的Ubuntu 20.04 LTS。高版本的gcc编译套件可能会报错，Arch系已踩过坑，降级gcc版本够折腾的😂。以及Rocky Linux 10可能会遇到ninja编译报错。
首先更新源，然后更新系统，接着安装所需依赖。
1 | sudo apt update -y |
然后就可以克隆代码，开始搭建交叉环境了。
1 | git clone --depth=1 https://github.com/hanwckf/immortalwrt-mt798x.git |
如果是Win11的话，可以在WSL设置中将网络模式设置为镜像mirror，然后在Windows上开启代理，那么WSL中就能连上GitHub。
接着cp -f defconfig/mt7981-ax3000.config .config拷贝配置文件，然后make menuconfig进一步编辑，如下图。
确保上面三项正确。
我们只需要编译一下SDK(不需要生成tar.xz文件)搭建交叉环境就可以了，不需要编译整个Immortalwrt for MT 798x因此执行make toolchain/install -j$(nproc)而不是make -j$(nproc)。
dogcom编译
交叉编译环境搭建好后，开始编译dogcom.
在Immortalwrt根目录执行:
1 | git clone https://github.com/mchome/openwrt-dogcom.git package/openwrt-dogcom |
1 | git clone https://github.com/mchome/luci-app-dogcom.git package/luci-app-dogcom |
前者为dogcom本体，后者为图形化界面。然后再次:
1 | make menuconfig |
分别添加dogcom和luci-app-dogcom
然后save，重新生成.config配置文件，然后就可以开始编译这两个软件包了（dogcom及其luci图形界面）。
执行：
1 | make package/openwrt-dogcom/compile |
这两个软件包都不大，很快就编译好了，生成的ipk文件（适用于OpenWRT及其衍生发行版的安装包文件）在/bin/packages目录下：
安装dogcom及图形化界面
拿到ipk文件后就可以安装了。使用ssh登录Immortalwrt。使用winscp上传两个文件，OpenWRT及其衍生版的包管理器为opkg。
配置校园网
具体的设置可以参考：https://www.bilibili.com/opus/1063006213972688900#reply273467151024
比较麻烦的是配置静态IPV4地址，很容易配置错误导致连不上认证服务器。
参考资料
https://www.bilibili.com/opus/1063006213972688900#reply273467151024
https://shenshichao.notion.site/OpenWrt-JLU-fb5132d707114e22a4486b0e657421f0
https://www.right.com.cn/forum/thread-215978-1-1.html
https://github.com/mchome/dogcom?tab=readme-ov-file
https://deconf.xyz/2019/04/13/openwrt-dogcom/#0x03-%E5%AE%89%E8%A3%85dogcom%E5%AE%89%E8%A3%85%E5%8C%85-%E5%90%AF%E5%8A%A8%E8%AE%A4%E8%AF%81%E7%A8%8B%E5%BA%8F
使用最新版的Miniconda3，python环境的版本是14.x，顺手pip install numpy安装的numpy版本过高，使用vaspkit尝试将.cif文件转换为POSCAR时报错:
1 | AttributeError: `np.mat` was removed in the NumPy 2.0 release. Use `np.asmatrix` instead |
错误原因在于当前的NumPy版本太新（>=2.0），而vaspkit的cif2pos.py脚本中使用了已被移除的np.mat函数，最简单的方式是降级numpy版本：
1 | pip install numpy==1.26.4 |
pi-pi作用的定义和说法比较乱，Sob老师认为：pi-pi作用是相距较近的两个片段上彼此朝向相对的pi电子之间的独特的色散吸引作用。
- pi-pi作用可以是分子间的也可以是分子内的，满足以上定义即可
- 诸如苯分子里面的pi电子之间的作用不叫pi-pi作用，那属于共享电子作用
- 两套pi电子的分布必须近乎彼此相对才可能算pi-pi作用。而诸如两套pi电子近乎肩并肩挨着就不能算pi-pi作用，像是不能说环丁二烯里面两套近乎定域的pi电子之间是分子内pi-pi作用
- “相距较近”一般是在4.0-4.5埃以内。也不是说更远距离就完全没有pi-pi作用了，只不过由于色散作用随作用距离r呈$-\dfrac{1}{r^6}$快速衰减行为，因此距离稍微一远pi-pi作用就非常弱了，就不太值得一提了。但相互作用的片段间距离也不能太近，若显著小于相接触的原子的范德华半径之和，则显著的位阻互斥作用会远大于pi-pi吸引作用，使得pi-pi作用也相对不值得一提
- pi-pi作用最常出现在碳原子间，因为碳最容易带显著的pi电子。碳的Bondi和CSD范德华半径分别为1.70和1.77埃。如果没有特殊的因素影响相互作用距离的话，在平衡结构（势能面极小点结构）下，C-C间的pi-pi作用出现在3.4-3.6埃左右是最常见的
- pi-pi作用属于弱相互作用和非共价相互作用范畴。虽然名义上算作弱相互作用，但实际强度可大可小，直接取决于作用面积，详见后文
- 有一个词叫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堆积作用的本质都是色散作用。
注：
- 还有个作用叫交换互斥，它起到互斥作用，不管算什么类型相互作用都要考虑这个，==只有在较近距离（比如小于两个原子的范德华半径和）的时候才开始凸显出来，且原子间距离越近此作用越大==。正是这个作用使得弱相互作用的势能曲线总是有极小点，而不会令原子间距离因为静电吸引和色散作用而无限减小。
- “范德华作用”这个词描述的是电中性的两个原子之间的非共价相互作用，故本质包含了交换互斥和色散吸引作用两个部分。这个词可以囊括比如两个氩原子间的相互作用的全部内涵。这个词也可以用来描述无极性分子间的全部相互作用（但实际上，没有绝对严格意义的无极性分子，比如就连一般说的无极性分子氮气也有四极矩，可以认为是局部极性，因此氮气之间也有静电相互作用，讨论见http://sobereva.com/209。但本文就不咬文嚼字了）。至于一些初级的教材里说的范德华作用包括“诱导力、色散力、取向力”，这种歪曲的说法应该彻底从量化研究者的脑海中删除
- “范德华作用”这个词描述的是电中性的两个原子之间的非共价相互作用，故本质包含了交换互斥和色散吸引作用两个部分。这个词可以囊括比如两个氩原子间的相互作用的全部内涵。这个词也可以用来描述无极性分子间的全部相互作用（但实际上，没有绝对严格意义的无极性分子，比如就连一般说的无极性分子氮气也有四极矩，可以认为是局部极性，因此氮气之间也有静电相互作用，讨论见http://sobereva.com/209。但本文就不咬文嚼字了）。至于一些初级的教材里说的范德华作用包括“诱导力、色散力、取向力”，这种歪曲的说法应该彻底从量化研究者的脑海中删除
- 《透彻认识氢键本质、简单可靠地估计氢键强度：一篇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}} $$ 其中ρ是电子密度，可有量子化学程序计算出，也可由高分辨晶体衍射确定
| Around nuclei | Around chemical bonds | Weak interaction regions | Boundary of molecule | |
|---|---|---|---|---|
| $\lvert \nabla \rho(\bold r) \rvert$ | Large | 0 ~ Minor | 0 ~ Small | Very small ~ Small |
| ρ(r) | Large | Medium | Small | 0 ~ Small |
| RDG($\bold r$) | Medium | 0 ~ Minor | 0 ~ Medium | Medium ~ Very large |
在原子核附近，低能级上的电子高速运动，电子密度大，周围的电子密度变化迅速（即$\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)$与λ2符号的乘积，其中λ2代表$\rho (\bold r)$的Hessian矩阵的第二大特征值。用于在RDG等值面上映射$\operatorname{sign}\left( {\lambda }_{2}\right) \rho (\bold r)$的常用色标如下图所示。由于化学体系中相互作用的复杂性和多样性，图中的标注仅适用于一般情况。
范德华(vdW)相互作用区域的ρ通常非常小，因此图2中对应的sign (λ2)ρ颜色为绿色。此外，π − π堆叠具有与范德华吸引部分相同的物理性质，即色散效应(dispersion effect)。因此，π − π堆叠区域的ρ也相当小。对应明显的空间位阻效应或吸引性弱相互作用的区域(例如，中等强度的H-键和X-键)具有较大的ρ(通常为0.02-0.05原子单位)，因此在图2中分别对应红色和蓝色区域。虽然离子键的强度通常比弱相互作用大一个数量级，==但它们相互作用区域内ρ的大小是相当的==。例如，Na-Cl键的BCP(临界点)处的ρ仅约为0.03原子单位。因此，RDG等值面不仅可以揭示弱相互作用，还能显示离子键相互作用。相反，具有显著共价性的相互作用，尤其是共价键，==其在相互作用区域的ρ通常明显大于0.05原子单位==。因此，它们在采用密度截断的标准NCI(非共价相互作用)图中不可见。显然，NCI本质上是一种旨在揭示各种非共价相互作用的方法
在绘制NCI(非共价相互作用)图之前，波函数分析代码需要计算大量在矩形区域内均匀分布的网格点的RDG和sign (λ2)ρ值，区域应合理设置以确保所有感兴趣的RDG等值面都位于其中。可以绘制RDG与sign (λ2)ρ的==散点图==，以更好地理解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 (λ2)ρ的网格数据质量非常敏感，这主要体现在网格点之间的间距上。网格间距过大会导致等值面出现难看的锯齿边缘；但如果网格间距过小，则会带来过高的计算成本。请注意，如果空间范围固定，计算网格数据所花费的时间与网格间距的三次方成反比。根据经验，网格间距为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) $$
- 物理意义
- 正曲率（凸起）→ 电子密度在某个方向“向外膨胀”（如原子核附近）
- 负曲率（凹陷）→ 电子密度在某个方向“向内收缩”（如化学键或分子间区域）
λ2是其Hessian矩阵的第二本征值，在AIM理论中：
- 键临界点处sign(λ2)=-1
- 环、笼临界点处sign(λ2)=+1
- 接近临界点的区域其值与临界点处一般相同
可以将sign(λ2)函数用不同色彩投影到RDG等值面上，用来表现某一个区域的相互作用类型，==这便是用RDG方法分析弱相互作用==，但不如==IRI方法==全面，IRI与RDG在展现弱相互作用方面效果是等同的，而IRI有个关键性的好处是可以把化学键作用区域一起直观地展现出来！因此可以一目了然地图形化考察体系当中存在的所有类型相互作用，更有价值。考察体系的弱相互作用时会绘制电等值面图，IRI方法绘制的等值面图横坐标是sign(λ2)ρ(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等值面上，用来表现某一个区域的相互作用类型。
