实验部分

球磨

• 球磨的原理是什么？
• 球磨的作用是什么？
• 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

什么是电子能谱？

电子能谱就是用一定能量的==电子束或光子与样品表面相互作用==，使样品表面原子中不同能级的电子激发成自由电子，并研究这些自由电子的电子强度（电子数目）按其能量变化的分布曲线。根据使用的激发源的不同，电子能谱又分为==X射线光电子能谱（XPS）、紫外光电子能谱（UPS）、俄歇电子能谱（AES）==

电子能谱的基础是光电效应：当能量为hν的单色光照射到样品上时，样品中的原子或分子吸收能量从而使电子脱离样品成为自由电子（即光电子）逸出，与此同时还能产生==俄歇电子和荧光X射线等== hν = E_b + E_k + E_r 其中E_b为原子的始态和终态能量差，可以看成发射的光电子的结合能；E_r为原子的反冲能量，可忽略

X射线光电子能谱

X射线光电子能谱 （XPS），也称为化学分析电子能谱 （ESCA）。价电子对X光子的光电效应截面远小于内层电子，所以XPS主要用于研究内层电子的结合能。由于内层电子不参与成键，保留了较多的原子性质，所以常用于表面组分的化学分析，也可以进一步应用于确定这些==元素的化学或电子状态==。它允许以非破坏性方式测定样品的原子组成，由于X射线仅能穿透样品 5 – 20 Å，因此可以进行表面特异性分析。

电子在暴露于足够能量的电磁辐射（波长在X射线范围的高能光子）时会从材料表面射出，即光电效应： KE = hν − E_b − φ 其中E_b是电离能（也叫电子结合能），φ 是光谱仪的功函数。

不同的原子具有不同的电子壳层和不同的结合能，因此XPS可用于识别样品表面的元素组成。计算发射的具有各种E_k的光电子的数量，并生成XPS光谱，这有助于识别除H和He之外的所有状态的所有元素。

如上图所示，X射线光电子能谱仪器的主要部件包括以下部分：

• X射线源
• 超高真空（UHV）
• 电子透镜
• XPS能量分析仪
• 检测器
• 监视器

首先，高速电子撞击金属阳极使其发生内壳层跃迁，辐射出特征X射线。特征X射线脱离射线源部分后进入到UHV超高真空腔，在这里，约10⁻⁹~10⁻¹⁰ torr甚至更低的压强可以减少XPS光谱中的噪声量并避免样品表面污染产生任何额外峰。接着，X射线在UHV环境中照射到样品表面，样品释放出光电子，光电子被电子透镜所汇聚并减速。大量减速后的光电子接着进入XPS能量分析仪，这是一个两端施加有电压的半球腔体。多通道检测器可以检测出电子的动能，类似于胶片。高性能计算机用作监视器，用于控制镜头和分析仪的电压，它还会产生 XPS 光谱——计算机自动绘制E_k强度与结合能的关系生成的。

定性分析

原子的电子能级可分为两种类型——与原子核结合的核能级和仅弱结合的价能级。原子的价能级是与其他原子的价能级相互作用以在分子和化合物中形成化学键的能级。它们的特性（价键特性）和能量通过这个过程发生变化，成为形成的新物种的特性。

原子的核心能级电子的能量几乎与原子结合的化学物质无关，因为它们不参与键合过程。因此，核心水平结合能的识别因此提供了元素的独特特征。元素周期表中的所有元素都可以通过这种方式识别（除了H和He，它们没有核心能级）。

定量分析

定量分析需要测量相对峰强度，峰的强度取决于表面膜中存在的元素的原子浓度 （N）。

化学状态分析

化学位移是同一分子或两种不同化合物中两种化学不同状态的同一原子之间的结合能差异。它是样品电离原子周围的原子性质和原子数的特征。它可以用来查看样品表面原子的化学键是什么样的。

• 它涉及化学键性质、官能团的性质和原子的氧化态。
• XPS 的一个重要功能是它能够区分不同化学状态的相同元素，因为电子的 BE 会因氧化引起的元素周围电荷密度的变化而不同

深度分析

对薄表面薄膜中元素深度分布的估计称为深度表面分析。这样做是为了知道薄膜是否均匀。

优势与局限性

优势：

• 可以分析所有样品，包括固体、气体和液体样品，但大多数是固体样品。
• 定量效果好
• 出色的化学状态测定能力
• 它适用于从生物材料到金属的各种材料。
• 这是一种非破坏性技术。

局限性：

• 缺乏良好的空间分辨率，即空间分辨率差
• 仅中等绝对灵敏度
• 它无法检测氢和氦
• 它不适用于痕量分析。

为什么常常用Al作为靶材

X射线衍射（XRD）实验中选择的是铜靶或者钼靶，为什么这里用的是铝靶呢？

XPS对X射线源的光子能量和能量分辨率（单色性/线宽）有非常特殊的要求。XPS的核心是精确测量样品表面原子发射的光电子的动能，进而计算其结合能。结合能的微小偏移（< 1 eV）往往对应着元素的化学态变化（如氧化态、化学键）。==如果X射线本身的线宽很宽，就会导致光电子谱峰变宽，使得相邻的、化学位移很小的峰无法分辨，严重影响化学态分析的精度==

XRD测量的是晶面间距引起的衍射峰位置（角度）。衍射峰的自然宽度通常比X射线源的线宽要宽得多（受晶体尺寸、应力等因素影响）。因此，X射线源本身的线宽（如Cu Kα约2.5 eV）对衍射峰位置测量的影响相对较小，强度更重要。

Al Kα的能量是1486.6 eV，这个能量足够高，能够有效地激发周期表中绝大多数元素（从轻元素如Li到重元素）的核心能级电子（如C 1s, O 1s, N 1s, Fe 2p, Au 4f等），这些能级是XPS分析化学态的关键。能量不是越高越好。能量过高（如Cu Kα的8048 eV）会产生更强的韧致辐射本底和更多的二次电子，增加谱图的背景噪音。同时，高能光子的光电离截面（电离效率）在特定能量范围外会降低。1486.6 eV是一个在激发效率和背景噪音之间取得良好平衡的能量点。

为什么光电方程中有功函数这一项？

什么是功函数呢？功函数 Φ 是指将一个电子从==固体材料==内部（==费米能级处==）移动到真空中静止状态（真空能级）所需的最小能量。想象电子在固体内部处于一个“能量坑”的底部（费米能级）。要把这个电子从坑里拉出来扔到坑外平坦的地面上（真空能级），你需要克服坑壁的阻挡做功。这个“坑壁的高度”就是功函数 Φ（单位：电子伏特）。

为什么光电方程中有功函数这一项？理解的关键在于区分 ==“原子内结合能 (Eb)”==和==“电子逸出样品进入真空并被检测器测到”== 这两个过程所涉及的能量参考点。

Einstein的原始光电方程是针对于自由原子或者气体的，Eb 是电子相对于原子真空能级的结合能（即把电子从原子轨道移到离原子无穷远处所需能量）：

• 对于孤立原子，真空能级是明确的参考点
• 测得的光电子动能直接反应原子内电子的结合能E_b

但是对于==固体==，情况变得更加复杂了。固体中，原子不是孤立的，电子处于==周期性势场中==。固体中的电子能量通常参考费米能级 (Fermi Level, EF)——绝对零度下电子占据的最高能级。总结一下固体中的光电效应：内壳层电子吸收光子能量跃迁到无穷远轨道处（依然处在固体体相中），这一步需要克服结合能，从无穷远轨道处跃迁到真空层（脱离表面）还需要克服一部分能量，这部分能量被称为功函数φ。

所以功函数补偿了谱仪测量参考点（其费米能级）与理论计算参考点（真空能级）之间的差异。它确保了测量动能 KE 与样品内结合能 Eb（相对于样品费米能级）之间的定量关系成立。

XPS怎么确定元素的氧化态？

1
2
3
4
5
6

graph TD
A[氧化态升高] --> B[原子周围电荷密度降低]
B --> C[原子实有效正电荷增加]
C --> D[内壳层电子受核束缚增强]
D --> E[结合能E_b升高]
E --> F[XPS峰向高结合能方向移动]

XPS的测试项目

• 全谱
  • 得出元素的成分信息，也可以对所测元素进行半定量计算
• 精细谱
  • 判断化学态，确定不同化学态的百分比含量
• 刻蚀/溅射
  • 去除样品表面杂志后进行分析
  • 深度剖析，探测不同深度下的样品成分
• 价带谱
  • 主要测量价带顶，也可以分析价带结构，探测外层价电子信息，也可以用来区分化学环境价态
• 俄歇谱
  • 主要测量价带顶，也可以分析价带结构，探测外层价电子信息，也可以用来区分化学环境价态
• Mapping
  • 表征不同元素成分在特定区域的分布情况，分辨率〈EDS，XPS深度10nm左右，EDS深度1um左右

样品制备

• 粉末
  • 样品量20mg以上，压片制样，粘到双面胶带上分析；颗粒越细越好，且分散均匀；在压片样品表面分析时尽量选用平整测试区域，越平整信号越强

Comment and share

本文仅为学习记录，并非教程

Multiwfn是一个波函数分析程序，功能非常强大，也能绘制ELF等值面图，它支持多种类型的输入文件，既有量子化学计算程序的输出文件，又有第一性原理计算程序的输出文件

以下是Multiwfn加载ELFCAR绘制等值面图的过程

• step1：用ELFCAR绘图其实是基于格点数据文件的绘图，需要修改根目录下的配置文件settings.ini，将参数iuserfunc改为-1或者-3
• step2：加载ELFCAR文件
• step3：4 Output and plot specific property in a plane
• step4：100 User-defined function (iuserfunc=
• step5：5 Shaded relief map with projection
• step6：1: XY plane
• step7：1a

如果将Multiwfn的根目录添加到PATH环境变量以方便调用，那么画出来的图可能会遇到如下如题

如果直接在根目录下双击.exe文件运行的话就能正常绘图

Comment and share

起因是我想合成各式各样的材料，需要用到各种药品，计算起来比较麻烦，尤其是非化学计量比的化合物比如Li_5.3PS_4.3Cl_0.7Br，遂用AI写了一个python脚本，个人觉得还是比较实用，脚本文件和相对原子质量的文本文件需要在同一文件夹内

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

from decimal import Decimal, getcontext
import re
from sympy import symbols, Matrix, linsolve, Rational
from collections import defaultdict

# Improve Decimal precision
getcontext().prec = 100

class FormulaParser:
    @staticmethod
    def parse(formula: str) -> dict:
        elements = {}
        pattern = re.compile(r"([A-Z][a-z]?)([\d.]*)?")
        idx = 0
        while idx < len(formula):
            m = pattern.match(formula, idx)
            if not m:
                raise ValueError(f"Invalid formula at: {formula[idx:]}")
            elem, num = m.group(1), m.group(2)
            coeff = Decimal(num) if num else Decimal(1)
            elements[elem] = elements.get(elem, Decimal(0)) + coeff
            idx = m.end()
        return elements

class AtomicMassReader:
    @staticmethod
    def read(filename="atomic_masses.txt") -> dict:
        masses = {}
        try:
            with open(filename, 'r') as f:
                for line in f:
                    line = line.strip()
                    if not line or line.startswith('#'):
                        continue
                    parts = line.split()
                    if len(parts) != 2:
                        raise ValueError(f"Line format error: {line}")
                    masses[parts[0]] = Decimal(parts[1])
        except FileNotFoundError:
            raise FileNotFoundError(f"File {filename} not found")
        if not masses:
            raise ValueError("Atomic masses file empty")
        return masses

class SynthesisCalculator:
    def __init__(self, atomic_masses: dict):
        self.atomic_masses = atomic_masses
        self.target = {}
        self.target_mm = Decimal(0)
        self.reagents = []
        self.reagent_formulas = []
        self.molar_masses = []

    def input_target(self):
        f = input("Enter target formula: ").strip()
        self.target = FormulaParser.parse(f)
        for e, c in self.target.items():
            if e not in self.atomic_masses:
                raise ValueError(f"Mass of element {e} not defined")
            self.target_mm += c * self.atomic_masses[e]
        print(f"Parsed: {self.target}")
        print(f"Target molar mass = {self.target_mm.normalize()} g/mol")

    def input_reagents(self):
        n = int(input("Enter number of reagents: "))
        for i in range(n):
            f = input(f"Reagent {i+1} formula: ").strip()
            self.reagent_formulas.append(f)
            parsed = FormulaParser.parse(f)
            mm = sum(parsed[e] * self.atomic_masses[e] for e in parsed)
            self.reagents.append(parsed)
            self.molar_masses.append(mm)
            print(f"{f} molar mass = {mm.normalize()} g/mol")

    def validate(self):
        tgt = set(self.target)
        prov = set().union(*(r.keys() for r in self.reagents))
        miss = tgt - prov
        extra = prov - tgt
        if miss:
            print(f"Error: Missing elements: {miss}")
            exit()
        if extra:
            print(f"Warning: Extra elements: {extra}")

    def solve_basis(self):
        # Construct A, b corresponding to 1 mol 
        elems = sorted(set(self.target) | set().union(*(r.keys() for r in self.reagents)))
        A = [[r.get(e, Decimal(0)) for r in self.reagents] for e in elems]
        b = [self.target.get(e, Decimal(0)) for e in elems]
        # Convert to rational number to avoid floating point error
        A_rat = [[Rational(str(val)) for val in row] for row in A]
        b_rat = [Rational(str(val)) for val in b]
        vars = symbols(f'x0:{len(self.reagents)}')
        sol = linsolve((Matrix(A_rat), Matrix(b_rat)), vars)
        if not sol:
            print("Error: No exact solution for 1 mol target.")
            exit()
        tup = next(iter(sol))
        # Convert rational to decimal exact representation
        base = {}
        for i, r in enumerate(tup):
            base[vars[i]] = Decimal(r.p) / Decimal(r.q)
        return base, elems

    def calculate(self):
        t = input("Input type (0 mass g, 1 moles): ").strip()
        if t not in ['0','1']:
            print("Invalid type")
            exit()
        v = Decimal(input("Amount: ").strip())
        n = v / self.target_mm if t=='0' else v
        print(f"Target moles = {n.normalize()}")

        base, elems = self.solve_basis()
        print("\nReagent requirements:\n" + "="*40)
        total_mass = Decimal(0)
        for i, var in enumerate(sorted(base.keys(), key=lambda x: int(str(x)[1:]))):
            moles = base[var] * n
            mass = moles * self.molar_masses[i]
            total_mass += mass
            print(f"Reagent {i+1} ({self.reagent_formulas[i]}):")
            print(f"  Moles: {moles.normalize():.10f} mol")
            print(f"  Mass:  {mass.normalize():.10f} g\n")
        print(f"Total mass = {total_mass.normalize():.10f} g")

        #Verify element balance
        print("Verification:")
        for e in elems:
            actual = sum((base[symbols(f'x{i}')] * n) * r.get(e, Decimal(0))
                         for i, r in enumerate(self.reagents))
            target_amt = self.target.get(e, Decimal(0)) * n
            print(f"{e}: target {target_amt.normalize():.10f}, actual {actual.normalize():.10f}")


def main():
    print("High-Precision Compound Synthesis Calculator")
    print("="*60)
    try:
        masses = AtomicMassReader.read()
        print(f"Loaded elements: {', '.join(masses.keys())}")
    except Exception as e:
        print(e)
        exit()
    calc = SynthesisCalculator(masses)
    calc.input_target()
    calc.input_reagents()
    calc.validate()
    calc.calculate()

if __name__ == "__main__":
    main()

附一张使用截图，后面也会放到Github上面

Comment and share

本文仅为学习记录，并非教程

这次记录文主要侧重于应用，参考了Sob老师的博文和网络上的各大教程

以FCC Cu为例进行ELF计算与分析

采用FCC面心立方最密堆积的铜单胞，使用vasp计算，先进行结构优化，确认收敛后进行ELF计算

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

Global Parameters
ISTART =  1            (Read existing wavefunction, if there)
ISPIN  =  2            (Non-Spin polarised DFT)
MAGMOM = 4*0.6
# ICHARG =  11         (Non-self-consistent: GGA/LDA band structures)
LREAL  = .FALSE.       (Projection operators: automatic)
 ENCUT  =  530        (Cut-off energy for plane wave basis set, in eV)
# PREC   =  Accurate   (Precision level: Normal or Accurate, set Accurate when perform structure lattice relaxation calculation)
LWAVE  = .TRUE.        (Write WAVECAR or not)
LCHARG = .TRUE.        (Write CHGCAR or not)
ADDGRID= .TRUE.        (Increase grid, helps GGA convergence)
LASPH  = .TRUE.        (Give more accurate total energies and band structure calculations)
PREC   = Accurate      (Accurate strictly avoids any aliasing or wrap around errors)
# LVTOT  = .TRUE.      (Write total electrostatic potential into LOCPOT or not)
# LVHAR  = .TRUE.      (Write ionic + Hartree electrostatic potential into LOCPOT or not)
# NELECT =             (No. of electrons: charged cells, be careful)
# LPLANE = .TRUE.      (Real space distribution, supercells)
# NWRITE = 2           (Medium-level output)
# KPAR   = 2           (Divides k-grid into separate groups)
# NGXF    = 300        (FFT grid mesh density for nice charge/potential plots)
# NGYF    = 300        (FFT grid mesh density for nice charge/potential plots)
# NGZF    = 300        (FFT grid mesh density for nice charge/potential plots)

Electronic Relaxation
ISMEAR =  1            (Gaussian smearing, metals:1)
SIGMA  =  0.2         (Smearing value in eV, metals:0.2)
NELM   =  90           (Max electronic SCF steps)
NELMIN =  6            (Min electronic SCF steps)
EDIFF  =  1E-08        (SCF energy convergence, in eV)
# GGA  =  PS           (PBEsol exchange-correlation)

Ionic Relaxation
NSW    =  0         (Max ionic steps)
IBRION =  -1            (Algorithm: 0-MD, 1-Quasi-New, 2-CG)
ISIF   =  2            (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)
EDIFFG = -2E-02        (Ionic convergence, eV/AA)
# ISYM =  2            (Symmetry: 0=none, 2=GGA, 3=hybrids)
LELF = .TRUE.
NPAR = 1

在优化完的结构上进行自洽计算，添加LELF = .TRUE.这个参数，确保ELFCAR文件输出

官方wiki还提到必须显示开启NPAR = 1，然而网上的各大教程都没有提到这一点

If LELF is set, NPAR=1 has to be set explicitely in the INCAR file in addition

且有一篇讨论帖：https://wiki.vasp.at/forum/viewtopic.php?t=20020

不过我自己尝试的时候发现超算平台上添加NPAR = 1会报错，而在自己电脑上则不会报错，有点迷、、、

用Vesta载入ELFCAR文件

先设置一下等值面数值

Isosurface Level（等值面水平）是一个关键参数，用于控制可视化电子局域化程度的阈值，VESTA会绘制出ELF等于该值的等值面（isosurface），直观展示电子局域化程度高于该阈值的区

Properties-Isosurfaces-Isosurfaces level，设置完后可以按以下Tab键就会自动保存应用

这里如果设置0.75的话啥也看不到，因为铜单胞中存在着金属键，并没有很强的局域化电子，而是均匀分布的电子气

高定域性通常出现在共价作用区域、孤对电子区域、原子内核及壳层结构区域

然后就可以查看ELF图了

首先是3D的图

左上角-Edit-Lattice planes，新建一个图，根据想查看的面的来设置米勒指数hlk，以及设置距离d

然后是2D的图

左上角-Utilities-2D Data Display-Slice即可

可以看到，顶点Cu与顶点Cu原子之间、顶点Cu和面心Cu之间的电子的局域化程度比较高并且没有明显的偏向偏离

接下来看一下一维时的电子定域化情况

左上角-Utilities-Line profile即可

• 首先棱上的情况

Cu和Cu两原子中间存在均匀的非零的电子局域函数，虽然不足0.5但分布平坦，这是金属键合的特征

两个铜原子的ELF(r)值相近

• 面对角线上的情况

三个铜原子的ELF(r)值相近

离子晶体CsF

再看一下离子键的情况，以CsF为例

氟和铯之间出现深蓝色，这说明在这个区域电子的分布较为离域化，即氟和铯之间不存在共用电子，键的类型偏向离子键

再看看对角线上的电荷分布

两端为铯中间为氟，可以看到氟和铯ELF(r)值差别较大

共价晶体——硅

如箭头所指，硅原子和硅原子之间有较强的电子定域化作用，键的类型为共价键

如下图，体对角线上的ELF图也能很好证明，Si-Si之间的区域电子化定域程度极高，呈现明显的尖峰

共价晶体——灰锡

我仍记得高中时课本讲的一个稍微有些特殊的例子——灰锡，灰锡是一个共价晶体

可以看到，灰锡与晶体硅结构上类似，Sn($\frac{1}{4},\frac{1}{4},\frac{1}{4}$)与Sn(0,0,0)与以及Sn(0.5,0.5,0)成键，Sn($\frac{3}{4},\frac{3}{4},\frac{1}{4}$)与Sn(0.5,0.5,0)和Sn(1,1,0)成键。成键电子局域化程度较高，为共价电子

Comment and share

2025年5月Intel平台VASP.6.4.3的编译安装

WSL的安装

选择Arch Linux，个人比较熟悉，仅供学习使用不作为稳定的生产力环境，配置如下，两颗E5 2676 v3

打开powershell或者终端，输入wsl --install archlinux以自动安装，这一步需要使用代理(MAGIC)，主机上使用Clash-Verge进行代理（运气好也能连上），终端中输入：

1
2

set http_proxy=http://127.0.0.1:7897 #自己使用的端口是7897，http代理
set https_proxy=http://127.0.0.1:7897 #自己使用的端口是7897，https代理

如果没有条件使用代理，也可以手动安装，参考Arch Wiki，https://wiki.archlinuxcn.org/wiki/%E5%9C%A8_WSL_%E4%B8%8A%E5%AE%89%E8%A3%85_Arch_Linux，手动下载.wsl镜像，然后安装

笔者计算机C盘性能太差，于是想要将wsl2迁移到另一个分区

安装成功

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

[root@NAS ~]# fastfetch
                  -`                     root@NAS
                 .o+`                    --------
                `ooo/                    OS: Arch Linux x86_64
               `+oooo:                   Host: Windows Subsystem for Linux - archlinux (2.4.12)
              `+oooooo:                  Kernel: Linux 5.15.167.4-microsoft-standard-WSL2
              -+oooooo+:                 Uptime: 5 hours, 16 mins
            `/:-:++oooo+:                Packages: 155 (pacman)
           `/++++/+++++++:               Shell: bash 5.2.37
          `/++++++++++++++:              WM: WSLg 1.0.65 (Wayland)
         `/+++ooooooooooooo/`            Terminal: xterm-256color
        ooosssso++osssssso+`           CPU: Intel(R) Xeon(R) E5-2676 v3 (48) @ 2.39 GHz
       .oossssso-````/ossssss+`          Memory: 1.09 GiB / 15.57 GiB (7%)
      -osssssso.      :ssssssso.         Swap: 304.95 MiB / 4.00 GiB (7%)
     :osssssss/        osssso+++.        Disk (/): 18.67 GiB / 1006.85 GiB (2%) - ext4
    /ossssssss/        +ssssooo/-        Disk (/mnt/c): 159.05 GiB / 237.01 GiB (67%) - 9p
  `/ossssso+/:-        -:/+osssso+-      Disk (/mnt/d): 13.16 TiB / 14.55 TiB (90%) - 9p
 `+sso+:-`                 `.-/+oso:     Disk (/mnt/e): 1.70 TiB / 3.64 TiB (47%) - 9p
`++:.                           `-/+/    Disk (/mnt/f): 49.45 GiB / 476.94 GiB (10%) - 9p
.`                                 `/    Disk (/mnt/g): 25.16 TiB / 29.10 TiB (86%) - 9p
                                         Local IP (eth0): 172.31.36.99/20
                                         Locale: en_US.UTF-8

编译器的安装

1
2
3
4
5
6
7
8
9
10
11

pacman -S make
pacman -S cmake
pacman -S gcc
pacman -S gcc-fortran
####发现gcc和gfortran编译器太新，会编报错(C23)，故降级软件包
#添加中科大源来安装yay
pacman -Syy yay
yay -S downgrade #安装降级工具
downgrade -S gcc-libs #选择12.1.0版本，安装gcc的依赖
downgrade -S gcc #安装gcc，选择12.1.0版本
downgrade -S gcc-fortran #安装gforran，选择12.1.0版本

前往Intel官网下载oneAPI BaseTool kithttps://www.intel.com/content/www/us/en/developer/tools/oneapi/base-toolkit-download.html?packages=oneapi-toolkit&oneapi-toolkit-os=linux&oneapi-lin=online

1
2

chmod +x xxxxxxxxxx.sh #赋予下载下来的sh文件可执行权限
xxxx.sh #运行sh文件

oneAPI BaseTool kit不需要全部安装，只需安装MKL部分，

HPC也只需安装一部分

此处博主踩了较多的坑：

• 编译器采用2022版，不可采用过新编译器
  • 2024年10月后，Intel删除了ifort编译器
• gcc编译器不可太新，否则编译某些fortran模块会报错，博主从14降级到了12.1.0

VASP.6.4.3的编译安装

源码包可以去各大私密性较强的论坛或者微信公众号获取，仅供个人学习使用，不用于学术、商业用途

选择arch/目录中的模板makefile.include.intel，复制到.(上一层目录)，并修改名称为makefile.include，修改目录如下

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

# Default precompiler options
CPP_OPTIONS = -DHOST=\"LinuxIFC\" \
              -DMPI -DMPI_BLOCK=8000 -Duse_collective \
              -DscaLAPACK \
              -DCACHE_SIZE=4000 \
              -Davoidalloc \
              -Dvasp6 \
              -Duse_bse_te \
              -Dtbdyn \
              -Dfock_dblbuf

CPP         = fpp -f_com=no -free -w0  $*$(FUFFIX) $*$(SUFFIX) $(CPP_OPTIONS)

FC          = mpiifort
FCL         = mpiifort

FREE        = -free -names lowercase

FFLAGS      = -assume byterecl -w

OFLAG       = -O2 -xhost  #添加-xhost标记
OFLAG_IN    = $(OFLAG)
DEBUG       = -O0

OBJECTS     = fftmpiw.o fftmpi_map.o fftw3d.o fft3dlib.o
OBJECTS_O1 += fftw3d.o fftmpi.o fftmpiw.o
OBJECTS_O2 += fft3dlib.o

# For what used to be vasp.5.lib
CPP_LIB     = $(CPP)
FC_LIB      = $(FC)
CC_LIB      = icc
CFLAGS_LIB  = -O
FFLAGS_LIB  = -O1
FREE_LIB    = $(FREE)

OBJECTS_LIB = linpack_double.o

# For the parser library
CXX_PARS    = icpc
LLIBS       = -lstdc++

##
## Customize as of this point! Of course you may change the preceding
## part of this file as well if you like, but it should rarely be
## necessary ...
##

# When compiling on the target machine itself, change this to the
# relevant target when cross-compiling for another architecture
VASP_TARGET_CPU ?= -xHOST
FFLAGS     += $(VASP_TARGET_CPU)

# Intel MKL (FFTW, BLAS, LAPACK, and scaLAPACK)
# (Note: for Intel Parallel Studio's MKL use -mkl instead of -qmkl)
FCL        += -qmkl=sequential
MKLROOT    ?= /opt/intel/oneapi/mkl/2022.1.0 #一定要添加MKLROOT的目录，否则会报错


/home/storm/intel/oneapi/mkl/2022.1.0
LLIBS      += -L$(MKLROOT)/lib/intel64 -lmkl_scalapack_lp64 -lmkl_blacs_intelmpi_lp64
INCS        =-I$(MKLROOT)/include/fftw

# HDF5-support (optional but strongly recommended)
#CPP_OPTIONS+= -DVASP_HDF5
#HDF5_ROOT  ?= /path/to/your/hdf5/installation
#LLIBS      += -L$(HDF5_ROOT)/lib -lhdf5_fortran
#INCS       += -I$(HDF5_ROOT)/include

# For the VASP-2-Wannier90 interface (optional)
#CPP_OPTIONS    += -DVASP2WANNIER90
#WANNIER90_ROOT ?= /path/to/your/wannier90/installation
#LLIBS          += -L$(WANNIER90_ROOT)/lib -lwannier

然后开始编译

发现使用命令make all -j48会编译报错

于是改为官方wiki的命令make DEPS=1 -j48，编译成功

测试

make test进行测试，有一些报错，懒得管了、、、能跑就行

Comment and share

安装、配置好vaspkit之后输入vaspkit可启动软件包

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

(base) storm@DESKTOP-HE4FQ8Q:~/my-learn/trash-vaspkit$ vaspkit
            \\\///
           / _  _ \         Hey, you must know what you are doing.
         (| (o)(o) |)       Otherwise you might get wrong results.
 o-----.OOOo--()--oOOO.------------------------------------------o
 |         VASPKIT Standard Edition 1.5.1 (27 Jan. 2024)         |
 |         Lead Developer: Vei WANG (wangvei@icloud.com)         |
 |      Main Contributors: Gang TANG, Nan XU & Jin-Cheng LIU     |
 |  Online Tutorials Available on Website: https://vaspkit.com   |
 o-----.oooO-----------------------------------------------------o
        (   )   Oooo.                          VASPKIT Made Simple
         \ (    (   )
          \_)    ) /
                (_/
 ===================== Structural Utilities ======================
 01) VASP Input-Files Generator    02) Mechanical Properties
 03) K-Path for Band-Structure     04) Structure Editor
 05) Catalysis-ElectroChem Kit     06) Symmetry Analysis
 07) Materials Databases           08) Advanced Structure Models
 ===================== Electronic Utilities ======================
 11) Density-of-States             21) Band-Structure
 23) 3D Band-Structure             25) Hybrid-DFT Band-Structure
 26) Fermi-Surface                 28) Band-Structure Unfolding
 31) Charge-Density Analysis       42) Potential Analysis
 44) Piezoelectric Properties      51) Wave-Function Analysis
 62) Magnetic Analysis             65) Spin-Texture
 68) Transport Properties
 ======================== Misc Utilities =========================
 71) Optical Properties            72) Molecular-Dynamics Kit
 74) User Interface                78) VASP2other Interface
 84) ABACUS Interface              91) Semiconductor Kit
 92) 2D-Material Kit               95) Phonon Analysis
 0)  Quit
 ------------>>

以下是中文翻译对照：

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

===================== 结构工具 ======================
01) VASP 输入文件生成器
02) 力学性质
03) 能带结构的 K 路径
04) 结构编辑器
05) 催化 - 电化学套件
06) 对称性分析
07) 材料数据库
08) 高级结构模型
==================== 电子工具 ======================
11) 密度态
21) 能带结构
23) 三维能带结构
25) 混合 DFT 能带结构
26) 费米面
28) 能带结构展开
31) 电荷密度分析
42) 势分析
44) 压电性质
51) 波函数分析
62) 磁性分析
65) 自旋纹理
68) 输运性质
====================== 杂项工具 ========================
71) 光学性质
72) 分子动力学套件
74) 用户界面
78) VASP 转其他格式接口
84) ABACUS 接口
91) 半导体套件
92) 二维材料套件
95) 声子分析
0) 退出

我目前还是初学者小白，使用得最多的功能还是生成输入文件功能尤其是生成INCAR文件，

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

 ------------>>
01
 ==================== VASP Input Files Options ===================
 101) Customize INCAR File
 102) Generate KPOINTS File for SCF Calculation
 103) Generate POTCAR File with Default Setting
 104) Generate POTCAR File with User Specified Potential
 105) Generate POSCAR File from cif (no fractional occupations)
 106) Generate POSCAR File from Material Studio xsd (retain fixes)
 107) Reformat POSCAR File in Specified Order of Elements
 108) Successive Procedure to Generate VASP Files and Check
 109) Submit Job Queue

 0)   Quit
 9)   Back
 ------------>>
  ------------>>
101
 +---------------------------- Tip ------------------------------+
 |          WARNNING: You MUST know what wou are doing!          |
 |Some Parameters in INCAR file need to be set/adjusted manually.|
 +---------------------------------------------------------------+
 ======================== INCAR Options ==========================
 ST) Static-Calculation            SR) Standard Relaxation
 MG) Magnetic Properties           SO) Spin-Orbit Coupling
 D3) DFT-D3 no-damping Correction  H6) HSE06 Calculation
 PU) DFT+U Calculation             MD) Molecular Dynamics
 GW) GW0 Calculation               BS) BSE Calculation
 DC) Elastic Constant              EL) ELF Calculation
 BD) Bader Charge Analysis         OP) Optical Properties
 EC) Static Dielectric Constant    PC) Decomposed Charge Density
 PH) Phonon-Calculation            PY) Phonon with Phononpy
 NE) Nudged Elastic Band (NEB)     DM) The Dimer Method
 FQ) Frequence Calculation         LR) Lattice Relaxation
 MT) Meta-GGA Calculation          PZ) Piezoelectric Calculation

 0)   Quit
 9)   Back
 ------------>>
 Input Key-Parameters (STH6D3 means HSE06-D3 Static-Calcualtion)

以下是中文对照：

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

+---------------------------- 提示 ------------------------------+
|          警告：你必须清楚自己在做什么！          |
|INCAR 文件中的一些参数需要手动设置/调整。|
+---------------------------------------------------------------+
======================= INCAR 选项 ==========================
ST) 静态计算            SR) 标准松弛
MG) 磁性性质            SO) 自旋轨道耦合
D3) DFT-D3 无阻尼校正  H6) HSE06 计算
PU) DFT+U 计算          MD) 分子动力学
GW) GW0 计算            BS) BSE 计算
DC) 弹性常数            EL) ELF 计算
BD) Bader 电荷分析      OP) 光学性质
EC) 静态介电常数        PC) 分解电荷密度
PH) 声子计算            PY) 使用 Phonopy 的声子计算
NE) 强制弹性带（NEB）   DM) 二聚体方法
FQ) 频率计算            LR) 晶格松弛
MT) Meta-GGA 计算       PZ) 压电计算

LR晶格弛豫就对应着结构优化，当然生成的INCAR文件可能仍需要进行一些修改，比如ENCUT可能需要调整为赝势中ENMAX的1.3倍

生成KPOINTS文件

对于非能带的计算，选择102) Generate KPOINTS File for SCF Calculation使用程序自动撒点即可，但是需要用户选择撒点方式和撒点密度，选择Gamma Scheme，使用0.04精度

关于生成 KPOINTS 文件

• 在 VASP 计算中，KPOINTS 文件用于指定 k 点网格的设置。对于非能带计算，比如进行能量计算、优化结构等，通常不需要像计算能带结构那样沿特定的高对称路径采样 k 点，而是可以在整个布里渊区均匀地分布 k 点，这就是所说的 “自动撒点”。

撒点方式

• Monkhorst-Pack 方式 ：这是最常用的自动撒点方式。它通过指定一组三个整数（nx, ny, nz），在三个倒易晶格方向上均匀地分布 k 点，形成的网格可以系统地逼近布里渊区的积分。例如，对于简单晶格，若指定 k 点网格为 3×3×3，程序就会按照 Monkhorst-Pack 方式在这三个方向各取 3 个点，生成一个均匀分布的 k 点网格。
• Gamma 方式 ：这种方式会在 k 点网格中包含 Γ 点（布里渊区的中心点）。对于金属等具有高对称性的材料，在 Γ 点附近采样很重要，这时候 Gamma 方式撒点更合适。比如指定 k 点网格为 4×4×4 并采用 Gamma 方式撒点，生成的网格会以 Γ 点为中心，向周围均匀分布其他 k 点。

K 点密度

• K 点密度指的是在布里渊区单位体积内分布的 k 点数量，它决定了 k 点网格的疏密程度。较高的 K 点密度意味着更精细地采样布里渊区，能够更精确地计算系统的物理性质，但也会增加计算量和资源消耗。
• 例如，对于一个计算资源充足的大型超算任务，为了获得高精度结果，可以选择较大的 K 点密度，如 8×8×8 甚至更高；而对于初步探索或计算资源有限的情况，可以适当降低 K 点密度，如 4×4×4。

生成KPOINTS的同时，POTCAR也会被自动生成，前提是在vaspkit的配置文件中正确设置了赝势的路径和赝势稳健的种类（默认type为PBE）

Comment and share

使用VASP计算离子电导率，似乎比较麻烦

一个能量较低、结构合理的起始三维构型是进行MD模拟的第一步

分子动力学模拟的基本步骤

• 确定起始构型
• 选用适当的力场和模拟软件
• 构建体系和能量最小化
• 平衡过程
• 数据采集过程
• 结果分析

根据系统的传播方式，分子模拟方法可以分为分子动力学MD和蒙特卡洛MC两类方法

• 分子动力学通过数值积分求解牛顿运动方程，生成系统轨迹，可研究结构、动力学和热力学性质。
• 蒙特卡洛方法基于概率规则生成新构型，侧重采样平衡态结构和热力学性质，不涉及时间动力学。

力场是一组描述分子中原子间相互作用的势能函数和参数的集合，用于计算系统的总能量和原子受力，从而驱动分子动力学（MD）模拟或蒙特卡罗（MC）采样：

• 将量子力学描述的分子相互作用简化为经典力学模型，降低计算成本。
• 通过参数化的势能函数，预测分子的结构、动力学和热力学性质（如键长、自由能、扩散系数等）。

固定电荷力场：固定原子电荷，不考虑环境引起的电子云变形（极化），计算效率高

• 生物分子力场
  • AMBER
    • 广泛用于蛋白质、核酸模拟，参数优化针对溶液环境，如 ff14SB、ff19SB
  • CHARMM
    • 广泛用于蛋白质、核酸模拟，参数优化针对溶液环境，如 ff14SB、ff19SB
  • GROMOS
    • 用于膜蛋白和聚合物，参数简洁，计算速度快
• 小分子力场
  • GAFF
    • 可参数化多种有机小分子，支持自定义分子
  • UFF
    • 覆盖周期表多数元素，适合无机材料和有机分子

极化力场：显示考虑原子极化率（电子云随电场变化），通过附加自由度（如 Drude 振子、诱导偶极）模拟极化效应，精度更高但计算成本显著增加

• AMOEBA：基于原子极化率和多极矩，适用于极性体系（如水、蛋白质 - 配体相互作用）。
• Drude 模型：每个原子附加虚拟振子模拟极化，如 CHARMM-Drude 力场。
• Polarizable AMBER (PAMBS)：在 AMBER 框架中引入极化参数。

粗粒化力场（Coarse-Grained Force Fields）：将多个原子简化为单个 “珠子”（如蛋白质的每个残基、脂质的疏水尾），降低计算复杂度，适用于大尺度系统（如细胞膜、聚合物）

• Martini：广泛用于膜蛋白、病毒颗粒和软物质，4:1 粗粒化（4 个原子→1 个珠子）。
• MARTINI 3.0：支持多尺度模拟，兼容全原子模型。

使用pymatgen提取数据

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

'''
分析AIMD结果，计算MSD 和 conductivity
'''
import os
from pymatgen.core.trajectory import Trajectory
from pymatgen.io.vasp.outputs import Xdatcar
from pymatgen.core import Structure
from pymatgen.analysis.diffusion.analyzer import DiffusionAnalyzer
import numpy as np
import pickle
import matplotlib.pyplot as plt

# 这一步是读取 XDATCAR，得到一系列结构信息
traj = Trajectory.from_file('XDATCAR')

# 这一步是实例化 DiffusionAnalyzer 的类
# 并用 from_structures 方法初始化这个类； 900 是温度，2 是POTIM 的值，1是间隔步数
# 间隔步数（step_skip）不太容易理解，但是根据官方教程:
# dt = timesteps * self.time_step * self.step_skip

diff = DiffusionAnalyzer.from_structures(traj,'Li',900,2,1)

# 可以用内置的 plot_msd 方法画出 MSD 图像
# 有些终端不能显示图像，这时候可以调用 export_msdt() 方法，得到数据后再自己作图
ax = diff.export_msdt("msd.dat")
#plt.show()

# 接下来直接得到 离子迁移率， 单位是 mS/cm
C = diff.conductivity

with open('result.dat','w') as f:
    f.write('# AIMD result for Li-ion\n')
    f.write('temp\tconductivity\n')
    f.write('%d\t%.2f\n' %(900,C))

画图

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

import numpy as np
import matplotlib.pyplot as plt

# 读取数据
data = np.loadtxt('msd.txt', comments='#')  # 假设你的文件名为 msd.txt

# 提取时间步长和 MSD 数据
time = data[:, 0]  # 第一列是时间步长
msd = data[:, 1]   # 第二列是总的 MSD

# 绘制 MSD 曲线
plt.figure(figsize=(10, 6))
plt.plot(time, msd, label='MSD', color='blue', linewidth=2)
plt.xlabel('Time (fs)', fontsize=14)
plt.ylabel('Mean Square Displacement (Ų)', fontsize=14)
plt.title('MSD vs Time', fontsize=16)
plt.legend(fontsize=12)
plt.grid(True)
plt.show()

1. System preparation
2. Minimization/Relaxation
3. Equilibration
4. Production

系统准备可以说是模拟中最关键的阶段，但在许多情况下却得到了最少的关注；如果你的系统准备有缺陷，这些缺陷可能会证明是致命的。最坏的可能结果是，准备好的系统并非你想要的（例如，它包含不正确的分子或质子化状态），但在化学上却是有效的，并且可以被你的力场很好地描述，因此没有错误地通过了剩余的步骤——实际上，这是系统准备问题中常见的后果。如果系统在随后的平衡步骤中行为良好，不应假设系统已正确准备；这里应该格外小心。

生产流程如下：

为什么需要有恒温器？

分子动力学模拟用于观察和获取某个研究系统的有趣性质。在很多情况下，为了模拟实验室条件下进行的实验，希望从正则（恒温）配分函数中进行抽样。一般来说，如果系统在模拟过程中必须保持温度不变，将采用某种温控算法。

分子动力学模拟的温度通常使用根据等分配定理定义的动能来测量： $$ \frac{3}{2}Nk_BT = <\sum_{i=1}^{N}\frac{1}{2}m_iv_i^2> $$ 带角标的括号表示温度是一个时间平均值。如果我们使用等分配定理来计算分子动力学模拟单时间点的温度快照，而不是时间平均，这个量称为瞬时温度。瞬时温度不一定等于目标温度；实际上，在正则配分函数中，瞬时温度应该围绕目标温度波动。

温控算法通过改变本质上微观可逆的牛顿运动方程来工作。因此，如果希望计算扩散系数等动力学性质，则最好不使用温控器；相反，在系统达到所需温度后，应关闭温控器。然而，虽然所有温控器都会给出非物理动力学，但有些研究发现它们对特定动力学性质的计算影响很小，并且它们在生产模拟中也常被使用。

一些恒温器

• Gaussian
  • 使用了高斯最小约束原理来确定维持瞬时温度所需的最小扰动力。这种温控器不采样正则分布；相反，它采样的是等速（恒动能）配分函数。但是，等速配分函数采样与正则配分函数相同的配置相空间，因此，使用任一配分函数都可以等效地获得位置依赖（结构）平衡性质。然而，速度依赖（动力学）性质在两个配分函数之间不会等效。这种温控器仅适用于某些高级应用。
• Simple Velocity Rescaling
  • 实施起来最容易的温控器之一,也是最不物理的温控器之一。这种温控器依赖于重新标粒子的动量，使模拟的瞬时温度恰好与目标温度相匹配，不建议使用
• Berendsen
  • 既不采样正则分布也不采样等速分布，它的使用可能导致模拟伪象，不建议使用
• Bussi-Donadio-Parrinello (Canonical Sampling Velocity Rescaling)
  • 正确采样了正则配分函数，也被称为“V-rescale”热浴方法，非常常用的热浴方法
• Andersen
  • 只能用于采样结构性质，因为动力学性质可能会受到突然碰撞的很大影响，不建议使用
• Langevin
  • 不了解
• Noés-Hoover
  • 不了解

Comment and share

ENCUT是用来指定平面波基组能量截止值的参数

以上是wiki对ENCUT的描述，POTCAR文件中已经包含了元素的ENMAX和ENMIN，这也是互联网上对ENCUT的选择方案对通常是1.3倍的ENMAX

侯柱峰老师提供了另外一种方法并指出通过该方法选择的ENCUT通常能满足1.3倍ENMAX，供参考使用

截断能指定了用于波函数展开的平面波基组的截断能量，此能量越大则用来描述波函数的平面波基组越多，精度越高，但计算也越耗时。

可以通过计算测试来选择合适的ENCUT。

以下是用于测试的脚本：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

#!/bin/sh
rm WAVECAR
for i in 150 200 250 300 350 400
do
cat > INCAR <<!
SYSTEM = Si-Diamond
ENCUT = $i
ISTART = 0 ; ICHARG = 2
ISMEAR = -5
PREC = Accurate
!
echo "ENCUT = $i eV" ; time vasp
E=`grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.6f \n", $5 }'`
echo $i $E >>comment
done

计算后会得到comment文件，它列出了在不同ENCUT下计算得到的总能量

1
2
3
4
5
6

150 -11.900655
200 -11.938864
250 -11.944599
300 -11.945248
350 -11.945503
400 -11.945622

==总能量变化在0.001 eV即可==，因此在这个例子中ENCUT可以选择250

Comment and share