计算化学合篇（一）

免费项目 2024-08-11 0

点击此处查看最新的网赚项目教程

————–特别注释：以下内容仅供闲聊———

一、初学者如何学习计算化学

就大三大四/研究生学习计算化学来说，阅读很多关于模拟的复杂数学书籍非常困难。而就直接去操作软件，像LAMMPS或CP2K等软件将仿佛又是个黑匣子。越接触越困惑。有时为了找一些支撑去检索《PR》或《JCP》的文献，又很难理解。

第一点，你要区分学习方法和学习如何在特定软件中实现该方法？要学的东西有很多，重要的有统计力学、热力学和分子模拟的知识，方式是可以听视频讲座和看相关书籍。自学学习这些东西，确实会有很多痛苦。有些模拟用到的方法，一旦你理解了原理，在操作时看到处理时就显得很轻松，但如果你不理解原理，就会感到困惑。如果先学工具，再学方法，那就是本末倒置，结果很难令人满意。

再者，鼓励新手不要急于计算成功你的job，而要花费一些时间从教科书中学习统计力学和分子模拟的原理。

第二点，MD/kMC等很像实验，实验也要准备药品，对于那些基本的实践技能和简单案例的练习是一个很好的初始准备。

第三点，必须接受的一件事是，不会有一本教科书或文章能够准确且仅包含按你要求的最好或最轻松方式的学习内容。你不得不去搜寻一些（或通过公开视频或从书本里），并且还必须学会如何寻找你需要学习的东西（检索到你需要的很重要）。后者与理解模拟技术一样是一项有用的技能。

第四点，你想更确切地了解，除了阅读源代码之外别无选择。但很有可能读不懂代码和相关的理论。需要花费很长时间成本。

二、关于系综NPT的一些问题

NPT：模拟中体系边界各向异性变化，能量不守恒。NVT：BOX体积不变，能量不守恒。相较于NVT而言，NPT系综似乎是更符合实际，而往往在使用时面临许多情况是难以解决。大多数运算错误是由于不良的初始结构、不良/不正确的力场参数、错误的拓扑、不良的模拟设置或它们之间组合一起引起的。如前面讲的原子丢失，一般是由于原子重叠造成的，这会导致高温，最终使原子被吹出盒子。而像H原子的轻质量元素最容易发生。初始结构会随时间的推移影响模拟结果。

1.热浴选择

首先，常见的热浴技术，Anderson是按步数随机选择盒子中的例子将粒子速度指派为按Maxwell分布的随机选取值。Anderson方法造成了轨迹不平滑，惯性所致。这里应当指出速度呈现Maxwell分布，是严格描述了理想气体条件下粒子的速度分布，对于具有粒子间相互作用且达到了热平衡的情况下也近似适用的。Berendsen不是将温度立即达到目标温度，而使热浴耦合效过有一定速率。虽然比原始速率调节方法好，但造成粒子的速度偏离Maxwell分布，强行乘以系数会使粒子速度之间差异增大。Nose-Hoover热浴是将热浴作为BOX体系的一部分，随体系一同演化，适合目标温度平衡体系，并旨在模拟较大体系所需的波动，以便可以对NVT统计力学系综进行采样。

2.控压失衡

NPT时，产生各种弹性波，过程中要比较下压力波动和温度波动的曲线。压力波动的幅度很大程度上取决于可压缩性、体系尺寸和时间步长的选择。凝聚态体系的可压缩性不是很好，固体体系的可压缩性更差。对于凝聚态甚至是固态体系压力的波动比温度的波动大得多。对于有限体系，尤其是较小的体系，压力和温度的波动也是有预期的。就特殊情况而言，如熔化淬火，压力值也受噪音太大、淬火太快、体系太小、过冷造成的超调等因素影响。对于具有低压缩性的固体，1 atm 和0 atm 压力之间的差异可以忽略不计。瞬时压力会产生显着波动，因此，保持统计上足够收敛的平均压力将需要大量的计算工作，但基本上没有任何增益。且使用经典势的体系误差要大得多。应当注意的是体系中的压力和温度仅在无限数量的粒子的限制下才是恒定的。在任何统计力学系综中，所看到的是系综平均值和波动（在粒子和时间上的），因此需要一个处于并保持平衡的体系。恒压器是会改变整个盒子的体积，从而影响盒子中的所有粒子，若想局部区域单独受压不容易实现。（1）体系问题：如果初始结构、力场、初始条件等合理模拟中，应对压力巨大波动，可采用建一个较大的盒子（体系越小，波动越大），或尝试使用更大的Pdamp值。

使用Langevin恒温器有助于抑制和随机化这些振动，从而减少压力波动。NPT是提供准确的平均值，若小体系，可能看到-1000 atm，然后是+1000 atm，平均值为0。（2）时间步长问题：时间步长影响压力的波动，原因是原子在MD单个timestep中移动得太远。通常可以通过减少时间步长来避免这种情况。对于非常轻的原子来说，有时1fs的时间步长可能太大，这里需要改变步长做验证。（3）力场问题：如果你的力场参数不合理导致非常大的力，比如力场中缺少排斥力，那么很容易遇到时间积分问题。原子越轻，力越大，需要的时间步长越短，否则，模拟（时间积分）将发散。即温度将以不受控制的方式增加。（4）阻尼系数问题：Lammps manual给出了Tdamp和Pdamp一般设置规则，100*timestep和1000*timestep，但是经验规则，具体还得测试。Nose-Hoover恒温器对于Tdamp任意值都表现欠佳，如果Tdamp太小，温度会剧烈波动；如果太大，温度将需要很长时间才能达到平衡。大多数选100。（5）截断半径：这里不单独讲截断的问题，截断另外说明，因为cutoff在初始设置、后处理都有很多应用，重点是根据模拟体系进行设置，需做有根据的测试。另外，NPT运算中可以输出势能的各自贡献，看哪个是导致体系行为异常的。

3.预平衡

NPT前使用NVT进行预平衡，可以减少体系的NPT的压力波动，引导体系达到平衡，直接NPT到达平衡往往会花费更长的时间。开始NPT-MD模拟时，有人也采用最小化模拟，有时可能消除一些极端力的点，可以做这样测试，但有时最小化可能解决不了问题，这时需要事先的NVT。注：在有限温度下，粒子在大多数时间内都不会停留在局部最小值，因此总势能会更高。若初始几何形状是亚稳定的，只有允许盒子尺寸发生变化时，体系才能进入较低的势能状态。当相反电荷的原子靠得太近时，使用不合适的力场参数可能会引起模拟崩溃，因此受到过大的力，然后移动得太快。

4. 0 atm和1 atm

设置0 atm 压力不意味着在真空下进行模拟，对于npt，目标压力的影响是：如果瞬时压力较小，积分器将缩小体系；如果瞬时压力较大，则积分器将扩大体系。因此，0 atm目标压力将平衡到没有外部压力的状态。如前所述，在这种情况下，1个大气压和0个大气压之间的差异可以忽略不计。

5. 初始速度

模拟目标温度时，使目标温度初始速度按Maxwell分布，初始压力还取决于初始位置和单元格大小，模拟时可能产生负压的情况，负压代表体系想要收缩。

6.PBC

PBC周期性边界条件是当粒子跨过原始盒子边界，粒子会从原始盒子反方向边界回到盒子，也就保持盒子内粒子守恒。PBC对BOX形状几乎没要求，但生物分子、聚合物等高分子量注意一般是它们最长直径d，BOX的尺寸大小至少是1.5倍的d，模拟中需应对构象转动等大幅度变化。盒子最简单为立/长方体，然而有些情况不能随意转换晶格。NPT模拟中很容易通过周期性边界导致原子重叠，是时间积分不稳定，进而出现NaN错误。这要进行最小化，如果需要提高键、角度等力常数（最小化在解决初始原子重叠很有帮助）。若你的结构很合理，也可以跳过最小化。如结构处于平衡状态则不需要恒温器。

三、关于计算化学的几点解释

1.物理和化学的匹配

（1）计算化学是研究物理的还是研究化学的？

计算化学是用计算机技术将物理原理汇编成程序然后研究化学问题。但问到计算化学是研究物理的还是研究化学时，存在疑问“仿佛他可以算化学，也可以算物理”。然而，计算化学只能研究化学问题，能不能计算物理很难说。计算化学能够计算原子分子体积。对于熔点和沸点，化学是研究最外层到倒数第三层电子运动波函数的，最外层电子到倒数第三层电子能不能影响熔点和沸点？显然不能。所计算出焓、吉布斯自由能、熵，能不能用，必须根据实际情况。需看定义，如果能满足定义可以，不满足不能用。

（2）计算化学的条件

代数和几何中几乎所有公式是有条件的，用之前先看条件。而编好代码后，很容易忘记条件。理论问题容易变成操作问题。

（3）量子化学近似

量子力学中的多体问题确实很麻烦，因此有必要使用不同的近似。薛定谔方程是建立在多项假设基础上，而其假设又是在稳定状态的基础上的，中心力场条件下的近似，并没有非稳态的薛定谔方程求解。

（4）分子尺度和晶体/不定性尺度

解决化学问题要考虑尺度问题，如何划分分子尺度、晶体尺度具有挑战性？要考虑哪些问题需要用分子计算，哪些用晶体/不定性计算。

（5）量子化学是模糊理论，时间对空间相对论的模糊理论，区别于模糊数学。空间是能够用N阶行列式或矩阵表示出原数，模糊空间就是指四维及其以上。最简单的是三维坐标与一维的时间。时间和空间测不准，所以模糊。

2. 盒子的尺寸和边界条件

（1）带/不带周期性的盒子的尺寸原则上不影响原子结构的最小化、压力、或原子结构的平衡。如果并行运算，会将原子划分到不同的域，这会改变力的求和顺序，进而导致数值的些许差异。

（2）分子动力学中，使用界面层跑NVT，为隔绝上层和下层的作用，真空层应该建的多大，如，水分子层，需要水分子的直径*水分子层的个数=真空层的垂直高度。若盒子里放一个纳米颗粒，外面的真空区域并不是确保盒子左右原子没有相互作用，而是空间可容纳将所有粒子排成水平一列，粒子半径乘以粒子个数。

（3）双层结构在计算相互作用能时，需考虑尺寸效应，a和b值阶梯取值，如12埃*12埃、18埃*18埃、24埃*24埃等，看相互作用能在哪个尺寸时趋于平衡，若在18-32，可以近似使用18*18的计算，无需再更精确的计算。

（4）周期性盒子尺寸应该多大才合理，一方面，较小的周期性盒子单一模拟出的结果可能是伪象，不具备说服力，改变构型（一般采用蒙特卡洛方法）做3~5次取平均。另一方面，足够大的盒子模拟太耗时（由于尺寸大相当于变相取平均），但计算速度大打折扣，因为模拟时间和粒子数并不是线性关系，成指数增长。因此，建议建小模型多次取平均，这样花费的时间成本更少。

（5）为什么NVT有时评估压力很不准？

从宏观角度来看，对于NVT系综来说，密度和温度是确定的。假设系统处于热平衡，则可确定压力。然而，往往发现平衡时的压力变化很大。从应力计算公式（），认为“virial term”导致内应力发生了很大的变化。压力是一种统计（统计力学），需要足够的统计样本。分子动力学模拟并不容易实现。一般压力求解是采用临近盒子的分子对BOX壁的碰撞，是统计力学的计算，越小的体系压力越算不准。

3.力场开发

当你创建一个力场时，你必须先确定一个函数形式，然后将模型参数拟合到一些数据集。数据通常是一堆凭经验确定的结构，或反应能量，或其他什么，但很多时候这是有问题的，因为没有大量关于非平衡结构等的数据（像深度学习对数据数量要求太高，可以首先使用回归算法）。如果你想做动力学，需对非平衡态进行采样。如果训练集中没有这些参数，你怎么知道你的模型参数有多准确？另一种方法是使用量子力学计算一堆不同配置（包括非平衡态）的能量。可以通过更广泛的构象获取大量数据。但这样会使力场变得更好吗？这取决于QM计算的好坏，更重要的是你想实现什么？

注：增加数据集会使力场朝着数据预期结果方向运算，这可能是你所期望的，但这些数据集就是大量的约束条件，约束条件的同时很有可能造成其自身的化学本质/机制的部分丧失。

四、关于“第一性原理”和“从头算”的几点解释

1. ReaxFF是经验性的，其来自第一性原理或实验数据的经验拟合。事实上一些人也不将赝势 DFT 称为“从头算”。然而，对于分子力场MM（包括ReaxFF）中的库伦势（1/r）确实存在一些物理学原理，该势能从电磁拉格朗日函数中获得，或者可以从标准模型拉格朗日函数中导出模型核相互作用的Yukawa势，通常一旦获得拉格朗日密度，就可推导出所有的相互作用力。库伦势具有物理学原理，且其发展比薛定谔方程和波函数理论早得多，但它不足以使MM成为从头算或第一性原理。

2.ReaxFF的其余一些项如过配位项、欠配位项等属于自由参数，被认为是经验的，ReaxFF的库伦项也不意味ReaxFF可以归类于半经验方程，因为“半经验”意味着它会非常具体且具有物理意义。

3.ReaxFF势的经验势中的物理属于，如：范德华、库仑、成指数衰减的吸引/排斥，甚至键级。这些物理学概念都源自波函数理论或受波函数理论的启发。问题是，通过拟合还是训练这些参数，再现 QM的结果使它们成为经验性的力场。

4.“第一性原理”：电子、质子和中子被认为是构成所有物质的基本粒子（此处不包括无质量光子），薛定谔方程是通过将这些粒子的基本属性（质量、电荷、自旋）与其量子性质（波特性）相结合来计算/预测物质属性的第一次成功尝试（形式上的）。理论上（除了建立波函数的哈密顿量的近似形式），物质的量子力学模型中不存在自由参数。

5. 经验势通常倾向于避免明确描述电子及其波动特征。结果是一堆临时（更好或更差）的表达式/函数，它们在理想情况下合并参数以模拟电子在某些条件下的一般行为，从而模拟它们的经验特征。ReaxFF、REBO、BKS、OPLS 等涵盖了实现上述目标的不同尝试。以 ReaxFF 为例，每次想要研究一个新体系时，几乎都需要一个新的参数化。换言之，无需考虑 FF 电势文件中列出的参数组，怎能称为第一原理方法呢？

6. 所有经验模型的出现也遵循Morse’s law（半导体工艺上常用的，也适用于其它领域），ReaxFF等经验模型是实现较低计算成本（体系规模和时间等）应对复杂体系，而非等待计算机变得更快。另外，CPU向GPU 过渡使采用完整的第一原理方法成为潜在可能。

五、分子动力学时间步长（timestep）的几点解释

#####时间步长需要事先估计和测试，过大的时间步长造成NVE能量不准确，过小则计算耗时太长##########现在由于计算能力的提高，大多数不采用测试“更快”的运算，而选择保守的已刊论文中惯用的数值（马太效应，不用反而质疑）########

1.时间步长如何影响结果？在一般的数值方法中，当我们选择足够小的时间步长时，如果继续减小时间步长，结果将是相同的。在lammps manual中给出了不同单位的默认值，实际上此默认时间步长不够小，需做测试。

2.时间步长并不是经验性的，而是取决于体系中移动最快的粒子。已知粒子的质量及其动能（即温度），可以估计合适的时间步长。粒子质量加倍相当于将时间步长除以 sqrt(2)。

3.MD模拟由于采用的积分器算法而产生的误差，虽然不同的数字将很快导致指数级发散的轨迹，但认为只要它们运行稳定并且离散化产生的误差，产生的结果是可以接受的，它们仍然对相同的统计力学系综进行采样。时间步长的选择通常是通过收集更多的模拟时间来改进统计采样。

4.事实上，缩短时间步长只会提高时间积分到某一点的精度（短期），截断的累积误差变得大于离散化的误差，因此进一步缩短时间步长不会再减少误差，但因此对计算机资源产生巨大浪费。

5.如果积分器是Verlet方法，则认为可能存在总能量存在不守恒得情况。ReaxFF测试时间步长过程是通过改变NVE系综的时间步长来验证能量守恒定律是否已经可靠地建立。时间步长可以取0.01 fs-0.1 fs上取递进数值进行验证，确定哪个时间步长下体系近似能量守恒，如0.1fs-0.3fs的体系的能量近似守恒，则选0.3 fs作为模拟的时间步长。

6.另，timestep与及积分器的选择有关，积分器采用的方法如Runge-Kutta method（龙格-库塔方法）、the sixth-order symplectic integrator。ReaxFF一般采用Verlet method的second order symplecticintegrator，不如sixth-order 准确且总能量数值误差累积小。

六、ReaxFF收敛错误信息-200步不收敛或丢失原子

WARNING: Fix qeq/reax CGconvergence failed after 200 iterations at 1 step and ERROR on proc 0:Non-numeric atom coords – simulation unstable

大多数是两个原因：

1. 体系中的几何形状不好，原子之间可能存在紧密堆积或重叠

2. 力场参数和相关设置有问题

即：这表明你的结构不好，或者力场不适合你的结构。

注：1.ReaxFF 是具有隐式键的力场，注意原子之间的初始距离设置

2. Ovito生成文件在处理非正交晶格时存在错误，建议采用VWD中topotools工具，大多数人采用Material Studio建模然后去VWD转换格式

3.键是隐式的并且由力场本身动态计算，有时可以使用“run 0”来调试初始设置

4. 有时首先使用 fix nvt 而不是 fix npt 以恒定体积平衡系统

七、ReaxFF的QEq电荷问题

#########ReaxFF必要时需要事先粗略分配电荷值，初始电荷分配预防或避免陷入局部最小值，让电荷非零的目的是打破对称性从而引导QEq计算器进入正确的计算方向，而初始给定电荷并不一定非常精确########

ReaxFF电荷采用QEq电荷，在迭代的每一步中进行更新。然而，在某些情况下，此过程可能会陷入非正常状态，而适当的初始猜测可以帮助获得正确的结果，并且显著加快第一步。

最简单的方法是从非极化经典力场（例如 CHARMM、Amber、OPLS 等）获取部分电荷。这个过程不必非常精确。关键点是打破其对称性，即确保原子具有预期的正电荷或负电荷。例如，对于水分子，氧原子应为负值，氢原子应为正值，但如果将它们分别设置为 -0.2 和 0.1，或者最初设置为 -0.8 和 0.4，应该不会产生太大差异。如果电荷值均为0.0，则 QEq 算法容易陷入错误的最小值。如果不用非极化经典力场，可以看看原子电负性，给所有比氢或碳负电性小的原子一个（小）正电荷，给所有电负性大的原子一个（小）负电荷，同时确保总电荷为0。

另外，非零电荷和零电荷有可能计算一样，但大多数体系可能不一样

非零初始电荷：

Step Temp Press KinEng PotEng TotEng

0 304.21579 25286.809 905.90267 -66506.359 -65600.456

10 1588.5172 15895.015 4730.3329 -70390.754 -65660.421

20 858.82628 11305.26 2557.438 -68289.665 -65732.227

30 701.53664 21764.819 2089.0563 -67835.757 -65746.7

40 1803.3259 21173.253 5369.9965 -71218.176 -65848.179

50 555.14885 14061.13 1653.1385 -67549.561 -65896.423

全零初始电荷：

Step Temp PressKinEng PotEng TotEng

0 304.21579 25286.809 905.90267 -66506.359 -65600.456