为了微观模拟体系能够反映宏观实验现象, 需要通过周期性边界条件对模拟对象体系进行周期性复制, 以避免在实际中并不存在的边缘效应(edge effects)。原则上,对于任何分子体系的理论研究都需要求解含时薛定谔方程。但在实际工作中,更关注的是原子核的运动轨迹,这样的轨迹可以利用波恩-奥本海默近似(Born-Oppenheimer approximation),通过求解经典力学运动方程获得。Alder和Wainwright曾表示计算机模拟实验会成为联系宏观实验现象和微观本质的重要桥梁,在他们首次进行分子动力学模拟实验之后10 年,法国物理学家Verlet提出了一套牛顿运动方程的积分算法,与此同时提出的还有另一套产生和记录成对近邻原子的算法,从而大大简化了原子间相互作用的计算。这两种算法至今仍以一些变形的形式在实践中被广泛应用[1,2]。
过去几十年开发了多种原子级模拟方法,包括晶格静力学、晶格动力学、蒙特卡罗和分子动力学等。其中,分子动力学特别适用于塑性变形的研究,它通过一些规定的原子间相互作用势函数的原子相互作用系统的牛顿方程的解,研究变形过程中的实时行为,并包括晶格的非简谐性、内应力的高度不均匀,以及系统的瞬态响应等方面的影响。
分子动力学主要依靠牛顿力学来模拟分子体系的运动,以在由分子体系的不同状态构成的系统中抽取样本,从而计算体系的构型积分,并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。它对原子核和电子构成的多体系统,求解运动方程,是一种能够解决大量原子组成的系统动力学问题的计算方法,不仅可以直接模拟物质的宏观演变特性,得出与试验结果相符合或相近的计算结果,还可以提供微观结构、粒子运动以及它们和宏观性质关系的明确图像,从而为新的理论和概念的发展提供有力的技术支撑。
分子动力学的对象是一个粒子系统,系统中的原子间的相互作用用势函数来描述,因此,正确选择势函数的类型及其参数,对于模拟的结果优劣具有重要作用。势能函数在大多数情况将描述分子的几何形变最大程度地简化为仅仅使用简谐项和三角函数来实现;而非键合原子之间的相互作用,则只采用库仑相互作用和Lennard-Jones 势相结合来描述。其中,对于原子间相互作用力的描述通常是经验或半经验的,这样虽然能够提高计算效率,但无法完全揭示电子键合的多体性质,尤其对于缺陷附近与自身结构和化学性有关的复杂自洽变分函数。Daw和Baskws的EAM(Embedded-atom model)势函数在某种程度上融合了电子键合的多体性质,将系统的总势能表示为:
式中:Fi是原子i的嵌入能函数;ρi是除第i个原子以外所有原子在i处产生的电子云密度之和;Φij是第i个原子与第j个原子之间的对势作用函数;rij是第i个原子与第j个原子之间的距离[1]。
势函数的可靠性主要取决于力场参数的准确性,而力场参数可以通过拟合实验观测数据和量子力学从头算数据得到。目前在生物大分子体系模拟中使用最为广泛的分子力场是CHARMM力场和AMBER力场,是早期研究生物大分子的分子力场。现有的力场参数仍在不断优化之中,并且涵盖的分子类型也在不断扩大。粗粒化(coarse-grained)模型在计算生物物理研究中越来越引起人们的关注,由于在该模型中定义了粗粒化粒子,对应于全原子模型中的若干原子或原子基团甚至分子,减少了体系中的粒子数,使得模拟的时间和空间尺度得以大幅度提高,但同时也将丢失原子细节信息。基于这种模型的分子动力学模拟适合于研究缓慢的生物现象或者依赖于大组装体的生物现象。
设计一个基本力场的根本原则就是要单位时间步长(time step)内计算能量的开销最小化,从而实现模拟尺度的最大化。这一点对于全原子力场尤为重要, 即便对于所谓的粗粒化模型也是一样。特别地,如果想要进行微秒甚至毫秒级时间尺度的模拟,这一原则极其重要。
图1显示了分子动力学的时间和空间尺寸的反比关系,图中从左至右依次为:(1)水,细胞的基本成分;(2)牛胰蛋白酶抑制体,一种酶,其“呼吸”行为可以在毫秒级时间尺度上进行考察;(3)核糖体,一个复杂的生物器件,可以对基因信息进行解码并生产蛋白质;(4)紫细菌光合膜片段,拥有2500 万个原子,图中显示的是嵌在磷脂双分子层上的捕光复合物及光化学反应中心[2]。
图1 经典分子动力学的时间的尺度关系
随着计算机处理器速度的快速增长以及大规模并行计算架构的发展,大规模并行化或专用的架构技术与可扩展分子动力学程序的结合,计算机模拟包含从位错到基于晶界的变形机制的整个晶粒尺寸范围,为探索材料体系的研究前沿领域开辟了新的途径。
例如,William Gonçalves等人通过选用Wolf BKS(van Beest, Kramer and van Santen)势函数来描述原子之间的相互作用力,使用大规模原子/分子并行模拟器LAMMPS(Large-scale Atomic/Molecular Massively Parallel Simulator)对二氧化硅气凝胶的弹性和强度进行了分子动力学方面的研究,他们使用velocity-Verlet算法和1.0 fs时间步长,并且在三个方向上均使用周期性边界条件[3]。
图2为模拟的超过7000000个原子的大体积样品3D示意图,以及20nm厚的样品切片和局部放大图(蓝色为氧原子,红色为硅原子),图3(a)是对803nm3气凝胶样品进行单轴拉伸试验,以获得300K的应力应变曲线,(b-d)是典型韧性断裂图像,(e)抗拉强度与样品体积的对数关系。他们分析认为,为了确保正确评估弹性等机械性能,模拟样品的尺寸至少为孔径的8倍,同时,表面积极高的二氧化硅气凝胶需要相对低的应变率以确保准静态条件。