1.周期性边界条件

1.什么是周期性边界条件(PBC)

• 将小体积原胞平铺成无穷大的体相物质，$A(x)=A(x+nL)$
• 消除表面效应，构造出一个准无穷大的体积来更精确地代表宏观系统

2.周期性边界条件基本特点

• 平移对称性
• 边界固定，原子可通过边界发生相互作用
• 原子从模拟单元中跑出时，等效于从另一边界回来，体系原子数不变

3.最小镜像约定

• 周期性边界条件导致模拟体系有无限多个镜像，某个粒子所受到的力，理论上是所有镜像中粒子的作用力
• 最小镜像约定(Minimum Image Convention)：当将势场截断与周期性边界条件结合使用时，通常要求每个原子经历的净相互作用只包括来自系统中的一个(原子/镜像)
• $L<r_{cut}$，可能出现自己和自己的相互作用
  $L<2r_{cut}$，可能出现与另一个粒子相互作用计算两次的情形
  $L>2r_{cut}$，满足最小镜像约定，只可能与另外一个粒子(或者它的镜像)作用一次
• 分子动力学中，通常要求晶胞边长 $L>2r_{cut}$

4.Python实现

• Particle类：用于表示粒子，封装了粒子位置信息，并提供了获取和设置位置信息的接口

  class Particle:def __init__(self,position):self._position = np.array(position,dtype=float)@propertydef position(self):return self._position@position.setterdef position(self,value):self._position = np.array(value)
  

• Box类：用于表示盒子，封装了盒子的上界、下界以及长度信息，并可通过字符串方法获取盒子的长度

  class Box:def __init__(self,lower,upper):self.lower = np.array(lower,dtype=float)self.upper = np.array(upper,dtype=float)@propertydef length(self):return self._upper - self._lower
  

• BoundaryCondition基类：用于表示边界条件，继承该类的子类都必须实现apply和displacement函数

  class BoundaryCondition:def apply(self,particles,box):raise NotImplementedErrordef displacement(self,position1,position2,box):raise NotImplementedError
  

• PeriodicBoundaryCondition子类：用于表示周期性边界条件

  class PeriodicBoundaryCondition(BoundaryCondition):def apply(self,particles,box):for p in particles:p.position = (p.position - box.lower) % box.length + box.lowerdef displacement(self,position1,position2,box):dr = position1 - position2dr = dr - np.rint(dr / box.length) * box.lengthreturn dr
  

2.势场的有限距离截断

1.原子间相互作用力

$$V(R)=\sum _{i<j}v(\left|r_i-r_j\right|)$$

• 假设分子是刚性的，处于非简并的基态；原子间相互作用很弱，所以分子内部结构受环境影响较弱
• 原子核之间存在短程的排斥作用：exp(-r/c)
• 二体势主要特点：近距离处排斥，键长附近以及更远的距离相吸，实际应用中通常会取截取半径
• 未截断的计算复杂度：O(N²)，截断后计算复杂度：O(m×N)

2.势场截断的理论基础

• Lennard-Jones势(以及类似的二体势)随距离迅速衰减
      截断前：$V(r)=4\epsilon [(\frac{\sigma }{r}{)}^{12}-(\frac{\sigma }{r}{)}^6]$
      截断后：$V(r)=\{\begin{array}{cc}4\epsilon [(\frac{\sigma }{r}{)}^{12}-(\frac{\sigma }{r}{)}^6]& r\le r_c\\ 0& r>r_c\end{array}$
• 通过忽略截断半径以外的相互作用可以节省大量的计算时间
• 势能截断引起的误差项：
  $$U_{tail}=\frac{N\rho }{2}{\int }_{r_c}^{\infty }drV(r)4\pi r^2$$    当势场$V(r)$衰减慢于$r^{-3}$时(比如库伦势$V(r)\sim 1/r$)，截断引起的误差项将发散
  $$V(r)\sim \frac{1}{r^2}\to U_{tail}\sim r$$ $$V(r)\sim \frac{1}{r^4}\to U_{tail}\sim \frac{1}{r}$$
• 截断后的势场：$$V(r)=\{\begin{array}{cc}4\epsilon [(\frac{\sigma }{r}{)}^{12}-(\frac{\sigma }{r}{)}^6]& r\le r_c\\ 0& r>r_c\end{array}$$     在$r=r_c$处，$V(r)$不连续，受力$F=-\frac{dV(r)}{dr}$发散！分子动力学不稳定！

3.势场截断方法

$$V_0(r)=\{\begin{array}{cc}V(r)-V(r_c)& r<r_c\\ 0& r\ge r_c\end{array}$$

• 在$r=r_c$时，一阶导数不连续！即能量的一阶导数力不连续！
  同时移位后的势阱深度也偏离了原势函数，引起热力学性质的计算偏差
• 引入开关函数 $S(r)$：$\tilde{V}(r)=V(r)S(r)$
  对 $S(r)$ 的要求：$S(r_s)=1S^{{}^{\prime }}(r_s)=0S(r_c)=1S^{{}^{\prime }}(r_c)=0$
• cos开关函数：$$t=\frac{r_c-r}{r_c-r_s}=\{\begin{array}{cc}1& r=r_s\\ 0& r=r_c\end{array}$$ $$S(t)=\frac{1}{2}(1-cos\pi t)$$
• 多项式开关函数：$$t=\frac{r_c-r}{r_c-r_s}=\{\begin{array}{cc}1& r=r_s\\ 0& r=r_c\end{array}$$ $$S(t)=t^3(10-15t+6t^2)$$

3.近邻列表构筑与更新

• 原则上，原子间相互作用无论间隔有多大，原子对之间的作用具有非零值
  实际上，对于像Lennard-Jones这样的势函数，一旦间距超过 $\sigma$ 的几倍，原子间相互作用近似为零
• 所以，与给定原子存在相互作用的近邻原子通常有限，可用一个近邻列表来存储

1.近邻算法：VerletList法

• 对于每个原子，保留 $r_s=r_c+\delta$ 范围内的相邻原子列表。如果 $r_{ij}<r_c$，则计算力。每隔几步更新近邻列表
• 如何选取合适的 $\delta$：取决于体系的温度、扩散速度、密度等，通常约10-20个时间步更新
• 更新VerletList的判断条件：体系中原子的最大位移大于壳层厚度 $\delta$ 的一半
• VerletList类提供了获取和设置近邻列表的接口。build函数用来构建近邻列表

  class VerletList(NeighborList):def build(self):self.neighbor_list = {}self.previous._positions = {p: np.copy(p.position) for p in self.particles}for p1 in self.particles:self.neighbor_list[p1] = []for p2 in self.particles:displacement = self.boundary_condition.displacement(p1.position,p2.position,self.box)if p1 != p2 and np.linalg.norm(displacement) < (self.cutoff + self.skin_depth):self.neighbor_list[p1].append(p2)
  

• update函数用来判断并且决定是否需要更新近邻列表

  class VerletList(NeighborList):def update(self):max_displacement = max(np.linalg.norm(self.boundary_condition.displacement(p.position,self.previous_positions[p],self.box))for p in self.particles)if max_displacement > self.skin_depth / 2:self.build()
  

2.近邻算法：区间排序法(CellList)

• 计算近邻列表时，将框划分为单元格大小 $L=r_c$
• 构建列表时，只需要考虑粒子所处的26个相邻单元格
• 近邻列表的构造复杂度是O(N)

4.温度与压强控制算法

1.温度控制

• 体系的温度可以用粒子平均动量进行衡量，$T(t)=\frac{1}{N}{\sum }_{i=1}^N\frac{1}{2}m{v}_i(t)\cdot v_i(t)$
• 在分子动力学模拟的过程中，体系的温度会波动。为了实现对体系温度的控制，需要引入额外的运动控制，也就是热浴
• 温度调控机制可以使系统的温度维持在给定值，也可以根据外界环境的温度使系统温度发生涨落
• 一个合理的温控机制能够产生正确的统计系综，即调温后各粒子位移发生的概率可以满足统计力学法则

2.温度控制算法

• 随机方法：Anderson热浴，Langevin热浴，Dissipative Particle Dynamics
• 确定性方法：直接速度标定法，Berendsen温控，Nose-Hoover稳定
• 直接速度标定法，通过调整粒子的速度使系统温度维持在目标值
      引入速度标定因子$$T(t)=\sum _{i=1}^N\frac{m_i{v}_i^2(t)}{k_B{N}_f}$$ $$\lambda =\sqrt{\frac{T_c}{T(t)}}$$     每隔一定积分步，以 $\lambda v(t)$ 代替 $v(t)$ 进行速度标定，从而使系统温度在目标值附近小幅波动
      优点：原理简单，易于编程
      缺点：模拟系统无法和统计力学的系综对于起来，突然的速度标定引起体系能量的突然改变，与真实结构的平衡态相差较远
• Nose-Hoover机制，引入一个额外的“虚拟”自由度(热浴)，这个额外自由度与系统粒子相互作用，吸收或释放能量，使系统温度趋近于目标值
      Nose-Hoover热浴的基本思想是通过改变时间步长来调整系统中的粒子速度和平均动能，因此，在Nose-Hoover方法中，引入了新的变量s以重新调整时间单位
      将一个按照微观正则演化的虚拟系统映射到一个按照正则系综演化的实际物理系统
      证明了虚拟系统中的微观正则分布等同于真实系统中的(p’,r’)变量的正则分布

3.压强控制

• 分子动力学中，压强定义为在容器内壁受的单位面积上的平均力。压强的微观来源来自单位面积上的动量传输$$P=\frac{1}{3V}[\sum _i{m}_i{v}_i^2+\sum _i\sum _j{r}_{ij}\cdot f_{ij}]$$
• MD中，压强控制通常通过改变模拟盒子的尺寸或形状来实现压力的调节
• 常用的压力调节算法包括Berendsen压力耦合、Parrinello-Rahman方法、以及Nose-Hoover压力耦合等
• 要求在压力调节方向上表面具有周期性边界条件(PBC)，以减少表面效应并保证模拟的准确性
• 系统的原子或分子经常会重新排列以适应不同的压力条件，这种重新排列有时可以导致新的稳定或者亚稳定的晶体结构的形成
• 压力的控制没有考虑相互作用的局限性，具有一定的非物理的动力学特征

5.原子间作用势

• 原子间相互作用势：描述原子间相互作用关系的函数或模型，直接决定了原子模拟结果的可靠性

  1.Lennard-Jones势

  $$V(r)=4\epsilon [(\frac{\sigma }{r}{)}^{12}-(\frac{\sigma }{r}{)}^6]$$其中，$(\frac{\sigma }{r}{)}^{12}$为短程排斥作用的Pauli exclusion principle，$-(\frac{\sigma }{r}{)}^6$为远程吸引作用的London dispersion
• 原子间平衡间距
  $$\frac{\partial V}{\partial r}=4\epsilon [-\frac{12}{r}(\frac{\sigma }{r}{)}^{12}+\frac{6}{r}(\frac{\sigma }{r}{)}^6]$$ $$=48\frac{\epsilon }{r}[-(\frac{\sigma }{r}{)}^{12}+0.5(\frac{\sigma }{r}{)}^6]$$ $$当r=\sqrt[6]{2}\sigma 时,V_{min}(r)=-\epsilon$$

1.吸引项的物理解释

• 由两个相互作用的振子构成的系统的哈密顿量为 $$\mathcal{H}={\mathcal{H}}_1+{\mathcal{H}}_2+\Delta {\mathcal{H}}_{12}$$
• 微扰项是两个偶极子的相互作用 $$\Delta {\mathcal{H}}_{12}=\frac{{\vec{p}}_1\cdot {\vec{p}}_2-3(\vec{n}\cdot {\vec{p}}_1)(\vec{n}\cdot {\vec{p}}_2)}{r_{12}^3}$$
• 从一阶微扰理论，可以计算出能量的变化 $$\Delta E\simeq \sum _n\frac{|⟨{\psi }_n|\Delta {\mathcal{H}}_{12}|{\psi }_0⟩{|}^2}{E_0-E_n}\propto -\frac{1}{r_{12}^6}$$

2.斥能项的物理解释

• 原子间作用力：分子内化学键(键长项、键角项、扭转项等等)，分子间范德华力、氢键、极化作用等等
• 原子与原子间的斥力主要来源：
      范德华排斥力：当两个原子过于接近时，它们的电子云之间的排斥效应变得显著
      库伦力：电子与电子之间、质子与质子之间的排斥力
• 为什么是 $(\frac{1}{r}{)}^{12}$ 尚无科学解释

3.约化单位

• 对于Lennard-Jones势，通常会使用模型内在的长度和能量单位来刻画相关的物理量。用 $\sigma$ 来衡量长度，用 $\epsilon$ 来衡量能量，把这样的单位制叫做LJ unit $$r_i^{{}^{\prime }}=r_i/\sigma V^{{}^{\prime }}(r)=V(r)/\epsilon$$
• 原始公式 $$V(r)=4\epsilon [(\frac{\sigma }{r}{)}^{12}-(\frac{\sigma }{r}{)}^6]$$
• 约化公式 $$V^{{}^{\prime }}(r)=4[\frac{1}{r^{{}^{\prime }12}}-\frac{1}{r^{{}^{\prime }6}}]$$
• 国际单位运动方程 $$m\frac{d^2}{d{t}^2}{\mathbf{r}}_{\mathbf{i}}=\sum _{j}48{\epsilon }_0[0.5(\frac{\sigma }{r_{ij}}{)}^6-(\frac{\sigma }{r_{ij}}{)}^{12}]({\mathbf{r}}_{\mathbf{ij}}/r_{ij}^2)$$ $${\mathbf{r}}_{\mathbf{i}}^{{}^{\prime }}={\mathbf{r}}_{\mathbf{i}}/\sigma t^{{}^{\prime }}=t/\tau \tau =\sigma \sqrt{m/\epsilon }$$ $$\frac{dr}{dt}=\frac{dr}{d{r}^{{}^{\prime }}}\cdot \frac{d{r}^{{}^{\prime }}}{d{t}^{{}^{\prime }}}\cdot \frac{d{t}^{{}^{\prime }}}{dt}=(\frac{\sigma }{\tau })\frac{d{r}^{{}^{\prime }}}{d{t}^{{}^{\prime }}}$$ $$\frac{d^{2}r}{d{t}^2}=\frac{d}{dt}(\frac{dr}{dt})=\frac{d}{d{t}^{{}^{\prime }}}\cdot \frac{d{t}^{{}^{\prime }}}{dt}(\frac{dr}{dt})=(\frac{\sigma }{{\tau }^2})\frac{d^2{r}^{{}^{\prime }}}{d{t}^{{}^{\prime }2}}$$
• LJ单位运动方程 $$\frac{d^2}{d(t^{{}^{\prime }}{)}^2}{\mathbf{r}}_{\mathbf{i}}^{{}^{\prime }}=\sum _{j}24(r_{ij}^{-6}-2r_{ij}^{-12})({\mathbf{r}}_{\mathbf{ij}}^{{}^{\prime }}/r_{ij}^{{}^{\prime }2})$$     此微分方程与具体的惰性气体分子性质无关

2.其他二体势

1.简谐势

• 势能： $V(r)=\frac{1}{2}k(r-r_0{)}^2$
  相应的力： $\vec{F}(r)=(r-r_0)\frac{\vec{r}}{r}$
• 简谐势在描述分子振动时通常适用于小振幅振动，即原子围绕平衡位置进行小幅度的振动
• 利用简谐势可以描述两个共价键合原子的键伸缩

2.硬球模型

• 硬球模型将分子看作是无限小但具有硬核半径的硬球
• 分子之间的相互作用仅由排斥力来描述
  势能函数： $V(r)=\{\begin{array}{cc}0& r\ge 2R\\ +\infty & r<2R\end{array}$
• 硬球模型的势能函数通常为0；当两球相接触时，势能为无穷大

3.软球模型

• 实际上，当分子互相靠近时，电子云能够有一定程度上的重叠，因此互相可以有一定的穿透，此为软球模型
• 势能函数： $V(r)=\epsilon (\frac{\sigma }{r}{)}^n$
• 硬球模型和软球模型缺乏对分子间吸引作用的考虑，因此不能用于模拟真实流体，但可用于考察排斥作用对流体结构和液固平衡的影响

4.Buckingham势

$$U(r)=A{e}^{-Br}-\frac{C}{r^6}$$ 其中，$A{e}^{-Br}$ 为短程排斥力，$-\frac{C}{r^6}$ 为远程吸引力

• 与Lennard-Jones势相比，Buckingham势的排斥项更加灵活，贴近实际，但计算量也大很多
• 当 $r\to 0$ 时，因为指数项收敛于一个常数，因此，当势函数随着原子间距变小趋向发散时，Buckingham势会呈现不稳定的现象

5.Born-Mayer势

$$U(r)=A{e}^{-\beta r}-\frac{k{q}_i{q}_j}{r}$$ 其中，$A{e}^{-\beta r}$ 为短程排斥力，$-\frac{k{q}_i{q}_j}{r}$ 为库伦吸引力

• Born-Mayer势能用于描述晶体格子中的离子相互作用。与简单库伦模型相比，Born-Mayer模型是一个更现代、更准确的离子相互作用表示方式，特别是它考虑了由于排斥Pauli不相容原理而产生的短程排斥作用

6.Morse势

• Morse势是一种对于双原子分子间势能的简易解析模型 $$U(r)=-D_e+D_e[1-e^{-\alpha (r-r_e)}{]}^2$$
      如果对Morse势在 $r_e$ 附近作Taylor展开，可以得到 $$V(r)\approx \frac{1}{2}k(r-r_0{)}^2$$
• $D_e$ 是势能的深度，表示势能的最低点处的值
• $\alpha$ 是一个调节势能曲线陡峭程度的参数，可拟合材料的其他参数，比如内聚能、晶格常数
• $r_e$ 是势能曲线的平衡距离

7.库伦势

• 原子电子亲和力和电离电位的差值导致原子间的电荷转移
      静电能(假设点电荷)：$$V_{elec}=C_{unit}\sum _{i<j}\frac{q_i{q}_j}{R_{ij}}$$
• 库伦势的计算方法：
      使用周期性边界条件的情况下(无限系统)：
      库伦和是条件收敛的，结果取决于加和的顺序，结果取决于表面的电荷排列
      需要使用特殊的求和技术和边界条件：Ewald方法，Particle-mesh Ewald等技术
• 如何获得原子电荷？ 形式电荷、偶极矩或极化、电子结构计算、自洽电荷平衡

8.Yukawa势

• Yukawa势用以描述核子之间的短程相互作用 $$V(r)=-g^2\frac{e^{-amr}}{r}$$
• g：幅度缩放常数，即电势的振幅
• m：粒子的质量
• r：到粒子的径向距离
• 如果 $m=0$，则Yukawa势退化为库伦势 $V(r)=-g^2\frac{1}{r}-g^2=\frac{q_1{q}_2}{4\pi {\epsilon }_0}$

9.Aziz势

• Aziz势是描述低温和中等温度下惰性气体的相互作用，如氦、氖等 $$V(r)=\epsilon [A{e}^{-\alpha x}-F(x)\sum _{j=0}^2{C}_{2j+6}/x^{2j+6}]$$ $$F(x)=\{\begin{array}{cc}{e}^{-(\frac{D}{x}-1{)}^2}& x<D\\ 1& x\ge D\end{array}x=\frac{r}{r_m}$$
• Buckingham势和Aziz势的差异主要在于原子间距离较小的部分，Aziz势与Buckingham势相比具有更高的排斥力

10.Tang-Toennies势

• Tang和Toennies提出了Tang-Toennies势来描述稀有气体(He到Rn)之间的范德华相互作用 $$V(r)=A{e}^{-br}-\sum _{n=3}^N{f}_{2N}(bR)\frac{C_{2N}}{R^{2N}}$$     其中，$A{e}^{-br}$ 为分子间库伦排斥势，$-{\sum }_{n=3}^N{f}_{2N}(bR)\frac{C_{2N}}{R^{2N}}$ 为分子间相互作用的色散能和电子云重叠效应
• 衰减函数 $$f_{2N}(x)=1-e^{-x}\sum _{k=0}^{2n}\frac{x^k}{k!}$$

3.二体势的局限

1.对势的局限

• 通常一个体系的总能量可以根据个体、对势、三体势等的坐标划分为不同的贡献项：$$V=\sum _i{v}_1({\vec{r}}_i)+\sum _i\sum _{j>i}{v}_2({\vec{r}}_i,{\vec{r}}_j)+\sum _i\sum _{j>i}\sum _{k>j>i}{v}_3({\vec{r}}_i{\vec{r}}_j{\vec{r}}_k)+...$$
• 对势近似给出了对体系总能量非常好的描述，因为通过定义“有效”对势可以部分地包括平均三体效应。计算机模拟中采用的通常是这种有效对势，代表了部分多体效应 $$V\approx \sum _i{v}_1(\overrightarrow{r_i})+\sum _i\sum _{j>i}{v}_2^{eff}(r_{ij})$$
• 所以为了再现实验数据所需的有效对势可能依赖于材料的密度、温度等，而真正的原子间对势则不然
• 在实际材料体系中，原子间相互作用不仅仅取决于键长，而且其键能取决于化学环境
• 二体势的问题在于其本身的形式，而不在于其参数
• 对势模型：原子间的键能不受周围其他的成键情况；实际情况：当中心原子与更多原子成键时，平均键能降低 $$对势预测：E\propto Z$$ $$实际体系：E\propto \sqrt{Z}$$

2.空位形成能的对势计算

• 原子配位数为Z，仅限于最近邻原子作用 $$E_{total}(N)=\frac{1}{2}NZ{E}_{bond}$$
• 内聚能 / 结合能 $$E_{cohesive}(N)=\frac{E_{total}}{N}=\frac{1}{2}Z{E}_{bond}$$
• 空位形成能 $${\epsilon }_{vac}=E_{vac}(N-1)-\frac{N-1}{N}{E}_{total}(N)$$   代入 $$E_{total}(N)=\frac{1}{2}NZ{E}_{bond}$$ $$E_{vac}(N-1)=E_{total}(N)-Z{E}_{bond}$$   求得 $${\epsilon }_{vac}=-\frac{1}{2}Z{E}_{bond}$$
• 对势近似下，空位形成能数值上等于结合能(符号相反)，与对势的具体参数无关！
• 但实际材料体系中，空位形成能与结合能比值 ~ 1/3

	Ecoh(eV)	Evac(eV)	Evac/Ecoh
Al	3.39	0.75	0.22
Ni	3.516	1.6	0.46
Cu	3.615	1.2	0.33
Ag	4.086	1.15	0.28
Pt	3.924	1.4	0.36
Au	4.079	0.95	0.23

3.一维无穷长链的对势计算

• 相邻原子间距为a，以Lennard-Jones势为例 $$E=\frac{1}{2}\times 2\times 4\epsilon \sum _{n=1}[(\frac{\sigma }{an}{)}^{12}-(\frac{\sigma }{an}{)}^6]$$ $$=4\epsilon [1.00(\frac{\sigma }{a}{)}^{12}-1.02(\frac{\sigma }{a}{)}^6]$$
• 能量极值 $$a=(\frac{2}{1.02}{)}^{1/6}\sigma \approx 1.118\sigma$$     与平衡位置相比处于压缩状态$${\epsilon }_c=1.04\epsilon$$
• 原子表面受力：$$F=4\epsilon \sum _{n=1}^{\infty }[\frac{12}{r_n}(\frac{\sigma }{r_n}{)}^{12}-\frac{6}{r_n}(\frac{\sigma }{r_n}{)}^6]\hat{x}$$ $$=\frac{4\epsilon }{a}\sum _{n=1}^{\infty }[\frac{12}{n^{13}}(\frac{1.02}{2}{)}^2-\frac{6}{n^7}(\frac{1.02}{2})]\hat{x}$$ $$=0.1456\frac{\epsilon }{a}\hat{x}$$ $$r_n=na=n(\frac{2}{1.02}{)}^{1/6}\sigma$$ $$(\frac{\sigma }{r_n}{)}^6=(\frac{1.02}{2})\frac{1}{n^6}$$
• Riemann zeta函数 $$\sum _{n=1}^{\infty }\frac{1}{n^{13}}=1.00024$$ $$\sum _{n=1}^{\infty }\frac{1}{n^7}=1.00835$$
• 表面原子受到向外的净力，对势作用下将向表面外侧弛豫；而实际情况下，表面原子通常将向内侧弛豫

4.多体势与机器学习势

1.EAM势 MEAM势

• 原子的能量非线性的依赖于周围的近邻原子(原子数目与距离)，用“电子密度”作为对局域化学环境的一种衡量
• EAM势 $$Ecoh=\underbrace{\sum _i{F}_i({\rho }_i)}+\underbrace{\frac{1}{2}\sum _i\sum _{j\ne i}V(R_{ij})}$$ $${\rho }_i=\sum _{i\ne j}f(R_{ij})$$     由原子核在电子气的嵌入能和原子核之间的相互斥能组成，适用于：金属及合金原子间相互作用，如Fe、Cu、Ni、Pt、Au等
• EAM势中的电子密度不具有方向性
• Baskes, Michael I. “Modified embedded-atom potentials for cubic materials and impurities.” Physical review B 46.5 (1992): 2727.
• MEAM势与EAM势方法类似，但是电子密度项中保留了不同的 s，p，d 分量 $$Ecoh=\underbrace{\sum _i{F}_i({\rho }_i)}+\underbrace{\frac{1}{2}\sum _i\sum _{j\ne i}V(R_{ij})}$$ $${\rho }_i=\sum _{j,k}f(R_{ij})\cdot f(R_{ik})\cdot g(cos{\theta }_{ijk})$$     角度项对于过渡金属元素和共价体系尤为重要

2.分子力场

• 成键作用：成键作用通常代表共价键的性质，键的拉伸、弯折、扭转可以表示为胡克定律
• 非成键作用：非成键作用通常来自于静电作用、范德瓦尔斯作用等
• 化学键的作用包括二体、三体、四体相互作用 $$U=\frac{1}{2}\sum _{i,j}V(R_i,R_j)+\frac{1}{6}\sum _{i,j,k}V(R_i,R_j,R_k)+...$$
  二体相互作用：简谐势、Lennard-Jones势、Morse势、Buckingham势等等
  三体相互作用：$$E_{\theta }=\frac{1}{2}\frac{k_{ijk}}{si{n}^2{\theta }_{ijk}^0}[cos\theta -cos{\theta }_{ijk}^0{]}^2$$
  四体相互作用(二面角) $$E(\Phi )=k_{ijkl}[1-cos(n_{jk}(\Phi -{\Phi }_0))]$$
• 分子力场的能量：简谐势描述的键振动和弯曲+键的扭转弯曲+库仑力与范德瓦尔斯作用$$U=\underbrace{\sum _{bonds}\frac{1}{2}{k}_r(r_{ij}-r_0{)}^2+\sum _{angles}\frac{1}{2}{k}_{\theta }({\theta }_{ijk}-{\theta }_0{)}^2}+\underbrace{\sum _{torsions}\sum _n{k}_{\varphi ,n}[cos(n{\varphi }_{ijkl}+{\delta }_n)+1]}+\underbrace{\sum _{non-bondedpairs}[\frac{q_i{q}_j}{4\pi {\epsilon }_0{r}_{ij}}+\frac{A_{ij}}{r_{ij}^{12}}-\frac{B_{ij}}{r_{ij}^6}]}$$

3.反应力场

• 反应力场的特点：
      1.仅元素种类不足以确定原子间相互作用，例如：sp3 碳和 sp2 碳的行为不同
      2.传统的分子力场中，需要指定力场类型来建立原子的特性，在MD模拟之前确定力场类型和原子连接性，力场类型和原子间连接在运行期间保持固定
      3.反应力场(如Reaxff和REBO)纯粹根据原子位置描述相互作用，并允许原子改变配位和环境
• Stillinger-Weber势，广泛应用于半导体材料 $$U=\sum _i\sum _{j\ne i}{f}_2(r_{ij})+\sum _i\sum _{j\ne i}\sum _{k\ne i,k\ne j}{f}_3(r_{ij},r_{ik},cos{\theta }_{ijk})$$ $$f_2(r)=\{\begin{array}{cc}A(B{r}^{-p}-r^{-q})exp[(r-a{)}^{-1}]& r<a\\ 0& r\ge a\end{array}$$ $$f_3(r_i,r_j,r_k)=h(r_{ij},r_{ik},{\theta }_{jik})+h(r_{ji},r_{jk},{\theta }_{ijk})+h(r_{ki},r_{kj},{\theta }_{ikj})$$ $$h(r_{ij},r_{ik},{\theta }_{jik})=\lambda exp[\gamma (r_{ij}-a{)}^{-1}+\gamma (r_{ik}-a{)}^{-1}]\times (cos{\theta }_{jik}+\frac{1}{3}{)}^2$$
• Tersoff势：
      1987年，Biswas和Hamann发现SW势能无法描述Si的非四面体构型，Tersoff意识到三体势的本质缺陷，基于Abell的理论工作提出，Abell, G. C. “Empirical chemical pseudopotential theory of molecular and metallic bonding.” Physical Review B 31.10 (1985): 6184. $$E_b=\sum _i\sum _{j>i}[a_{ij}{V}^R(r_{ij})-b_{ij}{V}^A(r_{ij})]$$     $b_{ij}$ 为键级参数(bond order)，Si体系势能很大程度上取决于键级 $$b_{ij}=(1+{\beta }^n{\zeta }_{ij}^n{)}^{-1/2n}$$ $${\zeta }_{ij}=\sum _{k\ne i,j}{f}_c(r_{ik})g({\theta }_{ijk})exp[{\lambda }_3^3(r_{ij}-r_{ik}{)}^3]$$ $$g(\theta )=1+c^2/d^2-c^2/[d^2+(h-cos\theta {)}^2]$$     Tersoff势可以准确模拟多种共价体系的亚稳态结构，可扩展到SiC和SiO₂等体系，Tersoff, Jerry. “New empirical approach for the structure and energy of covalent systems.” Physical review B 37.12 (1988): 6991.
• 键级决定了两个原子之间的局部键能，原子环境➡键级➡共价键
      1.键长拉伸、角度和扭转项取决于所涉及的键级，例如区分 sp3 和 sp2碳
      2.随着键级变为零，所有涉及的键能项也趋于零，键级是所有作用项的前置因子
      3.对过成键/欠成键的惩罚，限制一个原子可以成键的数目
• REBO势和AIREBO势
      Tersoff potential考虑了碳的不同杂化形式，而在Si晶体的模拟上取得了巨大成功，但无法处理石墨的非局域共轭键。比如，它无法区分石墨中的C-C共轭键及分子中的C=C键，在模拟石墨以及碳氢有机系统时遇到了困难
      Brenner改进了 $b_{ij}$ 的解析形式，考虑了C-H相互作用，改善了对自由基的描述，开发了reactive empirical bond order(REBO)
      为了描述复杂的CH大分子，Stuart等人在原REBO势函数的基础上引入范德华力和四体扭转力，开发了Adaptive Intermolecular Reactive Empirical Bond Order(AIREBO)
      AIREBO解决了石墨的非局域共轭键的模拟问题，REBO家族被广泛用于碳纳米材料的模拟

4.机器学习势

• 半经验势特点：假设某种给定函数形式(如二体或三体)，能够找到实验中的一些数据，使用理论+模拟来确定参数
• 机器学习势特点：以机器学习模型代替解析的函数形式，使用第一性原理数据确定参数
• ML势准确可靠，但比从头计算要快得多