水分子的力場模型

地球上的水的存在範圍和狀態非常廣泛，從數千米、甚至數萬米的地下高溫高壓狀態，到地球兩極及其其它高寒地區的固體狀態，再到地球表面的常溫常壓狀態，以及大氣層中的氣體狀態；從海水的無機鹽溶液，到生物體內的生物分子溶液等。如此廣泛的存在範圍和狀態，使水錶現出異常複雜、多樣的性質。

在科學技術高度發達的當代社會，水的作用不但沒有減小，反而比以往任何時候更加重要。生活中，我們不但需要每天飲用水以補充身體中的水份，而且需要用水洗滌我們的身體、衣物和其它物品。農牧業生產中，我們需要水灌溉農田、餵養家畜。工業生產中，我們需要水作為溶劑、洗滌劑、冷卻劑等。目前，隨著人口數量的增加、人們生活水平的提高、工農業生產和第三產業的發展，人們的生活和生產用水需求量迅速增加，缺水已經成為許多城市和地區面臨的最嚴重的問題。同時，環境污染引起的對水資源的破壞、氣候變化引起的水的分布的不平衡，進一步加劇了水資源的緊張。但是，地球上的水並沒有因為我們的合理、或不合理的使用而有任何減少，我們也無法人工製造出一滴水。為了保護水、更加有效地利用水，我們需要從各個方面研究水、研究水的結構與性質。

1. 水分子力場的研究歷史

早在1800年，英國化學家Nicholson通過電解水發現，水由氫和氧兩種元素組成。1805年, Gay-Lussac和von Humboldt證明水由兩份氫元素原子和一份氧元素原子組成，可以表示成H₂O的形式。在Loschmidt和von Hofmann的水分子結構圖和球棍模型中，一個水分子被描述成由一個氧原子和兩個氫原子組成的V形結構。在V形結構中，氧原子通過化學鍵與兩個氫原子相連，但兩個氫原子之間沒有直接的化學鍵連線，H–O–H三個原子之間形成一個夾角。但是，由於缺乏足夠的結構測定手段，在提出水分子V形結構模型後的很長時間內，都無法用實驗方法直接驗證水分子的結構。後來，X-射線衍射技術的發明，為分子結構測定提供了強有力的手段。但由於氫原子質量小，周圍電子密度低，X-射線衍射技術並不能解決水分子的結構問題。因此，直到1932年，才首次通過波譜學方法證明水分子的V形結構。

1933年，即實驗直接驗證水分子的V形結構模型的次年，Bernal和Fowler提出了兩個水分子力場。他們還用所提出的水分子力場，估算水分子的多種性質，與實驗相比效果良好。目前，Bernal和Fowler的水分子模型，雖然已經沒有任何實際使用價值，但仍具有歷史意義。此後，Rowlinson在1949和1951年先後提出了三個水分子力場 。在Rowlinson提出水分子力場的近20年後，Baker進行首次水分子模擬時，選用了Rowlinson的一個力場模型。在分子模擬技術發展以前，我們難於評估不同水分子力場的質量和精確程度，研究水分子力場雖然具有理論意義，但實用價值很小。

1953年，Metropolis等首次提出了後來被稱為Metropolis採樣的算法，實現分子體系的Monte Carlo模擬，標誌計算機模擬技術開始用於研究分子體系的結構和性質。除Monte Carlo模擬外，MD模擬方法不久也被建立，用於模擬硬球體系的相變現象。隨著計算機運算速度的提高，分子模擬方法和程式的開發和完善，以及可以獲得計算機計算服務的科研人員的範圍的擴大，分子模擬方法不久就被用於研究分子體系。1969年，Baker和Watts首次實現了對水分子體系的Monte Carlo模擬研究。分子模擬方法用於研究水分子及其其它含水分子的體系，對水分子力場提出了更新更高的要求，促進了水分子力場的發展。1976年，Matsuoka、Clementi和Yoshimine利用量子化學計算方法，首次得到兩個水分子之間的相互作用勢能函式（MCY模型），標誌水分子間相互作用勢能函式不再是經驗的猜想 。1981年，Berendsen等開發了一個可以精確描述液體水的力場，即SPC模型（simple point charge model），取得巨大成功。目前，SPC模型仍是被最廣泛使用的水分子力場之一，並得到不斷的改進和發展 。1993年，Laasonen等首次利用密度泛函理論，實現對液體水的第一性原理MD模擬。第一性原理MD模擬方法的發展和推廣使用，為計算機模擬方法研究化學反應提供了廣闊的發展空間，是目前發展迅速的重要研究前沿 。

目前，公認的O–H鍵長的實驗值為0.95718 Å，H–O–H鍵角為104.523^o，水分子的三個振動模式的振動波數分別為3755.79，3656.65，1594.59 cm^-1。

2. 水分子力場的種類

與其它任何分子相比，水分子具有無與倫比的重要性。同時，水分子的經典力場的數量和種類也無與倫比，已有不下百種水分子力場。與研究任何其它分子的力場一樣，可以從不同層次對水分子進行近似，得到不同類型的力場。目前，水分子經典力場大致包括下面四個層次（或類型）：剛性力場（rigid force fields）、柔性力場（flexible force fields）、可極化力場（polarizable force fields）、可離解力場或反應力場（dissociable force fields or reactive force fields）。

剛性力場 剛性力點模型是經典力學對水分子的最初級近似，反映了水分子結構和性質的最本質屬性。

水分子由一個氧原子和兩個氫原子組成，最直接的水分子力場，是用三個分別與水分子的三個原子對應、位置相對固定的力點（即剛體）近似水分子，這就是所謂的三力點模型。在實際的三力點模型中，有的模型的三個力點的相對位置與氣態孤立水分子的三個原子的實際位置一一對應，也有的模型不與水分子的三個原子一一對應。事實上，不同狀態下的水分子的結構參數（即鍵長、鍵角等幾何參數的實驗數值）有所不同；因此，在不同水分子力場中選取不同的力點位置，不足為奇。由於三個力點的相對位置固定，組成一個剛性的水分子模型，只需要考慮水分子間的相互作用，不需要考慮分子內的相互作用，這也是各種剛性水分子模型的共同特點。三力點水分子模型的分子間相互作用力，由靜電相互作用和van der Waals相互作用兩個部分組成。其中，與氧原子對應的力點帶負電，分別與氫原子對應的兩個力點帶正電。靜電相互作用由分布在力點上的點電荷決定，通常根據庫倫定律計算。值得注意的是，在ST2和BNS等力場中，靜電能並不完全根據庫倫定律計算。van der Waals相互作用利用球形對稱的Lennard-Jones勢函式近似，作用中心位於氧原子上。因此，水分子體系的總勢能可以表示為，

其中，前面一個雙重求和遍及所有電荷點，後面一個雙重求和遍及所以van der Waals力點，但不包括非van der Waals力點。因此，只要確定水分子模型的力點相對位置和勢參數，就可以計算水分子體系的總勢能。

如果在三力點模型的基礎上增加更多的力點，就可以得到具有更多力點的水分子模型。例如，在四力點模型中，為了模擬負電荷中心不與氧原子重疊的情況，與氧原子對應的負電荷中心被沿H–O–H鍵角的平分線移到更靠近氫原子的位置M，而van der Waals力點的位置仍然不變。在五力點模型中，負電荷被分解為兩部分，分別與氧原子的兩對孤對電子對應，van der Waals力點位置不變。除最常用的三力點、四力點、五力點模型外，還有更多力點的水分子模型，甚至還有單力點、雙力點的水分子模型。更詳細的內容，讀者可以參考Robinson的《Water in biology, chemistry, and physics: experimental overviews and computational methodologies》一書的第五章，或其它有關文獻。

柔性力場剛性水分子模型是對水分子的一級近似，不管在模型中引入多少個力點，都不能模擬振動光譜等與水分子內部運動相關的性質。這就需要引入分子內自由度，也就是柔性水分子模型。開發柔性水分子模型最常用的方法，是在剛性水分子模型的基礎上引入O–H鍵的伸縮自由度，H–O–H鍵角彎曲自由度等。但是，在引入內部自由度後，原剛性力場中的幾何參數需要略作調整，以得到更好的模擬結果。常用的柔性力場有TIP4P/FQ模型、柔性SPC模型等。

可極化力場在早期的分子力場中，電荷點和其它力點一樣，被賦予固定的位置、電荷數或偶極矩。這樣的模型有一個很大的缺點，不允許分子中電荷的重新分布，不能描述分子的極化現象。這種不考慮分子中電荷重新分布的不可極化力場，雖然取得了巨大的成功，但無法描述水的許多實驗現象。事實上，當一個水分子與另一個水分子相互靠近時，將發生極化效應，引起水分子中電荷的重新分布。在分子力場中計入極化效應的分子力場，被稱為可極化分子力場。同時，還可以根據極化過程中產生誘導電荷、偶極矩、多極矩等不同的極化效應，把可極化分子力場分為多種類型。可極化分子力場較不可極化分子力場複雜，計算量也大大增加。可極化力場的具體效果不但與力場的精度有關，也與實際研究體系的本質有關。最早的可極化水分子力場是PE力場。

我們知道，在分子模擬中廣泛採用的兩體勢或對勢，實際上是一種有效的兩體勢；不但包括了兩個指定分子或原子之間的相互作用，還包括了平均化了的周圍分子（處於平均位置）對這兩個指定分子的相互作用，這就是所謂的多體相互作用。如果以量子化學方法精確計算兩個水分子之間的相互作用勢函式，並以此模擬大量水分子的性質，模擬結果將顯著偏離實驗結果。相反，利用經仔細最佳化、包括多體相互作用的有效兩體勢進行類似的模擬，模擬結果較利用精確兩體勢得到的模擬結果更接近實驗值。這種精確的兩體勢與有效兩體勢之間的差別，就是分子的極化效應引起的結果。本質上，多體相互作用是一種極化相互作用：兩個指定分子之間的相互作用與周圍分子的存在相關；當周圍分子靠近兩個指定分子時，兩個指定分子之間的相互作用因極化效應而加強，顯著偏離周圍沒有其它分子存在時的相互作用。由於極化效應的結果總是增強分子之間的純兩體相互作用，因此，有效兩體勢不可避免地大於純兩體勢。這也是經驗水分子力場中，模型的偶極矩往往顯著地大於氣體水分子的偶極矩，模擬得到的第二維里係數顯著大於實驗得到的第二維里係數的原因。

據分析，有效兩體勢的氫鍵勢阱應比實際勢阱深25%，偶極矩比氣相偶極矩大25%，這樣才能有效包括多體勢的相互作用。

可離解力場或反應力場 我們知道，與AIMD模擬相比，經典MD模擬方法最嚴重的缺點是無法模擬化學反應，也就是化學鍵的斷裂與生成。事實上，這不是經典MD模擬的缺點，而是經典分子力場的缺點。如果在經典MD模擬中，採用可離解力場或反應力場；那么，就可以用經典MD模擬方法研究化學反應。

反應力場是經典力場的一種很有前途的擴展，成功實現了在經典MD模擬的框架內模擬化學鍵的生成和斷裂，卻比AIMD模擬方法具有更低的計算量、更高的模擬效率。

其它力場水分子是被最廣泛研究、模擬的分子。利用前述的各種水分子經典力場，可以模擬水團簇、液體水、水溶液、固體水等的結構和性質。這些模型的一個顯著特點是水分子在模型中顯式出現，被統稱為水的顯式模型（explicit model）。除此之外，MD模擬中還用到水分子的粗粒度模型（coarse-grained model）、隱式模型（implicit model）或連續溶劑模型（continuum solvent model）等，用於模擬作為溶劑的水分子對大分子、生物分子等的影響，達到降低計算量，提高模擬效率的目的。

3. 發展水分子力場的方法

現代分子力場的開發，依賴於一種稱為“循環擬合”（recursive fitting）的方法，以擬合力場參數與實驗數據。具體方法是：首先確定力場參數的初值；其次，利用力場參數的初值對水分子體系進行Monte Carlo模擬或MD模型，得到水分子的性質；第三，把模擬結果與有關實驗結果進行比較，修正力場參數，得到改進了的力場參數；第四，根據新的力場參數再次進行相關模擬，並進一步修正力場參數；如此反覆，直到模擬結果與實驗結果之間的差值小於預先設定的數值，得到收斂結果。目前，用於循環擬合法確定力場參數的常用實驗數據包括晶體結構數據（如劍橋結構資料庫，Cambridge Structural Database）、液體的徑向分布函式、分子的偶極矩、極化率、遷移性質、蒸發熱等。最近，也有人利用水的相圖數據，特別是密度反常數據擬合水分子力場參數，取得理想效果。

以上這種利用實驗數據擬合力場參數的方法在數學上也稱為逆問題方法（inverse problem）。用逆問題方法研究勢函式，與勢函式的形式無關，也與勢函式種類無關。任何包括分子內成鍵相互作用、分子間和（或）分子內氫鍵相互作用、van der Waals相互作用、靜電相互作用等的體系，都可以實現。

除了上述經驗方法外，分子力場參數還可以通過量子化學方法計算得到。具體是，先用量子化學方法直接計算水分子的二聚體、三聚體、或多聚體等的勢能面；然後，用適當的勢函式擬合勢能面，得到水分子力場的勢參數。以此方法得到的力場有MCY ，MCHO，NCC 等。值得注意的是，由於量子化學計算方法中沒有計入極化效應對有效兩體勢的貢獻，由此得到的勢函式雖然可以精確地描述少數水分子體系的結構和性質，但用於模擬由大量水分子組成的巨觀水的性質時的效果較經驗勢函式不理想。比如，利用MCY模型模擬得到的水的徑向分布函式，第一峰的位置是合理的，但第二峰的位置太近，總體與水的高溫行為（具有較弱的分子間相互作用）類似。

4. 常用水分子力場

早期的水分子力場 Bernal 和 Fowler 提出的BF模型由三個點電荷和一個van der Waals相互作用中心共四個力點構成。由於該模型詳細指定了各個力點的幾何位置、所帶電荷量、Lennard-Jones勢參數，成為一個完全確定的力學模型，是第一個經典力學意義上的水分子模型（圖a）。由於BF模型的簡潔性、以及對水分子性質的預報能力（雖然定量效果不理想），該模型成為許多後繼模型的原型和範本。比如，TIPS2模型中的力點的數量和位置與BF模型完全一致，只是把點電荷電量調整為0.535。

(a)Bernal-Fowler 水分子模型

(b) Rowlinson 水分子模型

考慮到氧原子有兩對孤立電子，Rowlinson在BF模型的基礎上把負電荷中心分解為兩部分，分別位於氧原子的上下0.2325 Å的位置，並將單個正電荷的電荷量調低至0.3278，同時稍微調整其它力場參數（圖b）。這樣，Rowlinson模型有五個力點，可以看作是五力點模型的原型。Rowlinson模型的最大缺點是Lennard-Jones勢參數中碰撞因子s 太小，僅為2.8 Å，大致與冰中氫鍵的最近鄰O–O距離對應。第一個被廣泛套用於MD模擬的五力點模型是ST2模型。ST2的原型是BNS模型。與BNS模型一樣，ST2模型中靜電相互作用不純粹由庫倫定律計算，而是在庫倫定律的基礎上經開關函式調整後計算得到，

簡單水分子模型 雖然早期的水分子力場在精確程度等方面存在不足，但它們的形式卻在以後發展的力場中保留了下來，成為一種規範。不同的是，不同水分子模型的電荷點的數量、位置、電量不同，van der Waals力點的數量、位置、力場參數不同。一般地，由於優先考慮的性質不同，不同的水分子模型的電荷點的數量和位置既可以與水分子中氫、氧原子對應，也可以不對應。通過調整電荷點的數量、位置、電量，可以調整水分子模型的偶極矩，並與水分子的實驗偶極矩進行比較。由於水分子的氣相和液相偶極矩分別為1.855 和2.22 D，只通過調整電荷點數量、位置、電量，一般無法使模型的氣相偶極矩和液相偶極矩與水分子的相應實驗值的結果同時符合。水分子模型的第二要素是van der Waals力點的數量、位置、勢參數。由於氫原子的van der Waals半徑遠小於氧原子的van der Waals半徑；H–H、H–O的van der Waals相互作用範圍被O–O間相互作用範圍所掩蓋，大多水分子模型中只考慮O–O間的van der Waals相互作用。受van der Waals相互作用影響最大的性質是水的密度和蒸發熱，通過調整O–O間van der Waals相互作用的Lennard-Jones勢函式的碰撞因子s ，可以調節水分子的模擬密度；通過調節Lennard-Jones勢函式的能量因子e ，可以很容易地調整水分子的模擬蒸發熱。

目前，大多簡單水分子模型由三到六個力點組成。如TIP3P模型（transferable interatomic potential with three points model）、SPC模型（simple point charge models）等由三個力點組成，分別與水分子的三個原子對應。其中，一個負電荷點與氧原子對應，兩個正電荷點與氫原子對應。不同的是，TIP3P模型中正電荷點與氫原子的位置重合，但SPC模型中正電荷點的位置與水分子中的氫原子並不重合。我們知道水分子的O–H實驗鍵長為0.9572 Å，H–O–H實驗鍵角為104.52^o；但在SPC模型中正負電荷點之間的距離為1.0 Å，略大於O–H鍵長的實驗值，三個電荷點之間的夾角為109.47^o，與sp雜化軌道一致，但與H–O–H實驗鍵角不一致。

(c) 水分子模型的力點及其分布

除三力點模型外，常用的還有四力點模型、五力點模型、六力點模型等。簡單水分子模型的力點分布情況如圖c所示。

雖然隨著力點數目的增加，力場參數也相應增加，力場也更加複雜；但力場的精確程度似乎與力點數量的關係不大。比如，具有四面體結構的五力點水分子模型，模擬得到的水的結構過分有序，不如平面結構的四力點水分子模型。此外，SPC、TIP3P、TIP5P、SPC/E等力場，不能描述液體水和凍的相圖，但TIP4P模型可以產生定性正確的相圖。對TIP4P模型進行重新最佳化參數後得到的TIP4P/2005模型可以產生定量正確的相圖，以及用於AIMD模擬的TIP4PQ/2005模型 。

研究水分子力場的目的和意義，不僅僅是為了模擬水的結構和性質。相反，研究水分子的力場的主要目的和意義基於這樣一種理念，一個可以精確預報水的相圖等各種性質的水分子經驗模型，肯定可以更好地預報水的結構、特別是水溶液中各種複雜生物分子等與水結合時的結構。這些複雜生物分子等在水溶液中的結構，在科學研究和實際套用中具有重要的意義，也是現代實驗技術難於表征的重要性質。

參考文獻