软物质的分子建模与模拟

作者简介

内容简介

第 1 章 绪 论
More is different！(多者异也！)
Philip Warren Anderson (菲利普 沃伦 安德森)
1.1 软物质的定义和特征
　　目前还缺乏关于软物质公认的准确定义，但它的特征是明确的。软物质可以是既不属于简单液体又不属于晶体的凝聚态物质的统称，其各项性质也正好介于简单液体和晶体之间。典型的软物质包括胶体溶液、气溶胶、高分子、乳液、玻璃、凝胶、液晶、生物聚合物、生物机器等。相对于硬物质 (原子排列整齐的固体——晶体) 而言，软物质中多项相互作用有类似强度，它们互相以及和热运动竞争，导致非常丰富的软物质性质、特征以及对外界条件的敏感性依赖。在软物质中，相比能量，熵效应通常更加重要，从而使软物质有较大的构象涨落，通常表现在介观 (介于微观和宏观之间) 的时间、空间尺度 (如毫秒、微米量级)。正是这些特征使得软物质得到了广泛的应用。
　　软物质大多对分子结构和相互作用的微观细节非常敏感，因此基于连续介质假设或平均场近似等的宏观理论方法在软物质的研究中的作用是有限的，而可以较精确刻画分子结构和相互作用的分子建模和模拟方法成为软物质理论研究中必不可少的工具。另外，软物质通常又具有介观的时间和空间特征尺度，往往超出了常规的分子建模和模拟所能达到的时空范围，因此必须对常规的分子建模和模拟方法进行扩展和改进，通过保留软物质特定性质相关的最少必要和结构相互作用细节并简化近似其余，才能有效应用于研究软物质的性质。这种近似通常需要同时考虑多个不同的时间和空间尺度特征，相应的方法称为多尺度建模和模拟。
　　自由能计算是软物质理论研究的关键。软物质自由能的状态及对外部条件的依赖决定其结构、性质和行为。将分子建模和模拟方法应用于软物质研究的重要目标之一也是计算软物质系统的自由能，包括计算系统自由能对热力学条件的依赖性，软物质不同状态、结构间的自由能差等。由于软物质的性质通常由熵效应主导，软物质在不同条件及状态下的自由能差别一般比较小，因此较准确的自由能差值的计算对于软物质分子建模和模拟显得尤其重要。
1.2 软物质模型和模拟
　　与其他宏观系统相同，对软物质的计算研究原则上可以从基于量子力学的第一性原理 (量子力学) 出发。第一性原理计算直接从组成系统的所有原子核和电子的相互作用出发，根据量子层面的物理学基本原理 (如薛定谔方程)，使用计算机数值计算系统的结构、性质以及演化。第一性原理计算的优点在于方法普适、结果较精确，但是受限于现有的计算能力，目前第一性原理计算通常只能对小系统 (几十到几百个原子) 计算较短时间 (皮秒量级) 的演化特征。这对于研究作为多体分子体系的软物质系统在有限温度下的统计性质，特别是结构和动力学性质，往往是远远不够的。因此，基于在较小系统中使用第一性原理计算的数据，结合实验测量，人们建立了全原子模型以及相应的分子模拟方法。该方法忽略了原子内部的自由度，把每个原子看作一个质点，考虑原子间经验力场的相互作用所决定的体系中所有原子的行为，使用统计物理原理得到系统 (在有限温度下) 的统计性质及动力学特征。
　　全原子分子模拟方法主要分为蒙特卡罗模拟方法和分子动力学模拟方法两大类。蒙特卡罗模拟方法基于平衡态统计物理理论，使用随机过程作满足玻尔兹曼统计分布的重要性采样，近似求解构象函数的系综平均值。其主要优点是实现简单、计算量小，并且不受积分步长的限制；缺点是原则上只能模拟平衡态体系，并且得不到真实的动力学信息。分子动力学模拟方法则通过求解多体牛顿运动方程模拟真实多体分子体系的时间演化。其优点是既可以模拟平衡态体系，也可以模拟非平衡态过程，而且不但可以获得结构信息，还可以获得动力学信息；缺点则是受积分步长的限制 (一般几飞秒)，可达到的模拟时间一般小于微秒量级。另外限于计算机的计算能力，目前全原子分子模拟可以研究的体系空间尺度最多到几十纳米。
　　为了研究较大的软物质系统，通常对全原子模型作更进一步的近似处理，即粗粒化建模和模拟。在粗粒化模型中，系统中紧密结合的一些原子团用单个粗粒化粒子表示，通过建立粗粒化粒子间的近似 (平均) 相互作用，即粗粒化力场，人们可以计算更大的系统。同时，因为粗粒化粒子间的相互作用一般比原子间相互作用更加平滑，特征时间尺度较长，所以分子动力学模拟可以使用更长的积分时间步长，以达到更长的模拟时间尺度。从全原子模型简化到粗粒化模型的思路与由第一性原理模型到全原子模型类似，都是通过忽略一些不太重要的内部自由度，使用近似的相互作用演化系统的部分特征。这种粗粒化程序可以在不同的层面实施，得到系统在不同时间和空间尺度的性质。需要指出的是，通过仔细构造力场，粗粒化建模和模拟一般在重建系统的较大尺度的静态性质和结构方面能给出较好的结果，但在描述系统的动力学特征等方面，随着粗粒化程度的增加，准确程度一般呈减少趋势，需要越来越多的努力以抓住系统的主要特征。
　　在软物质系统中，由于分子间相互作用较弱，往往显现出不同层级间关联较强的特征，例如软物质自组装往往是关联较强的多层级、多步骤的过程，而很多生物蛋白酶催化位点的作用与蛋白酶本身的结构有很大关系。这种不同时间和空间尺度间的耦合是软物质建模和模拟的中心难点，使得通常简单的粗粒化建模等方法不够充分。模拟研究这类软物质问题时，人们逐渐发展了多尺度建模和模拟的思想，即将不同层级、不同尺度的模型结合。一个较常用的例子是QM/MM(quantum mechanics/molecular mechanics，量子–分子力学杂化方法) 和全原子/粗粒化杂化模型，这类多尺度计算模拟方法是解决软物质中许多问题不可或缺的手段。QM/MM 对于需要电子自由度的部分 (比如蛋白酶的催化反应位点)进行第一性原理计算，其他部分则用全原子分子力场计算或模拟。全原子/粗粒化方法则是对需要全原子分辨率的部分进行全原子模拟，其他部分用粗粒化力场进行模拟。这几个较成熟的杂化模型方法取得了较大的成功，但限制也非常明显，它们要求预先设置系统中各个部分的分辨率，且在模拟过程中不再改变，这仅对一些特殊问题适用。对许多软物质中的问题，由于体系活性位点随时间变化，所以不能预先指定。例如，当研究蛋白质构象变化时，蛋白质分子中特殊位点与其近邻水分子间相互作用细节一般要求全原子尺度的描述，而远离这些位点的水分子可用粗粒化模拟描述。但水分子运动改变了它与蛋白质位点的距离，因此在不同时刻需要不同的模型。目前，一些采用自适应的杂化模型，在模拟过程中改变体系不同部分的模型精度的研究，取得了初步成果。发展合适的多尺度模型对模拟研究软物质系统至关重要，关于这方面的研究目前还比较初步，需要做大量工作。
1.3 增强抽样和加速模拟
　　计算宏观系统的自由能是对其进行理论和模拟研究的核心。目前已经有多种基于分子模拟计算系统自由能的方法。统计物理告诉我们，系统自由能 F 对其哈密顿量 H 中所含参数 λ 的导数等于哈密顿量对此参数导数在相应系统中的热力学系综平均值。因此，为了得到系统取两个不同参数值时的自由能差，可以计算这个热力学系综平均值在此参数区间的积分，此即热力学积分方法。
　　软物质体系沿一个或几个集体变量 (简记为 x) 的自由能变化曲面 F(x) 通常能对其许多性质，如亚稳态稳定性、过渡态的初步确定、亚稳态间转变反应动力学等提供直接刻画。原则上，分子模拟样本的分布 P(x) 直接给出这个自由能曲面 (仅差一个常数，对应自由能零点选取的任意性)，F (x) = c.kBT ln P (x) 。这里 kB 是玻尔兹曼常量 (Boltzmann constant)，T 是体系温度。kBT 是热力学下自由能的自然单位。实际上，由于通常模拟能给出的这个分布的数值范围不超过几个量级，从而仅能得到变化范围很小 (几个 kBT) 的自由能区间，不能满足要求。
　　这实际上就是使用分子模拟方法时常常遇到的模拟产生的样本不足的困难，与前文谈到的分子模拟能达到的时间尺度在研究许多软物质问题时不够长相类似，是目前应用分子模拟研究软物质问题的主要障碍。近几十年发展出了许多新的模拟方法，试图克服这个限制，以得到更大参数范围内的自由能曲面。这些方法的具体实现可能看起来十分不同，但它们的基本原理大多类似，即修改体系的哈密顿量，H (x) → H (x).ΔU (x) ，使得修改后的系统的自由能在整个 x 空间的变化范围较小，能被分子模拟充分访问从而较准确估计，而原始系统的自由能则通过所使用的有偏势重建，显然，如果选择偏置势 ΔU (x) 近似等于 F(x) (差别不大于几个 kBT 时)，则分子模拟修改后系统的所有 x 处都能充分访问，从而可以得到较准确的样本分布eP (x)，以给出对原系统自由能 F(x) 的较准确估计。在实际模拟中，由于 F(x)未知，人们一般使用迭代方法设置有偏势 ΔU (x)，使其逐渐接近 F(x)，在此过程中模拟抽样的样本所能达到的范围也逐渐扩大，直到所有的 x 都有较充分的抽样。
　　目前，依赖于不同的具实现 (主要是迭代设置有偏势方法) 发展出了许多模拟方法，如多正则系综模拟、Wang-Landau 算法、tadynamics 模拟、伞形抽样 (umbrella sampling) 方法、Hyperdynamics 方法等。另外，如果不在 x 空间使用有偏势，还可以对其热力学对偶控制参数 (如能量的对偶参数是温度) 进行调控来增加分子模拟在 x 空间的访问范围，同样能增强模拟抽样的效率，更有效地计算自由能曲面。这类方法包括副本交换、温度积分、温度行走模拟等，可以统称为广义系综模拟，即用其他分布代替标准的玻尔兹曼分布进行构象空间的抽样模拟，以扩展分子模拟在构象空间的访问范围，而系统满足玻尔兹曼分布的统计平衡性质可以从这个广义系综重建得到。
　　目前，这类增强抽样分子模拟方法取得了很大成功，但是也仍然存在一些不足，其主要困难是 (集体) 变量 x 必须是慢变量，而其他独立于 x 的变量是快变量，而在实际系统中慢变量的选择本身是一个困难问题。目前，这些不同的增强抽样模拟及自由能计算方法在具体问题中的使用各有优劣，一般根据问题特征选择合适的模拟方法克服困难。
　　在增强抽样模拟方法中，有一类不仅要求能从改变后的系统模拟中重建原始系统的热力学 (静态) 性质，还要求能够重建系统的动力学 (dynamics)，至少是重建其动理学 (kinetics) 行为。这类方法的一个典型是超动力学 (hyperdynamics)模拟，它抬高构象空间中那些能被充分访问的区域 (势能面盆地，更准确地说是亚稳态) 的能量，但不改变那些被较少访问 (势能面鞍点，或者说过渡态) 处的能量。这样系统动力学的关键部分 (其在过渡态区部分) 不被有偏势影响，而其在亚稳态区的很长时间的等待则被此有偏势统计地缩短，而这种时间缩短效应可以用亚稳态内的局部平衡假设理解与消除。超动力学方法在如固体表面原子迁移等小系统中得到了很好的应用，但是对一般系统，识别亚稳态与过渡态从而设计合适的有偏势是很困难的，从而不得不使用一些近似处理，所得结果也不再能保证重建体系的动力学特征，不过仍对重建系统的平衡性质有帮助。对于软物质，归因于系统中更为主导的熵效应，超动力学的有效使用需要更多努力，如需要对特定系统有微观层面的详细理解，以明确构建合理的有偏势。一般而言，这种理解的唯一来源正是分子模拟。在此基础上可以设计迭代方案，从分子模拟得到的数据出发，分析得到系统的亚稳态与过渡态，从而设置超动力学有偏势，扩展模拟的范围，然后再分析更多的亚稳态并设置新的有偏势。重复这个过程多次，有望扩展到需要的构象空间范围及想要达到的时间尺度。