這篇繼續講古典物理、熱力學和動力學模擬方法 首先我推薦想鑽研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/m.aspx?n=bbs/Physics/M.1644743796.A.3B8.html ※ 編輯: caseypie (198.47.45.153 加拿大), 02/13/2022 17:32:22