GULP(General Utility Lattice Program)
GULP最初设计的目的是拟合立场，现在已经逐步发展成为模拟凝聚态物质的通用代码，可以模拟无机固体、团聚体、缺陷、表面、界面以及聚合物等。

简介      
GULP可以对三维周期性固体，气相团簇和大体积材料中的孤立缺陷进行各种模拟。通过使用壳模型，GULP还可以处理分子固体和离子材料。与其它模拟程序不同的是，GULP对固体使用对称性从而简化输入，加速计算。这会有一个数量级的改善。GULP现在已包含在Materials Studio软件包中。
功能      
1. 体系的类型：三维大体积材料，团簇，缺陷
2. 能量最小化
3. 过渡态
4. 由用户指定体系的空间群（230种）
5. 晶体特性：弹性常数，体积模数，扬氏模数，泊松率，静电绝缘常数，高频绝缘常数，折射指数，压电常数，声子频率，声子态密度，声子态密度投影，声子散射曲线，k空间的Patterson对称性，零点振动能，熵（恒定体积），热容（恒定体积），Helmholtz自由能，静电势，电场，电场梯度
6. 缺陷：Mott-Littleton方法，缺陷能量，缺陷位移的过渡态，缺陷频率
7. 拟合。弹性常数，静电和高频介电常数，晶格能，压电常数，梯度，频率，电势和结构的经验拟合
8. 分子动力学(NVE)
9. 结构分析：键长，键角，扭转角，密度和单元体积。

这样对离子进行能量上分析时涉及组构单元之间的互作用有四类：
(a)     每个离子自身的电子壳（s）与实心（c）之间的互作用：


(b)     不同离子的电子壳、实心之间的Coulomb作用：

其中i≠j 寓意同一个离子自身的电子壳和实心之间不计Coulomb作用。
(c)     不同离子的电子壳之间的短程 non-Coulomb排斥作用，其中正、负离子之间的non-Coulomb作用按Born-Mayer势形式描述：


负、负离子之间的non-Coulomb作用按Buckingham势形式描述：

这里A,ρ和C 是需要程序拟合的势参数，其具有不唯一性和可转移性。另外，考虑到正离子的半径小，将正、负离子之间的non-Coulomb作用忽略。

d)     由于氧化物晶体中存在共价结合性质的互作用，所以要增添涉及三体的弯曲键（bond-bending）作用项：

这里 表示所考虑的三离子均处于平衡位时的夹角， 是强度系数。

在实际分子和凝聚态的计算模拟中，原子间相互作用势的选择

主要取决于所研究的具体体系。

1.对于分子晶体， Lennard-Jones 势往往是合适的有效势；
2.对于共价晶体，反映共价键结合的原子间相互作用多体势是非常有效的；
3.对于过渡金属，镶嵌原子方法是一个理想的方案；
4.对于离子键占主导地位的陶瓷材料和金属氧化物，人们基于壳层模型发展了相应的半经验势，如库仑＋ Buckingham 势。

6.对于有机和生物大分子，通常采用分子力场方法。
Two-body Interatomic Potentials
Functional forms for two-body interatomic potentials incorporated (where r represents the distance between two atoms i and j).

(2)     整数和分数离子电荷模型
在常规计算点缺陷形成能时，为纳入固体化学对Schottky和Frenkel 缺陷生成的热力学讨论程式，离子电荷取电子电量e的整数倍。然而从理论物理概念上除整数离子电荷模型外，还可采用离子携带非电子电量整数倍的分数离子电荷模型：这是因为晶体中离子之间的电子云是交叠的，即一个电子在两个相邻离子周围都有出现的可能，一个电子的电量可按处在两离子周围的概率比例分配给两离子，于是按这种考虑，归属各离子的电量就可能不再是简单的整数：其具体数值可通过第一性原理（ab-initio）计算晶体电荷密度得到。事实上在拟合晶体的动力学谱问题中分子离子电荷模型早已被前人应用了。至于当系统失去（得到）一个离子时，它是否带走（带来）分数个电子电荷的问题我们可以这样看：只要在形成点缺陷的平衡方程两端都是电中性的，就说明晶体中某一处出现正离子空位缺陷，相应在另一处出现正离子的间隙或替代位，两者有效电荷相反，总电量相同，维持电中性。尽管晶体内部离子间出现非整数的电荷迁移，但不影响全局电中性的结果，这种物理图象是可以接受的。

§3.2.2 Mott-Littleton法

根据前一节分析，由于缺陷问题所涉及自变量数为离子数(~150)×离子独立组分(2)×坐标自由度(3)~1000个，因此为了节约计算机时，著名的Mott-Littleton法[74]基于Born离子晶格模型，将缺陷周围的晶格简化处理为三个区域，分别记为I，IIa和IIb区，如图3.2[17,73-74]所示，其中区域I按原子水平处理，区域内所有离子都将因缺陷的存在而发生驰豫，直到它们重新达到平衡（即它们的力学和电学合相互作用为零）为止。而区域II距离缺陷较远，受到来自缺陷的作用力相对小一些，故晶格驰豫的处理采用了相对“经济”的算法。在实际处理中，IIa区域内的离子被假定处在一个谐振子势阱中，它们对缺陷引起的作用力产生响应。对IIb区域内的离子只考虑亚点阵而不是每个离子的极化。为了保证能量的收敛，区域I和II的半径要有一定的要求。在我们的处理中，根据对PWO晶体的模拟经验，设定区域I的半径为10.0 Å，区域II 的半径为16.5 Å, 因为此时计算出来的值与增大半径时计算出来的值几乎相等，即保证了计算的正确性。

§3.2.3 离子间互作用势的确定
对于经验途径计算方法一个重要的环节就是确定离子间互作用参数，在经验途径研究中对互作用势模型的讨论是唯象的。所谓经验参数化方法是通过尝试误差法来确定既定电子壳模型中的参数（Y，K等）和势参数表达式中所涉及的参数（如Buckingham函数中的A、ρ、C等），即合理的一组参数选取，应使基于它们所计算出来的晶体性质与相应的实验数据有最佳符合。           

令人满意的势模型应能准确复现晶体实验结构并能较好的解释或预言以下四个主要的物理性质：①晶格能，②弹性系数，③介电性质，④晶格振动。随着更多方面实验消息的介入，势模型也在不断的改进和发展。
为此可以根据理论计算值与实验值的偏差来建立定量评估所设置参数**理程度的目标函数；随后通过某种相应的方法势参数集优化。一般有两种方法：原先的常规拟合方法和新近发展的驰豫拟合（relaxed fitting）方法。
常规拟合经验参数方法确定经验参数的拟合方法所使用的评估参数集优化质量的目标函数可以写作:


其中{Ti}是当晶体处在X-Ray实验结构时的内应力及外应力的**，对内、外应力计算的项数随所研究晶体种类不同而变化，其中主要看晶体所归属的空间群及晶体原胞中的离子数，一般对高对称、原胞离子数少的晶体对应彼此不相关项的应力项数少。{PjExpr}、{ PjCalc }分别表示晶体物理性质的实验值和计算值的**，所以公式第二项也就表示具有代表性物理性质的计算值同实验值的相对误差。另外由于两项的量纲不统一，需要引入一个参数μ来使量纲一致，并通过调节μ的大小来控制拟合的相对权重。接下来就是需要以某个经验参数集作为起始点通过一定的优化途径使目标函数F达到极小得到最佳的参数**。
与常规拟合经验参数方法不同的是驰豫拟合参数法包含了两个步骤：第一步是结构优化，由经验参数集根据能量最小原则获取与之相应的晶格（驰豫低能）结构。添加这个步骤的合理性可以从两个方面来理解，一方面以实验结构为前提，某经验参数集所对应的晶格能可能恰好处在能量的鞍点上，这时尽管应力很小，晶格却处在不稳定点上，所以常规方法单纯从应力大小来判断经验参数集复现晶体实验结构的能力有时是不恰当的。另一方面，许多晶体性质（如介电、弹性、晶格能等）计算值的准确性却以晶体处在能量最小的稳定状态为前提条件，而常规拟合方法往往不满足这个条件。驰豫拟合的第二步是以第一步已优化的结构（或称“驰豫结构”）为前提计算物理性质，建立定量描述优化程度的目标函数：