本文则利用约束的分子动力学模拟来研究LiCl离子对和NaCl离子对在ABEEM-7P水中的缔合。我们选择LiCl离子对和NaCl离子对作为我们的研究对象是由于它们在化学和生化等方面的重要性。我们计算了Li+-Cl-离子对和Na+-Cl-离子对在水溶液中的平均力和平均力势，以及接触离子对解离的结构和动力学性质。
约束动力学(constraintdynamics)是对单独的内坐标或特定的坐标进行约束，即是在模拟的过程中固定这些值，但不影响其它内部自由度。使用约束动力学的优点在于当与刚性自由度相联系的高频率振动消除时，能在MD算法中使用一个较长的时间步长。但如果体系中存在大量的约束，将变得困难或甚至不能解析求解二次型的约束方程。因此在实践中，人们通常用一个较简单的迭代方法来满足约束，这一方法称为shake方法。在模拟中，在一个迭代循环中逐一处理所有的约束，然后重复这一过程直到所有的约束收敛至指定的精度。在模拟反应速度时，约束动力学用来模拟位于可将反应物和产物隔离的自由能垒的顶端的体系。另外，通常可以通过计算一个体系变换到另一个体系的可逆功来计算在一个体系的两种状态之间的自由能的差别。此外，约束动力学模拟还可以用来计算此可逆功。
在气态反应时，u(r)为反应物R+和X-直接相互作用势，如果反应发生在溶液中，则要加上由于溶剂作用的非直接的有效势能(r)，即在溶液中反应物R+和X-有效相互作用势为w(r)u(r)(r)。图8-28给出了高密度的溶剂对反应物R+和X-，以及过渡态RX的影响，可以看出，溶剂的作用对过渡态的影响是最大的。下面我们以一个离子对在水溶液中发生缔合反应来研究一下溶剂对反应的具体影响。
图8-28溶剂对离子缔合反应的影响的示意图我们考虑一个体系包含两个溶质R+和X-(R+=Li+或Na+以及X-=Cl-)，以及214个水分子，采用约束分子动力学模拟。在分子动力学模拟的过程中，保持LiCl离子对的间隔和NaCl离子对的间隔固定，因此，增加的约束是溶质R+与水分子之间的相互作用力RF以及溶质X-与水分子之间的相互作用力XF可以在模拟过程中进行计算。溶剂对溶质R+和X-之间的平均力的贡献(r)可以通过对不同样本进行平均来计算，根据下式r是沿着RX方向的单位矢量。如果令()dFr表示溶质-溶质间的直接相互作用力，则两个溶质R+和X-之间的平均力可以写为：
()()()dFrFrr（8-29）这样，离子对R+和X-在水溶液中的平均力势W(r)可以通过积分得到，即()0Wr的选择是非常重要的，其要能使得离子对处于较大间隔时可以获得可靠的平均力势。一般的，()0Wr应与宏观的库仑势相近，即0ij0W(r)qq/εr，其中，iq和jq为阳离子和阴离子的电荷，为水的介电常数。在本文，LiCl离子对选用的0r值为7.0埃，NaCl离子对选用的0r值为8.0埃。模拟过程中的一些其它细节与前几章模拟金属离子水溶液和卤素离子水溶液的模拟细节相同。
沿着选定坐标的自由能表面被称为平均力势（thefreeenergysurfacealongthechosencoordinateisknownasapotentialofmeanforce，PMF），对于离子对在水中发生缔合反应，作为离子对间距离的函数的离子对自由能（又称为平均力势），它的形状揭示了离子对溶剂化过程的很多有趣的细节，而且是与离子反应的许多特点密切相关的。我们利用约束的分子动力学模拟和标准的积分技术计算了LiCl离子对和NaCl离子对在ABEEM-7P水中的平均力势。在每一步模拟中，固定离子对的相对位置可以计算沿着离子对轴方向的平均力，其为离子对间距离的函数，积分平均力就会得到平均力势。在这些计算中，LiCl离子对间的距离被固定的范围为2.0~7.0埃，NaCl离子对间的距离被固定的范围为2.5~8.0埃，中间间隔为0.25埃。每次模拟的总时间至少为200ps，在一些重要区域，如接触离子对CIP和它的解离势垒区域，我们增加了模拟时间。表5.1总结了LiCl离子对和NaCl离子对在不同间隔时的平均力。
图8-29给出了LiCl离子对和NaCl离子对在ABEEM-7P水中的平均力势（PMF），从图中可以看出，LiCl离子缔合和解离的势垒分别为4.04Kcal/mol和2.05Kcal/mol，即LiCl离子对缔合势垒高于离子对解离的势垒，而且，溶剂分离的离子对SSIP区域宽于接触离子对CIP区域，也就是说，在LiCl水溶液中，溶剂分离的离子对比接触离子对更稳定。最近一篇文献[294]研究LiCl离子对在水溶液中的缔合情况，其确定的SSIP→CIP势垒为5.68Kcal/mol，而CIP→SSIP势垒为1.60Kcal/mol，与我们的结果略有差异，但其确定在LiCl水溶液中溶剂分离的离子对的布居明显大于接触离子对的布居，这一事实与我们的模拟结果完全一致。NaCl离子缔合和解离的势垒分别为1.76Kcal/mol和1.83Kcal/mol，而且溶剂分离的离子对SSIP区域宽于接触离子对CIP区域，可以看出，接触离子对与溶剂分离的离子对都占有一定的优势，存在着相互竞争。对比以前的模拟计算，关于LiCl水溶液和NaCl水溶液中接触离子对与溶剂分离的离子对的相对稳定性，人们的研究结果并不完全相同，可能是由于采用不同的势能模型。例如，Smith和Dang模拟NaCl离子对在可极化水中的缔合，他们给出接触离子对CIP解离的势垒是1.8Kcal/mol，溶剂分离的离子对SSIP的最低能量仅比接触离子对CIP的最低能量高0.2Kcal/mol，与我们的计算结果很相近。
图8-29LiCl离子对和NaCl离子对在ABEEM-7P中的平均力势我们再关注一下桥水分子，如果离子对间的距离增加，桥水分子取向放松，趋向于更加自由、单一配位的分布。此外，在这些模拟中桥水分子的平均数目（Nb>发生连续的振动。在LiCl离子对间的距离固定在3.0埃（rLiCl=3.0埃），NaCl离子对间的距离固定在3.75埃（rNaCl=3.75埃）时，我们计算了桥水分子的平均数目（Nb>与沿着离子对轴的力（Force）的相关性，展示在图8-30中。当水分子逐渐地插入到离子对中间的过程中，一个平均力使得这个离子对逐渐分离，因此，模拟过程中的桥水分子的平均数目（Nb>有非常大的波动。
图8-30桥水分子的数目与沿着离子对轴的力的关系图接触离子对CIP和溶剂分离的离子对SSIP构象的互变对研究溶液中的离子对反应的动力学是非常重要的。经典的过渡态理论给出接触离子对CIP解离的速率常数仅是温度、离子对约化质量以及CIP区的平均力势的函数。过渡态理论给出一个反应的正向反应速率常数为：
其中，r是活化势垒的位置，1B)(Tk，I是离子对约化质量。考虑到对正向反应速率常数TSTfk的校正，我们获得校正后的速率常数为TSTfk，这里是转换系数（transmissioncoefficient）。转换系数可以通过反应通量（reactiveflux）来获得，归一化的反应通量[296]为：
[x]是Heaviside函数，如果x大于0，其值是1，否则为0。平均值可以通过计算模拟轨迹得到（如来自1000个储藏轨迹），离子对的距离约束在r（LiCl离子对：
r=3.0埃；NaCl离子对：r=3.75埃），接着约束被移走，新的速度来自于Boltzmann分布，粒子向前运动1.5ps，向后运动1.5ps。转换系数还可以利用其它理论来计算，如Grote-Hynes转换系数(GH)[297]和Kramers转换系数(Kr)，分别根据下面两式计算。
在公式(8-33)中，bω是过渡态处的势垒虚频（theimaginaryfrequencyatthetransitionstate），(t)是施加于反应坐标方向上的取决于时间的磨擦（thetime-dependentfrictionexertedonthereactioncoordinate），是磨擦常数（frictionconstant）。我们根据振动力f(t，r)（fluctuatingforce）的时间相关函数C(t；r)来计算(t)和，根据下面的公式：
归一化的反应通量和归一化的振动力的时间相关函数展示在图8-31中。计算的过渡态理论的速率常数kTST、各种理论计算的转换系数以及各种理论校正后的速率常数k都列于表8-17中。对于LiCl离子对和NaCl离子对在水中的缔合，我们计算的过渡态理论的速率常数TSTfk分别为0.44ps-1和0.33ps-1，计算的转换系数分别为0.062～0.26和0.20～0.37，相应的速率常数TSTfk分别为0.027～0.11和0.066～0.12。从我们计算的结果可以看出，溶剂的作用在一定程度上减小了离子缔合反应的速率常数。其它的模拟利用不同的分子间相互作用势和计算方法得到略有不同的结果，例如，Rey等研究NaCl离子对在水中的缔合，给出过渡态理论的速率常数TSTfk为0.37ps-1，转换系数为0.22（根据反应通量方法）、0.22和0.37（根据Grote-Hynes理论）、0.08和0.35（根据Kramers理论）。
8.7金属离子水溶液性质
饱和的烷烃分子是有机分子中最简单、最基础的体系，力场的最初发展也是从这样的简单分子开始的。1976年，N.L.Allinger发表了MM1力场，1977年发表的文章就系统地阐述了用MM2力场计算有机分子构象和结构的结果和前景，其后MM（1-4）、CFF、MMF等力场对烷烃分子的构象、结构和光谱性质进行了系统的研究，从简单的对角（谐振子）函数力场到耦合项力场，目前他们着重于利用力场方法对光谱性质的研究。但是，利用极化力场对烷烃体系的研究，尤其是系统的计算还没有报道。杨等人利用ABEEM方法对烷烃等有机分子的电荷分布进行了计算，能够非常好地重复从头计算Mulliken布居分析电荷。为了初步探讨我们模型的浮动电荷静电势能否与分子力场其它势能函数恰当地结合，所以，我们利用ABEEM/MM浮动电荷力场初步地对烷烃分子的构象进行了研究。
在计算ABEEM/MM模型下烷烃分子的构象能时，我们采用了实验和从头计算的初始结构，我们在表8-18中列举了利用ABEEM/MM浮动电荷力场计算的13个烷烃分子的构象能结果，表8-19中详细地列举出了四个烷烃分子在我们模型下的结构信息。我们发现，从乙烷到丁烷，随着链的增长，实验测得绕C-C键转动的势能垒高度从2.88kcal/mol增加到4.89kcal/mol，我们给出的结果为从2.90kcal/mol到4.90kcal/mol，几乎与实验值相等。我们对18个构象能对进行了统计，ABEEM/MM力场给出总的均方根偏差（Root-Mean-Square-Deviation，RMSD）为0.21kcal/mol；AMBER、MM3、CHARMM(MSI)力场给出的相应值分别为：0.58kcal/mol、0.50kcal/mol、0.77kcal/mol。通过以上同实验结果的比较，结果表明，加入了浮动电荷静电势的ABEEM/MM力场给出了更加接近实验值的烷烃分子构象能结果。下面我们就对各个分子构象的详细结果进行分析。
表8-18ABEEM/MM和其它力场，以及实验得到的烷烃分子构象能8.7.1直接烷烃
乙烷：我们模型给出的反式结构的C-C键长为1.526，重叠式结构键长为1.528，实验值为1.534，从头计算键长为1.527，我们的结果略小于实验值；C-C-H键角为111.2，和实验值相等；绕C-C键转动的势能垒值为2.90kcal/mol，也与实验的2.88kcal/mol符合得很好。
丙烷：我们模型计算得到的C-C键长和C-C-C键角为1.528和111.8，实验测得的C-C-C键角为112.0；C-H键长平均值为1.112，实验测得的平均值为1.107（5）；我们计算得到的三个碳上的氢都处于交错式（conf1），一对相临两个碳上的氢处于交错式另外一对处于重叠式（conf2）与都处于重叠式（conf3）之间的势能差为3.35kcal/mol和3.64kcal/mol，实验给出的结果为3.3kcal/mol和3.9kcal/mol，前者符合很好，较高的转动能垒值略小于实验值。
顺丁烷：我们模型计算的反式丁烷的末端和中间C-C键长度为1.529和1.530，实验只测得末端的值为1.531(2)；C-C-C键角为112.0，略大于对应的实验值113.8(4)；ABEEM/MM计算得到的交错式丁烷的C-C-C-C的二面角为64.27。从头计算（MP2/6-31G）优化的最终角度为63.67；对丁烷的构象能量分析一向是非常重要的，实验的数据也是最充分的，图8-32是简单地描述丁烷绕中心C-C键旋转的势能曲线，在0~180区间内，ABEEM/MM力场给出了清晰的两个势能垒，最高的势能垒为在0左右的顺式构象，相对最低的势能点反式构象的高度为4.90kcal/mol，次之的是邻位重叠式构象势能垒，相对交错式构象的势能高度为3.99kcal/mol，实验高度为4.1kcal/mol。
