分子動(dòng)力學(xué)模擬I

Gromacs中文教程淮海一粟分子動(dòng)力學(xué)（MD丿模擬分為三步：首先，要準(zhǔn)備好模擬系統(tǒng)；然后，對(duì)準(zhǔn)備好的系統(tǒng)進(jìn)行模擬；最后，對(duì)模擬結(jié)果進(jìn)行分析。雖然第二步是最耗費(fèi)計(jì)算資源的，有時(shí)候需要計(jì)算幾個(gè)月，但是最耗費(fèi)體力的步驟在于模擬系統(tǒng)準(zhǔn)備和結(jié)果分析。本教程涉及模擬系統(tǒng)準(zhǔn)備、模擬和結(jié)果分析。一、數(shù)據(jù)格式處理準(zhǔn)備好模擬系統(tǒng)是MD最重要的步驟之一。MD模擬原子尺度的動(dòng)力學(xué)過程，可用于理解實(shí)驗(yàn)現(xiàn)象、驗(yàn)證理論假說，或者為一個(gè)待驗(yàn)證的新假說提供基礎(chǔ)。然而，對(duì)于上述各種情形，都需要根據(jù)實(shí)際情況對(duì)模擬過程進(jìn)行設(shè)計(jì)；這意味著模擬的時(shí)候必須十分小心。丟失的殘基、原子和非標(biāo)準(zhǔn)基團(tuán)本教程模擬的是蛋白質(zhì)。首先需要找到蛋白質(zhì)序列并選擇其起始結(jié)構(gòu)，見前述；然后就要檢查這個(gè)結(jié)構(gòu)是否包含所有的殘基和原子，這些殘基和原子有時(shí)候也是模擬所必需的。本教程假定不存在缺失，故略去。另一個(gè)需要注意的問題是結(jié)構(gòu)文件中可能包含非標(biāo)準(zhǔn)殘基，被修飾過的殘基或者配體，這些基團(tuán)還沒有力場(chǎng)參數(shù)。如果有這些基團(tuán)，要么被除去，要么就需要補(bǔ)充力場(chǎng)參數(shù)，這牽涉到MD的高級(jí)技巧。本教程假定所有的蛋白質(zhì)不含這類殘基。結(jié)構(gòu)質(zhì)量對(duì)結(jié)構(gòu)文件進(jìn)行檢查以了解結(jié)構(gòu)文件的質(zhì)量是一個(gè)很好的練習(xí)。例如，晶體結(jié)構(gòu)解析過程中，對(duì)于谷氨酰胺和天冬酰胺有可能產(chǎn)生不正確的構(gòu)象；對(duì)于組氨酸的質(zhì)子化狀態(tài)和側(cè)鏈構(gòu)象的解析也可能有問題。為了得到正確的結(jié)構(gòu)，可以利用一些程序和服務(wù)器（如WHATIF）。本教程假定所用的結(jié)構(gòu)沒有問題，我們只進(jìn)行數(shù)據(jù)格式處理。二、結(jié)構(gòu)轉(zhuǎn)換和拓?fù)浠粋€(gè)分子可以由各個(gè)原子的坐標(biāo)、鍵接情況與非鍵相互作用來確定。由于.pdb結(jié)構(gòu)文件只含有原子坐標(biāo)，我們首先必須建立拓?fù)湮募?，該文件描述了原子類型、電荷、成鍵情況等信息。拓?fù)湮募?duì)應(yīng)著一種力場(chǎng)，選擇何種力場(chǎng)對(duì)于拓?fù)湮募慕⑹且粋€(gè)值得仔細(xì)考慮的問題。這里我們用的是GROMOS9653a6連接原子力場(chǎng)，該力場(chǎng)對(duì)于氨基酸側(cè)鏈的自由能預(yù)測(cè)較好，并且與NMR試驗(yàn)結(jié)果較吻合。

拓?fù)鋽?shù)據(jù)要與結(jié)構(gòu)吻合，這點(diǎn)很重要；這意味著結(jié)構(gòu)數(shù)據(jù)也需要在所用的力場(chǎng)中進(jìn)行轉(zhuǎn)換。為了轉(zhuǎn)換結(jié)構(gòu)并建立拓?fù)鋽?shù)據(jù)，可以用pdb2gmx程序，該程序?qū)μ囟ńY(jié)構(gòu)塊（如氨基酸丿建立拓?fù)鋽?shù)據(jù)。它需要使用結(jié)構(gòu)塊庫，如果結(jié)構(gòu)里面有庫中沒有的結(jié)構(gòu)單元，轉(zhuǎn)換過程就會(huì)失敗。輸入以下命令來進(jìn)行結(jié)構(gòu)轉(zhuǎn)換，選擇GROMOS53a6力場(chǎng)和SPC水分子模型.選項(xiàng)-gnh表示起始文件中的氫原子會(huì)先被除掉，然后在相應(yīng)的力場(chǎng)中重建。仔細(xì)閱讀屛幕提示，并檢查對(duì)于組氨酸質(zhì)子化所做的選擇，注意蛋白質(zhì)總電荷數(shù)；也要瀏覽輸入文件（protein.pdb,pdb格式丿和輸出文件（protein.gro,gromacs格式丿。注意兩者的區(qū)別；最后瀏覽拓?fù)湮募捌浣Y(jié)構(gòu)。三、結(jié)構(gòu)的能量最小化（真空中丿現(xiàn)在，一個(gè)適合于所選力場(chǎng)的正確結(jié)構(gòu)文件己經(jīng)建立好，同時(shí)還得到了相應(yīng)的拓?fù)鋽?shù)據(jù)。然而，由于結(jié)構(gòu)轉(zhuǎn)換中，牽涉到氫原子的刪除和重建，這可能會(huì)帶來一些相互作用力的變化，比如由于兩個(gè)原子的位置靠得太近引起的作用力。因此我們必須對(duì)結(jié)構(gòu)進(jìn)行一次能量?jī)?yōu)化，使之松弛一點(diǎn)。這其實(shí)就是一個(gè)模擬步驟，包括兩個(gè)過程：第一步，結(jié)構(gòu)和拓?fù)鋽?shù)據(jù)被合成在一起，同時(shí)還包含了一些控制參數(shù)。這個(gè)過程對(duì)產(chǎn)生一個(gè)文件，該文件作為輸入文件，用于第二步mdrun程序的動(dòng)力學(xué)模擬。把這個(gè)過程分為兩步的好處是，輸入文件可以傳送到專門的(高性能丿計(jì)算機(jī)網(wǎng)絡(luò)或者超級(jí)計(jì)算機(jī)，在那里完成模擬計(jì)算。然而，一般只在正式動(dòng)力學(xué)模擬中才會(huì)這么做。為了執(zhí)行第一步，需要用到一個(gè).mdp文件，該文件(文件名minim.mdp丿己經(jīng)被做好并放在你們的目錄下面。該文件包含一系列能量最小化的控制參數(shù)，請(qǐng)查看之。注意，文件以"integrator"行開始，表示指定的數(shù)學(xué)模型。本例中，被指定采用最速下降法計(jì)算5000步?，F(xiàn)在，用以下命令合并結(jié)構(gòu)和拓?fù)湮募?，以及一些控制參?shù)：grompp程序?qū)?huì)標(biāo)明此體系存在非零電荷，并將寫入有關(guān)系統(tǒng)和控制參數(shù)的其他信息。該程序還會(huì)產(chǎn)生一個(gè)額外的輸出文件(mdout.mdp)，該文件包含所有的控制參數(shù)設(shè)置.下一步，運(yùn)行mdrun：由于該蛋白只含有1500個(gè)原子，能量?jī)?yōu)化非?？臁＿x項(xiàng)-v將把每步的勢(shì)能和最大勢(shì)能力打印出來，以便于跟蹤能量?jī)?yōu)化過程。-deffnm選項(xiàng)設(shè)定了輸出文件的默認(rèn)文件名，而不需要對(duì)每個(gè)輸出文件都進(jìn)行命名，以減少選項(xiàng)設(shè)置?，F(xiàn)在，結(jié)構(gòu)松弛下來了，我們?cè)撛O(shè)定周期性邊界并添加溶劑。Whichmethodwasusedfortheenergyminimization?Howmanystepswerespecifiedintheparameterfileandhowmanystepsdidtheminimizationtake?Whatcouldcausetheenergyminimizationtostopbeforethespecifiednumberofstepswasreached?Whatisthefinalpotentialenergyofthesystem?ComparethestructuresbeforeandafterenergyminimizationbyloadingthemintoPyMol

四、周期性邊界條件設(shè)定Ererg/Winimized5tructureedirconf0S.truGtureW'ithPfiC在添加溶劑之前，需要選定模擬體系的大致外形。大多數(shù)常見的模擬都是在周期性邊界條件(PBC)下進(jìn)行的，這意味著要定義一個(gè)外形單元框，該框可用于空間填充堆積。這樣，就可以對(duì)一個(gè)無限、周期性系統(tǒng)進(jìn)行模擬，以避免由于邊框造成的邊界效應(yīng)。建立PBC時(shí)，僅有少數(shù)幾個(gè)通用的形狀。我們將采用rhombicdodecahedron,因?yàn)槠鋵?duì)應(yīng)于最優(yōu)堆積空間，因此是自由旋轉(zhuǎn)原子的最佳選擇。為了不產(chǎn)生周期性影像的直接接觸，我們?cè)O(shè)定了蛋白質(zhì)距離邊框的最小距離為1.0，這樣的話，兩個(gè)相鄰的堆積單元的距離就會(huì)大于2.0nm。周期性邊界條件用editconf命令設(shè)定：看看editconf的輸出，注意體積的變化。同時(shí)，也看看文件protein-PBC.gro的最后一行。在gromacs格式文件(.gro)中，最后一行表示溶劑盒子的形狀。我們常常用三斜晶系矩陣表示法，其中前三個(gè)數(shù)字表示對(duì)角線元素(xx,yy,zz)，最后六個(gè)表示非對(duì)角線元素(xy,xz,yx,yz,zx,zy)。Whatisthenewunitcellvolume?五、溶劑條件設(shè)定SQlvsntStructure.T~—SQlvsntStructure.T~—溶劑盒子設(shè)定好了之后，就可添加溶劑了。溶劑模型有好幾種，每種模型多少都與一種力場(chǎng)相對(duì)應(yīng)。GROMOS96力場(chǎng)常用于SPC水分子模型。溶劑添加不需要寫入拓?fù)鋽?shù)據(jù)，但是也可以在拓?fù)鋽?shù)據(jù)中包含溶劑模型。向盒子中添加溶劑用以下命令：genbox-cpprotein-PBC.gro-csspc216?gro-pprotein.top-oprotein-water.gro看看protein.top的結(jié)尾溶劑添加的情況，添加了多少溶劑？Howmanywatermoleculesareaddedtothesystemandwhatvolumedoesthatcorrespondto?Howmanyatomsareinthesystemnow?用以前學(xué)過的命令，將.gro文件轉(zhuǎn)換為.pdb文件。用Pymol程序讀入溶劑化的結(jié)構(gòu)文件(protein-water.pdb)oPymol中，可用“showcel”'命令顯示盒子形狀。showcell這將畫出一個(gè)框架線來顯示盒子的三維空間邊界。這個(gè)三斜晶圖形對(duì)應(yīng)菱形十二面體，但可能不太明顯。另一個(gè)值得注意的事情是，并非所有的溶劑都在盒子，但以磚形排布。非常類似地，蛋白質(zhì)也有部分伸出水外面。Whyisitnotaproblemtohavetheproteinstickingoutofthewater?

六添加相反的離子以平衡體系電荷目前我們己經(jīng)有了溶劑化的蛋白質(zhì)了，但是體糸的凈電荷數(shù)不為零。為了使體糸電中性，我們需要添加一定數(shù)量的相反離子。添加一定離子的步驟，是一個(gè)很好的練習(xí)，程序genion能完成這些任務(wù)，但是它需要一個(gè)輸入文件，例如，一個(gè)包含了結(jié)構(gòu)和拓?fù)鋽?shù)據(jù)的文件。就像前面一樣，我們可以用grompp來產(chǎn)生這樣的文件：grompp-v-fminim.mdp-cprotein-water.gro-pprotein.top-oprotein-water.tpr文件protein~water.tpr既可作為genion程序的輸入文件。我們?cè)O(shè)定了NA+/CL-(-pnameNA+~nnameCL-)的濃度為0.15M卜cone0.15)，并指定特定離子的加入以是體糸電中性卜neutral)。Genion程序?qū)⒃儐栍秒x子取代何種分子，選擇'SOL'。genion-sprotein-water.tpr-oprotein-solvated.gro-conc0.15-neutral-pnameNA+-nnameCL-注意添加離子的數(shù)量，并注意一個(gè)額外的氯離子被加入，以中和體糸電荷。通過將水分子置換為離子，體糸的拓?fù)鋽?shù)據(jù)文件］protein.top)將不正確了(你可以修改拓?fù)湮募瑴p少溶劑分子數(shù)目；同時(shí)加入一行命令，說明添加Na離子的數(shù)量并寫入另一行命令說明Cl離子的數(shù)量丿。

Editthetopologyfileanddecreasethenumberofsolventmolecules.AlsoaddalinespecifyingthenumberofNAionsandalinespecifyingtheamountofCL.Howmanysodiumandchlorideionsareaddedtothesystem?七、溶劑化體系的能量最小化SystEmConfiguration5y5temTopologySIMULATIONCqntroIPurarneTsrsEMSolventSystEmConfiguration5y5temTopologySIMULATION到此為止，模擬糸統(tǒng)己經(jīng)準(zhǔn)備好了。在進(jìn)行正式模擬前，還需要使體糸再次松弛。因?yàn)槿軇┑奶砑右约半x子的置換，可能產(chǎn)生不利的相互作用力，例如，原子重疊、同種電荷的接近，等，因此需要對(duì)體糸再次能量?jī)?yōu)化，其步驟與前面相同：先運(yùn)行g(shù)rompp，再運(yùn)行mdrun：grompp-v-fminim.mdp-cprotein-solvated.gro-pprotein.top-oprotein-EM-solvated?tprmdrun-v-deffnmprotein-EM-solvated-cprotein-EM-solvated.groHowmanystepswerespecifiedintheparameterfileandhowmanystepsdidtheminimizationtake?Whatisthefinalpotentialenergyofthesystem?八、溶劑和氫原子位置的松弛：位置限制性MD為了緩解體系的隨機(jī)分布的力，我們是溶劑與蛋白質(zhì)相適應(yīng)，比如，我們使溶劑自由移動(dòng)，而使得蛋白質(zhì)上面的非氫原子或多或少地固定在相對(duì)位置。這么做的目的是為了確保溶劑的構(gòu)象“適合”蛋白質(zhì)。這一步是正式MD的第一步。這一步的控制參數(shù)在pr.mdp文件中，執(zhí)行時(shí)需要被讀入。看看這個(gè)文件， 注意integrate行和define行。后者用于允許拓?fù)湮募牧鲃?dòng)控制（toallowflowcontrolinthetopologyfiledefine=-DPOSRES將導(dǎo)致定義關(guān)鍵字"POSRES"。查看拓?fù)湮募淖詈笠恍?，看看怎樣進(jìn)行位置限制。采用grompp和mdrun程序運(yùn)行模擬：Remark:editthefilepr.mdpandchangethevalueofgen_seedgrompp-v-fpr.mdp-cprotein-EM-solvated.gro-pprotein.top-oprotein-PR.tprmdrun-v-deffnmprotein-PRWhatisthelengthofthesimulationinpicoseconds?Howistheinelusionofthepositionrestraintshandledinthetopologyfile?你會(huì)很有趣地看到總能量、勢(shì)能和動(dòng)能的變化。在前一步，粒子沒有速度，所以沒有動(dòng)能，因此也沒有溫度?，F(xiàn)在，模擬一開始，原子就被賦予速度因此獲得動(dòng)能。模擬過程的能量信息被保存為不可讀的二進(jìn)制文件格式（.edr文件丿。該文件的信息可用gromacs工具g_energy程序提取。注：在下一步模擬運(yùn)行時(shí)，你能用以下命令并在執(zhí)行過程中回答有關(guān)問題，以節(jié)省時(shí)間。如果這一步不奏效，讓該程序繼續(xù)在這里運(yùn)行，而打開另外一個(gè)終端，運(yùn)行下一步的grompp和mdrun程序，但要注意兩個(gè)終端的都要在正確的文件夾下運(yùn)行。g_energy-fprotein-PR.edr-othermodynamics-PR.xvg你將要在一列能量參數(shù)中進(jìn)行選擇，選擇好的參數(shù)結(jié)果將會(huì)被輸出。輸入與溫度、勢(shì)能、動(dòng)能和總能量相對(duì)應(yīng)的數(shù)字，然后輸入0（零丿，回車。輸出文件（.xvg）可以用xmgr或xmgrace程序查看，這些程序需要預(yù)先安裝好。也可以用Python語言腳本xvg2ascii.py查看，它能在終端窗口上顯示出圖線。注：綠色曲線是總能量，黑色是勢(shì)能，藍(lán)色是溫度，紅色是動(dòng)能。你也可以點(diǎn)擊曲線并填寫“圖標(biāo)”窗口以改變圖標(biāo)。別忘了點(diǎn)擊'accept'以接受改變。Whathappenswiththetemperature?Whathappenstothepotential/kinetic/totalenergyandhowcanthisbeexplained?九、非限制性MD,打開溫度耦合系統(tǒng)現(xiàn)在應(yīng)該相對(duì)松弛了一些，我們就引入溫度并開始對(duì)系統(tǒng)進(jìn)行熱浴耦合吧。換句話說，將在打開溫度耦合的情況下運(yùn)行一小段時(shí)間，讓系統(tǒng)松弛到一個(gè)新的狀態(tài)。下載控制參數(shù)文件nvt.mdp,然后檢查一下其中的溫度耦合參數(shù)；最后用grompp和mdrun運(yùn)行模擬：Therearetwogroupswhichareseparatelycoupledtoaheatbath?Whichgroupsarethatandwhatdoyouthinkisineachofthesegroups?運(yùn)行結(jié)束后，查看一下能量和溫度。再次提醒，你可以在查看能量的同時(shí)，運(yùn)行下一步。數(shù)據(jù)提取方法同前：觀察能量圖。Whathappenstothetemperature?Whathappenstothepotential/kinetic/totalenergyandhowcanthisbeexplained?十、非限制性MD模擬：打開壓力耦合當(dāng)引入溫度并松弛后，就該引入壓力了。導(dǎo)入控制參數(shù)文件npt.mdp,短時(shí)間運(yùn)行有壓力下的模擬。先查看一下此控制參數(shù)文件，找到控制引入壓力的句柄。此處也應(yīng)該注意，因?yàn)槲覀冃枰匦沦x予離子速度，通過gen_veI行實(shí)現(xiàn)。在這行下面，設(shè)定了速度分布（Maxwell分布丿的溫度，并通過gen_seed產(chǎn)生初始值。gen_seed現(xiàn)在設(shè)定為9999，但是會(huì)被改為一個(gè)新的值。MD模擬是最重要的一步，所以每次都以相同的坐標(biāo)體系、速度和相同的參數(shù)開始模擬時(shí)，就會(huì)導(dǎo)致模擬的結(jié)果一模一樣，這不是我們所希望的。當(dāng)你設(shè)定好了控制參數(shù)后，運(yùn)行以下命令進(jìn)行模擬。Remark:editthefilenpt.mdpandchangethevalueofgen_seedgrompp-v-fnpt.mdp-cprotein-NVT?gro-pprotein.top-oprotein-NPT.tprmdrun-v-deffnmprotein-NPT運(yùn)行完成后，再次查看能量和溫度，同時(shí)注意觀察壓力變化。數(shù)據(jù)提取方法同前:g_energy-fprotein-NPT.edr-othermodynamics-NPT.xvg看看能量曲線、溫度和壓力的變化。Whathappenstothetemperature?Whathappenstothepotential/kinetic/totalenergyandhowcanthisbeexplained?最后一個(gè)問題與從有限數(shù)量的粒子系統(tǒng)中提取熱動(dòng)力特性直接有關(guān)。粗略地說，熱動(dòng)力學(xué)指的是大量粒子（例如，幾十億個(gè)而不是幾千個(gè)丿的行為。取大量粒子的平均特性能夠使數(shù)據(jù)波動(dòng)減少，相反，僅對(duì)少量粒子進(jìn)行平均，波動(dòng)就較大。由于我們引入了壓力，系統(tǒng)的密度會(huì)發(fā)生變化，所以從能量文件中提取密度數(shù)據(jù)，方法如下：Howdoesthedensitybehaveovertime?Whydoesthedensityofthesyste