背景：力場是經典模擬計算的核心，因為它體現了結構中每種類型的原子如何與其周圍的原子相互作用。對于體系中的每個原子，分配一個力場類型，描述原子的局部環境。力場包含描述各種屬性的信息，例如平衡鍵長度和力場類型對之間的靜電相互作用。通常，力場是通用的，以覆蓋許多體系，或針對特定問題進行參數化。一般的力場，如Dreiding，在不同環境中有許多原子的參數。但是，對于特定情況，它永遠不會像僅針對該情況參數化的力場那樣表現良好。因此，能夠編輯力場以用于尚未對其進行參數化的體系會很有用。

這可能在許多情況下發生：您可能希望參數化力場以模擬沸石框架中分子的吸附等溫線，您可能具有要用于提高多晶型預測計算準確性的實驗晶體數據，或者您可能希望更改力場以更好地描述聚合物的密度。在所有這些情況下，您都需要編輯基礎力場并觀察這些更改的效果。

編輯力場是一項復雜的任務，因為許多交互都是顯式或隱式耦合的。例如，旋轉扭轉角的能量曲線不僅取決于扭轉項，還取決于非鍵能。因此，獲取正確的參數可能需要相當長的時間，并且涉及許多不同的計算。最好的方法是進行小的體系性更改，并查看對體系的影響。與所有形式的參數化一樣，需要在獲得可靠結果和過度擬合消耗資源之間權衡。

(資料圖片僅供參考)

簡介：在本教程中，您將修改小分子的單個扭轉角，以改善優化模型和晶體結構之間的匹配。最初，您將使用Dreiding力場優化爆炸材料三氨基三硝基苯(TATB)的晶體結構。然后，您將修改Dreiding力場以改善晶體結構和力場優化模型結構之間的匹配。您將結合使用Conformers模塊來探測扭轉能關系，DMol3來提供從頭算數據，Forcite來優化晶體結構。

目的：說明了在Dreiding中編輯力場參數的工作流程，以改進計算晶體結構數據與實驗晶體結構數據之間的擬合。

本教程重要節點：

優化晶體結構-檢查扭轉角度-使用非鍵項-修正范德華項

1. 優化晶體結構

打開Import Document對話框。導航至ExamplesDocuments3D Model文件夾，雙擊TATNBZ.xsd文件。將結構重命名為TATNBZ_crystal.xsd。

在將修改后的力場應用于晶體結構以驗證這些更改之前，將對孤立的分子執行所有力場編輯工作。應該復制一個孤立的分子。

在Project Explorer中，右鍵單擊project目錄，從彈出的快捷菜單中選擇New | 3D Atomistic。

在TATNBZ_crystal.xsd，單擊其中一個分子中的一個原子，右鍵并選擇Select Fragment。按下CTRL + C鍵。在新建的原子文檔中按下CTRL + V鍵。在Project Explorer中，將新的原子文檔重命名為TATNBZ_molecule.xsd。

如果檢查晶體結構，會發現硝基在苯環所在的平面上。現在，將使用Forcite模塊利用標準Dreiding力場優化結構。

確保TATNBZ_crystal.xsd文件為當前活動文檔。單擊Modules工具條上的Forcite按鈕，在下拉列表中選擇Calculation，打開Forcite Calculation對話框。

將Task修改為Geometry Optimization。在Energy選項卡中，將Forcefield設置為Dreiding，并將Charges更改為Charge using Gasteiger。

可以使用更復雜的電荷計算方法，例如使用DMol3應用ESP電荷，但Gasteiger電荷將適用于自定義力場的開發。

單擊Run按鈕。計算任務完成后，打開優化后的TATNBZ_crystal.xsd文檔。

應該看到硝基不再和苯環在處于同一平面上。優化后的結構中硝基的離面性質會影響這些分子在晶體中聚集的方式，因此獲得平面硝基非常重要。

硝基和苯環之間有一個單一的扭矩，這將控制基團的平面度。修改控制此扭矩的參數將是本教程的重點內容。

2. 檢查扭轉角度

現在已經確定了要更改的扭轉角，可以開始更詳細地查看該扭轉角的能量分布。通過計算扭矩旋轉360°時的單點能，Conformers計算將提供該能量分布。

右鍵單擊TATNBZ_molecule.xsd，從彈出的快捷菜單中選擇Label，打開Label對話框。確保Object type已設置為Atom，從Properties列表中選擇Name。單擊Apply按鈕。

使用Measure/Change工具在硝基上定義原子序列為O5-N5-C5-C4的扭矩。

在Label對話框中，單擊Remove按鈕，關閉對話框。

將使用Conformers模塊系統性地旋轉扭矩，并在每點計算能量。

單擊Modules工具條上的Conformers按鈕，在下拉列表中選擇Calculation，打開Conformers Calculation對話框。單擊Torsions…按鈕，打開Conformers Torsions對話框。

將# STEPS修改為72，關閉對話框。

在Energy選項卡中，將Forcefield設置為Dreiding，并將Charges更改為Charge using Gasteiger。單擊Run按鈕，關閉對話框。

當計算結束后，即可以查看扭矩-能量圖，從而觀察能量分布。

在數據表中選擇B列和C列，單擊Quick Plot按鈕。

將圖表文件重命名為Dreiding_Original.xcd。

Dreiding力場計算硝基-苯扭轉角的扭矩-能量圖

從扭矩-能量圖中可以看出，在0°和180°處存在兩個絕對值較大的能量最大值，這為旋轉提供了約14 kcal mol-1的勢壘。通過檢查晶體結構，可以預期在0°和180°處存在能量最小值，而不是扭矩-能量圖所示的最大值。

也可以將這個勢能面與DMol3預測的勢能面進行比較。由于已經有一個包含構象的數據表，可以使用DMol3模塊來計算每個構象的能量。

注意：由于此計算需要精確的能量，將使用較高的精度設置，此計算可能需要一個多小時。ExamplesStudyTables文件夾中提供了包含計算結果的數據表。如果希望使用此文件，請跳過接下來的兩個步驟。

在數據表中，選擇A列。單擊Models按鈕，

打開Models對話框。定位至DMol3 Molecular Energy模型，雙擊并編輯。

對于該計算，需要精確能量。應使用高精度泛函和收斂標準。

在Model Editor – DMol3 Molecular對話框的Inputs選項卡中，將Functional修改為BLYP，Quality level修改為Fine。單擊Save按鈕，關閉對話框。在Models對話框中單擊Run按鈕。

高精度的精度等級不僅更改SCF收斂標準，同時也改變其他設置，如基組。該計算需要一定時間才能完成，建議遠程運行或隔夜運行。必須等待計算完成，才能進入下一步，也可導入示例數據表。

打開數據表，或從ExamplesStudyTablesTATNBZ_conformers.std導入示例數據表。

將觀察到Total Energy (DMol3 Molecular)的單位為Hartree，因此需將其轉換為kcal mol-1 (1 Hartree = 627.51 kcal mol-1)。假設總能量列為列D，則需要定義一個函數D×627.51。一旦完成單位轉換，應該像繪制原始構型圖一樣繪制能量圖。

單擊Define Function按鈕，

打開Define Function對話框。在Expression文本區輸入D*627.51。輸入DMol3 Total Energy (kcal/mol)作為Name，單擊OK按鈕。選擇包含扭轉的列和新創建的列，單擊Quick Plot按鈕。

注意：能量的絕對值將非常高，因為它們沒有使用DMol3進行優化。然而，這里只關心相對勢壘高度，而不關注能量絕對值。

DMol3單點能的二面角-能量圖

從圖表中可以看出，勢壘高度約為14 kcal mol-1。在本教程的后面部分，將嘗試從力場計算中獲得相同的近似勢壘高度。

現在已計算了由Conformers和DMol3預測的能量分布，可以將其與Dreiding中的扭矩項進行比較。Forcefield Manager用于創建標準力場的可編輯副本。

單擊Modules工具條上的Forcite按鈕，在下拉列表中選擇Forcefield Manager。在Standard Forcefields部分，選擇Dreiding。單擊>>按鈕。關閉對話框。

將在工程中創建Dreiding力場的副本，并打開Dreiding.off力場供編輯。力場文檔由四個選項卡組成：

Summary – 包含描述力場的摘要文本。

Types – 包含力場類型及其關聯的屬性。

Interactions – 包含主要的相互作用，如價鍵項和非鍵項。

Equivalences – 包含等效相互作用的定義。

編輯力場之前，可以更改力場文檔的名稱。

在Project Explorer中，將力場文檔的名稱更改為Dreiding_new.off。

現在已準備好查看力場文件。

在Dreiding_new.off文件中，選擇Types選項卡。

其中包含Dreiding中可用的力場類型，以及每個力場類型的類型、元素和范德華參數的說明。可以使用Forcefield Type Properties對話框上的復選框顯示其他屬性，例如雜化、電荷和氫鍵。

有許多力場類型描述每個元素的不同局部環境，有些元素具有多個力場類型。滾動瀏覽所有這些內容會很枯燥，因此可以按力場類型過濾顯示的內容。過濾器框在力場文檔中顯示為黃色。

在黃色的Type Filter文本框中，輸入C_*并單擊ENTER鍵。

現在，Types對話框中僅顯示與碳原子有關的原子類型。

將Type Filter更改回*。

也可以按3D模型文檔中存在的力場類型進行過濾。

勾選Filter by selection in，在下拉列表中選擇TATNBZ_molecule.xsd。

將看到一個警告對話框。這是因為分子中的原子沒有力場類型。可以使用Forcite Calculation對話框分配力場類型。

打開Forcite Calculation對話框。在Energy選項卡中，單擊與Forcefield關聯的More…按鈕，打開Forcite Preparation Options對話框。取消選擇Forcefield types的Calculate automatically復選框，將TATNBZ_molecule.xsd更改為當前文檔。單擊Calculate按鈕。再次勾選Calculate automatically。關閉所有Forcite對話框。

這計算了TATNBZ_molecule中原子的力場類型。需要刷新力場文檔中的力場類型。

將Dreiding_new.off更改為當前文檔。重新選擇TATNBZ_molecule.xsd。

此時可以觀察到代表文檔中的原子的四個力場類型，和一個虛擬力場類型。

Dreiding_new力場文件，顯示了力場類型信息

本例中所關注的是芳香族碳C_R和芳香族氮N_R之間的扭轉角。

打開Interactions選項卡。將Show interaction更改為Torsion。

扭矩項有六種不同的組合。可以更改Functional Form過濾器以查看不同的選項。

單擊Functional Form過濾器并選擇Dihedral。

顯示描述扭矩的二面體函數形式中使用的參數值。在線幫助中記錄了所有參數。

選擇顯示的行之一，然后按下F1鍵。

將顯示幫助，顯示不同的函數形式。二面體扭轉的函數形式為：

可以使用之前生成的Conformers數據表繪制純扭矩項的變化，并比較Conformers和DMol3得出的整個能量表達式的函數形式。對于力場相互作用，將上述方程中的Bj替換為25，d替換為1，nj替換為2。將使用數據表中B列的角度值，轉換為弧度。

使得TATNBZ_conformers.std為當前文檔，打開Define Function對話框。在Expression文本區，輸入(25*(1-1*cos(2*B*0.0174532925)))/2。將Name設置為Dihedral from Dreiding，單擊OK按鈕。選擇B和J列，單擊Quick Plot按鈕。

生成的圖表應與下面的圖表相匹配，顯示與DMol3圖相對應的最小值和最大值，但能壘高于預期。

使用X C_R N_R X扭矩力場編輯器得到的純二面角分項繪制二面角-能量圖

從上面的圖中，可以看出二面角扭轉看起來是合理的。這表明問題一定出在其他地方。在下一節中，將研究非鍵項，以了解它們如何影響扭轉能量。

3. 使用非鍵項

從前面的計算中，觀察到純二面角的能量貢獻與觀察到的扭矩旋轉時的能量變化顯著不同。可以使用力場編輯器關閉體系中不同類型相互作用的貢獻，以查看它們的效果。因為只改變了一個末端扭轉角，除了扭矩本身的能量，改變的貢獻是非鍵能。因此，將關閉非鍵能，以查看它們對扭矩-能量圖的影響。

由于硝基的氧原子很可能與相鄰的胺形成強氫鍵，因此首先禁用氫鍵項是有意義的。

打開Dreiding_new.off文檔。在Interactions選項卡中，單擊More…按鈕。

打開Forcefield Preferences對話框，在其中可以設置力場的的一般首選項。

在相互作用列表中取消勾選Hydrogen Bond復選框。

提示：在將力場文檔與Conformers或Forcite等模塊一起使用之前，無需保存該文檔。如果進行了大量更改，建議將其保存。

現在你需要用新的力場運行另一個Conformers計算，以觀察關閉氫鍵相互作用的效果。

在工程根目錄中打開TATNBZ_molecule.xsd文件。打開Conformers Calculation對話框，選擇Energy選項卡。從Forcefield下拉列表中選擇Browse…，打開Choose Forcefield對話框。選擇Dreiding_new.off，并確保Charges設置為Charge using Gasteiger。

在Forcefield文本框中將顯示Dreiding_new。反斜杠表示已從工程中選擇了力場，它不是標準的力場。

注意：此力場將在運行計算時傳輸到服務器，并將返回到新創建的計算結果文件夾中。

在Conformers Calculation對話框中，單擊Run按鈕。

當計算任務結束后，打開TATNBZ_molecule.std文檔，對B和C列繪制曲線圖。將其與Dreiding_Original.xcd進行比較。

應注意到禁用氫鍵相互作用對扭矩-能量圖的影響微乎其微。接下來，試著金庸靜電相互作用。

將Dreiding_new.off文檔改為當前文檔。在Forcefield Preferences對話框的相互作用列表中取消勾選Electrostatic復選框。將工程根目錄中的TATNBZ_molecule.xsd更改為當前文檔。打開Conformers Calculation對話框，選擇Energy選項卡。

修改過的Dreiding版本仍然為所選力場，因此可以僅運行計算任務。

運行Run Conformers計算。當計算任務結束后，繪制扭矩和能量關系的曲線圖，將其與Dreiding_Original.xcd進行比較。

可見，與打開靜電相互作用相比，能量較低的極大值在高度上有所降低。另外，總能量增加了幾kcal mol-1。

最后，關閉范德華相互作用以查看效果。

將Dreiding_new.off文檔改為當前文檔。在Forcefield Preferences對話框的相互作用列表中取消勾選van der Waals復選框。對TATNBZ_molecule.xsd文檔運行Conformers計算。在輸出數據表中，繪制扭矩和能量關系的曲線圖，將其與Dreiding_Original.xcd進行比較。

當關閉范德華相互作用項時，應在扭矩-能量圖上看到重大的變化。能量分布應該與通過二面角相互作用生成的曲線圖非常相似。

這表明范德華項在這個扭轉角的能量分布中起著主要作用。修改范德華項將是下一節的重點。

將Dreiding_new.off文檔改為當前文檔。在相互作用列表中勾選van der Waals、Hydrogen Bond和Electrostatic復選框，關閉Forcefield Preferences對話框。

去除范德華勢后產生的變化與硝基上的氧原子和胺基上的氫原子之間的近距離有關。通過修改范德華項，使其具有較短的平衡距離，從而在較短的距離內減少排斥力，應該看到，Conformers得到的圖開始與DMol3預測的更接近。

在工程根目錄中將TATNBZ_molecular.xsd設置為當前文檔。使用Measure/Change工具在氧原子和相鄰氫原子之間添加Distance測量。

距離應約為1.8 ?。在下一節中修改范德華項時，將使用此信息。

4. 修正范德華項

Dreiding中的范德華項是根據每種力場類型的參數組合自動生成的。由于參數存儲在每個力場類型中，因此它們將顯示在Types選項卡上的力場文檔中。還可以為原子對添加自定義范德華項，這些項將被用來代替自動計算的項。

在Dreiding_new.off文件中，選擇Types選項卡。單擊van der Waals Form的黃色過濾器框，然后選擇LJ_12_6。

注意：可以選擇顯示函數形式，也可以Ignore相互作用。如果選擇Ignore，則不會在能量表達式中計算相互作用。為van der Waals力場類型選擇Ignore，將忽略包含該類型的所有范德華相互作用。盡管這在本例中看起來過于激進，但如果想關閉單獨的相互作用，ignore可能很有用。

可以設置函數形式，也可以選擇忽略范德華參數。Dreiding使用簡單的Lennard Jones 12-6形式計算范德華參數，其中D0為阱深，R0為平衡距離。

記下H___A和O_2的D0和R0值。

可見與O_2相比，H___A的阱深非常淺。參數組合通常使用算術平均值計算。然而，由于氫的阱深非常小，在這種情況下，阱深參數的幾何組合可以更好地表示范德華相互作用。

幾何平均值計算為阱深乘積的平方根，對于H___A和O_2等于0.0031。

打開Interactions選項卡。將Show interaction更改為van der Waals。將Functional Form更改為LJ_12_6。

注意到，沒有為Dreiding定義特定的范德華項。要引入特定的范德華項，必須首先指定力場類型序列。

在Fi選框中單擊第一個空行，然后從下拉列表中選擇H___A。單擊Fj選框并從下拉列表中選擇O_2。在Functional Form選框中單擊并選擇LJ_12_6。將D0的值設置為0.0031。

對于該測試，應該使用算術平均值計算R0。該值為3.29。

將R0的值設置為3.29。在工程根目錄中的TATNBZ_molecule.xsd文件中運行Conformers計算。生成扭矩-能量圖。

可見極小值的高度已經改變了約0.5 kcal mol-1，但對曲線圖沒有太大的總體影響。改變函數形式的另一種方法是減小平衡距離的值。之前，計算了分子中氧原子到氫原子的距離。這約為1.8 ?，而R0為3.29 ?。這表明范德華勢的排斥項對扭矩勢有重要貢獻。要糾正此問題，應減小R0的值。雖然很容易將平衡值設置為測量距離的平衡值，但這可能會過度擬合相互作用，從而導致整體結構的變化。事實上，選擇正確的數字可能需要一些猜測，但最好先較高值開始逐漸降低，運行多次Conformers計算。在本教程中，使用了2.7 ?的R0值。

在Dreiding_new.off中，將R0的值更改為2.7。在TATNBZ_molecule.xsd中運行Conformers計算。生成扭矩-能量圖。

力場編輯器中調整X C_R N_R X扭矩的范德華參數后的二面角-能量圖

曲線圖顯示，改變平衡位置對扭轉角的能量分布有顯著影響，在0度和180度時能量最大值顯著降低。然而，總體能量勢壘現在明顯低于先前DMol3計算中預測的能量勢壘。為了補償這一點，應該增加扭轉能量勢壘。

在Dreiding_new.off中，將Show Interaction更改為Torsion，Functional Form過濾器更改為Dihedral。將B的X N_R C_R X值更改為50。在TATNBZ_molecule.xsd中，運行Conformers計算。產生扭矩-能量圖。

可見，新的扭矩-能量圖現在與從頭算計算預測的結果相似，扭轉勢壘的大小和位置也相似。在0度和180度的兩側仍有兩個低能極小值，這將導致結構變形。然而，改變H___A O_2范德華相互作用的值將不再有太大的影響。

這表明小的排斥項一定是由另一個原子與硝基上的氧相互作用引起的。

將工程根目錄中的TATNBZ_molecule.xsd文件打開為當前文檔。在氧原子和最近鄰胺基氮原子上添加一個Distance測量。

距離應約為2.5 ?。與氫-氧相互作用一樣，這小于力場中的平衡距離。該距離意味著范德瓦爾斯斥力將很大，這可能導致勢能面上的絕對值較小的最大值。

將Dreiding_new.off打開為當前文檔。選擇Types選項卡，記下N_R和O_2的D0和R0的值。

與氫和氧范德華項一樣，將使用幾何平均值計算新的D0項，并將R0設置為2.7 ?。

打開Interactions選項卡。僅顯示van der Waals相互作用。單擊黃色過濾器選框Functional Form，選擇LJ_12_6。單擊空白的Fi選框，選擇N_R。單擊空白的Fj選框，選擇O_2。在Functional Form選框中單擊并選擇LJ_12_6。將D0的值設置為0.086，R0的值設置為2.7。在TATNBZ_molecule.xsd中，運行Conformers計算。產生扭矩-能量圖。

不應存在額外的最小值，圖表的形狀應與之前生成的DMol3圖相似，最小值略為平坦。然而，當改變氧和氮的范德華參數的函數形式時，二面角的能量勢壘現在太高了。根據DMol3計算，希望其約為14 kcal mol-1。

將Dreiding_new.off文檔改為當前文檔。在Interactions選項卡中，選擇Torsion相互作用。在Functional Form中，選擇Dihedral。將B的X N_R C_R X值更改為35。在TATNBZ_molecule.xsd中，運行Conformers計算。產生扭矩-能量圖。

這將得到約為14 kcal mol-1的勢壘。可以繼續調整參數以進一步改善力場，但現在的力場應該足以測試晶體結構。

打開原始的TATNBZ_crystal.xsd文件。打開Forcite Calculation對話框。選擇Energy選項卡，選擇Dreiding_new.off作為力場。單擊Run按鈕。

當計算完成后，可見硝基幾乎在苯環平面上排列，就像它們在晶體結構中一樣。

從這里，可以使用新的力場運行多晶型預測Polymorph Prediction計算。