这篇继续讲古典物理、热力学和动力学模拟方法 首先我推荐想钻研molecular modeling，尤其是生物分子模拟方法的人 阅读这篇历史回顾文章： Pnina Dauber-Osguthorpe & A. T. Hagler, "Biomolecular force fields: where have we been, where are we now, where do we need to go and how do we get there?" Journal of Computer-Aided Molecular Design (2019) 33:133–203 https://doi.org/10.1007/s10822-018-0111-4 纯论物理理论的艰深程度，上一篇讲的各种量子物理方法都比这一篇要讲的有格调得多 古典分子模拟方法的原理说穿了只有两点：牛顿力学、玻兹曼分布，没了 但是卖量子模拟方法全家桶的Gaussian Inc.一年的营收据说是$4M（当然是镁） 卖分子动力学模拟的Schrödinger, Inc. 2020 Q4的营收则是...$33M 其中$25M来自纯粹卖软体（这家公司有上市，看得到财报，Nasdaq代码SDGR） 巨大的市场差异说明了需求所在：量子力学模拟实在太慢了，也跑不了太大的系统 大家想要模拟整条蛋白质、想要模拟整条染色体，还有神经病想模拟整个细胞 量子力学模拟说准是准（其实也常常不准），但算到宇宙终结都算不出个东西来，没用 最後大家还是充满怨气回来使用平凡无奇、拼凑、毫无美感的古典分子模拟工具 前言太长，以下正文 4.古典分子模拟，或者说分子建模，或者说force field 总之建立这类工具的想法可以追溯到量子力学出现之前的十九世纪末 当时的人对於不同化学分子有不同的构型有着模糊的概念 开始思考怎麽把这些分子「具现化」到人类能凭感官直接分析的尺度来做研究 当然这些要认真说都算是穿凿附会 真正的force field大约是1930年代有量子力学支持的光谱分析发展起来後的需求 分子的放射频谱透露出分子内部原子间的键结如何震荡、摆动，分子如何旋转 化学家自然开始想到：能不能根据这些资讯，建立一个分子模型 让我们能随意调整摆弄这个模型，了解分子结构变化时的各种性质？ 在那个年代，用量子力学方法来模拟分子当然是不可能的 所以建立在古典力学之上的force field就出现了 基本上就是把原子当成小球，化学键则是弹簧或杆子（看简化程度） 化学键的键脚、扭矩、弹性系数等等由光谱实验提供 再加上原子之间的静电力和凡德瓦力 就完成了一个建立在实验数据上的位能函数，也就是force field 这个位能函数的自变数就是每个原子的位置向量 决定了这个分子系统的所有力学性质 想要知道分子基态结构？最小化位能函数取得能量最低点的个原子位置 想要知道分子如何从激发态结构开始振动？给定初始结构然後用F=ma作动力学模拟 想要知道分子在不同温度压力以下的表现？连上thermo/barostat做Langevin dynamics 机器不够不想管动力学只做热力学？用波兹曼分布做蒙地卡罗搞定 整个古典模拟领域，无论你是做新材料研发、新药物研发、溶液系统改良 生物物理、生化、物化、化工、生技、药物动力学 这麽多领域里，只要应用的工具是force field，那就只是F=ma加上P=exp(-U/kT) 不同的领域不同的原子分子造成的差别就只是如何写下最适合你的force field 因为force field只是个effetive/empirical theory，不是真正的基础理论 它的参数来自实验或者量子力学模拟，甚至是低一阶的古典分子模拟 这就是所谓的coarse graining，它的hierarchy大致上是： 量子力学方法（前一篇文章任一种） ->全原子力场 ->复合原子力场 ->单体高分子力场 ------以下属於发废文用工具------ ->细胞力场 ->组织/水溶液扩散方程式 ->...... 4-1.全原子力场就是每个原子都被化简成独立小球，有各自独立的参数 有些力场宣称自己能够处理任何分子，这叫做generalized force field (GFF)， 有些力场则是针对特殊系统调整参数，可以被称作xxx force field 在这类特化力场里，同类原子在不同位置可以有完全不同的大小和带电量 例如在蛋白质力场里，C-alpha和C-beta的参数就不一样 甚至在不同氨基酸里的C-alpha的参数都可以不一样 由於有效理论的局限性，许多GFF也会进一步区分不同位置的同类原子 比如说AMBER通用力场里sp3和sp2混成轨域的碳原子就被算作是不同种类的原子 实际分类更复杂，只包含COHN四种原子的AMBER版本大概有几十种不同的「原子」 目前通行最广的全原子力场有三家： - 由Martin Karplus (2013诺贝尔化学奖)建立的CHARMM (Chemistry at HARvard Molecular Mechanics) https://www.charmm.org - 由Paul Weiner和Peter Kollman等人建立的AMBER (Assisted Model Building with Energy Refinement) https://ambermd.org - 由Wilfred van Gunsteren和Herman Berendsen等人建立的GROMACS (GROningen Molecular Simulation) https://www.gromacs.org 三家的原始程式码和物理模型都可以追溯到更早的Shneior Lifson的CFF上 (Consistent Force Field) 之间的关系错综复杂，恩怨情仇可以写个好几十集的连载小说 不过我也不真的算很圈里的人，所以不要问我详情 4-2.复合原子模型是全原子模型的简化版 早期电脑太弱，在Karplus团队拿原始版CHARMM模拟蛋白质时根本跑不动 於是他们把氢原子省略，当成其相邻重原子(C/N/O)的一部分， 并对参数作相应调整才成功搞定 现在电脑够强了，所以复合原子模型已经很少用了 4-3.单体高分子力场通常是化工和材料系拿来模拟多分子系统用的 比如说化工研究高分子溶液系统，需要同时模拟上百甚至上千条高分子 这时就不太可能使用全原子模拟，必须再做简化 最後常用的模型就是把高分子的每个单体都当成一个小球处理 小球之间用弹簧连接，再加上静电力凡德瓦力等相应作用力 最後各项参数用实验数据调整 最近也流行先跑一个小的全原子模拟，用全原子的数据来调整单体模型的参数 这有个很fancy的词叫multiscale modeling，说穿了也没什麽 现在（北美）化工界常用的这类力场的开发工具是HOOMD-blue https://glotzerlab.engin.umich.edu/hoomd-blue/ 说「开发工具」而不说是一个力场 是因为一般说的force field还包含了某套特定参数 而HOOMD-blue是不给出参数的，毕竟每个人想研究的系统都不一样啊 它只是架好了一个小球-弹簧的模型 弹簧的长度、弹性系数、小球体积，甚至小球之间的位能函数都能自由定义 4-4.更加coarse grained的模型就太广泛了， 任何解不出/懒得想/懒得解物理方程所以搞出来的简化模型都算 这里只说一个比较常见的通用手段：lattice space 这个lattice不是固态物理晶体模拟的意思，而是把整个空间切块！ 原本分子可以出现在空间的任意位置，但这里强迫分子只能占据某个晶格 於是系统大幅简化为在晶格空间中「推箱子」，计算速度大为提升 只是如此粗粒化後的模拟结果是否依然有效就看你怎麽解释了 除此之外，晶格系统的很多物理特性还能够用纯解析方法推出 比如说，由於lattice gas system其实就是一个总磁矩固定的Ising model， 很多相变的critical exponent因此是相通的 Pierre-Gilles de Gennes (1991年诺贝尔物理奖) 也很常用晶格模型研究高分子 搭配重整化群得出了很多高分子相变的重要结论 後来又由此衍生出了高分子场论这门奇怪的学科则是另一个故事了... （另外要提一下，de Gennes 2007年癌症过世时是倒在他家的书桌前的， 桌上是正在写作的论文和计算草稿，真的是字面上的一辈子做物理了） 5.这边提一下古典分子模拟的几种不同的方法 其中有很多值得深究的数学物理问题（毕竟也算是古典力学的分支） 推荐一本书： Molecular Dynamics: With Deterministic and Stochastic Numerical Methods by Ben Leimkuhler & Charles Matthews 5-1.NVE模拟，即系统的粒子数(N)体积(V)和总能量(E)守恒 整个系统严格依照F=ma运动，又称Hamiltonian模拟 力场给定了U，系统给定了各粒子的m，各粒子的a就可由-grad(U)/m得出 问题是，得到a之後如何得出速度v和位置x？ 这就变成了一个数值积分问题，求解的演算法被称为integrator 最无脑的方法当然是：求解系统由时间T到T+t(t为某微小时间间隔)的变化 则v(T+t)=v(T)+a(T)*t, x(T+t)=x(T)+v(T)*t+a(T)*t^2/2 这就是Euler method 这个方法有许多缺点 首先，他的Global truncation error是O(t)，太差了 其次，Euler method并不是symplectic 大致上可以理解为，如果把x(T)->x(T+t)的演化当成一个空间变换 这个变化并不具有某些很好的守恒性质 後一点可以靠某些变换式来解决，但是O(t)的误差实在是没办法 所以Euler method基本上没人用（除非系统很小或者误差非常不重要） 最被广泛使用的integrator是Verlet/velocity-Verlet method 它的逻辑是： (1)由x(T)推出U(x(T))进而得到a(T)，然後计算v(T+t/2) = v(T) + a(T)*t/2 (2)计算x(T+t) = x(T) + v(T+t/2)*t (3)由x(T+t)推出a(T+t)，然後计算v(T+t) = v(T+t/2) + a(T+t)*t/2 这个方法满足symplectic，并且误差是O(t^2) 所需计算量又事实上跟Euler method一致 (第三步获得的a(T+t)可以直接给T+t的第一步使用，不必重算一次) 因此广受欢迎 更高阶的integrator也有很多，最有名的应该是Runge-Kutta methods 要一一介绍那就真是个巨大无比的坑 这个领域就叫做Simulating hamiltonian dynamics 似乎是数学物理一个最近还算有热度的课题 5-2.NVE系统有许多值得分析的数学性质，但绝大多数分子模拟并不采用NVE 因为真实世界基本上不存在绝热系统啊 常用的模拟是NVT(等温)和NPT(等压等温)模拟 统计力学里对应到canonical ensemble和grand canonical ensemble 具体的操作上 NVT是藉由导入随机作用力来完成--也就是Langevine dynamics， 也可以藉由rescale粒子速度来完成，这些方法叫thermostat NPT则是除了随机作用力之外还不断调整容纳系统的「箱子」的大小 来确保平均反弹力守恒，这叫barostat 关於如何纳入随机力和如何调整箱子大小又有各种方法 关於如何定义温度和压力也有各种方法 （是的，比如说，可以由位能而非动能来定义温度） 如何把各种thermostat和barostat写入integrator又有各种方法 要深究也是个做不完的大坑 常用的古典方法就讲到这里 下一篇（如果有）讲包含量子力学的动力学模拟 --

※ 发信站: 批踢踢实业坊(ptt.cc), 来自: 198.47.45.153 (加拿大) ※ 文章网址: https://webptt.com/cn.aspx?n=bbs/Physics/M.1644743796.A.3B8.html ※ 编辑: caseypie (198.47.45.153 加拿大), 02/13/2022 17:32:22