簡介APBS 教程匯編是面向生物信息學(xué)、生物物理學(xué)及計算化學(xué)研究者的靜電計算入門與進(jìn)階指南。教程以“怎樣做”為主線系統(tǒng)講解如何準(zhǔn)備 PQR/XML 格式結(jié)構(gòu)文件、借助 PDB2PQR 完成 PDB 轉(zhuǎn)換并演示利用 VMD、PyMOL、PMV 觀察靜電勢分布隨后依次覆蓋溶劑化能、結(jié)合能、溶劑化力與 pKa 值計算并給出并行與異步時序等降內(nèi)存方案以及將 APBS 接入 MM/PBSA 分子動力學(xué)流程和通過 Gemstone 在線運行的實踐方法。全文配有玻恩離子等實例、方程與圖像清單便于讀者對照操作。資源為單個 PDF 文件大小約 3.37MB閱讀輕量、便于離線保存。目前已有 107 人學(xué)習(xí)下載適合希望系統(tǒng)掌握 APBS 靜電計算全流程、需要解決大分子體系內(nèi)存瓶頸或想將靜電分析融入 MD 工作流的中高級科研用戶。1. 從 PDB 到 PQRAPBS 靜電計算的第一步APBSAdaptive Poisson-Boltzmann Solver解決的是生物大分子周圍和內(nèi)部的靜電性質(zhì)問題。它的核心是求解連續(xù)介質(zhì)下的泊松-玻爾茲曼方程輸入是帶電荷與半徑的分子結(jié)構(gòu)輸出是空間格點上的靜電勢、溶劑化能和受力。真正攔住新手的不是求解器本身而是輸入數(shù)據(jù)PDB 文件只有坐標(biāo)和元素名沒有電荷和半徑。APBS 需要原子點電荷作為源項同時用范德華半徑構(gòu)造溶質(zhì)與溶劑的介電邊界。因此幾乎所有 APBS 管線都從 PQR 或 XML 參數(shù)文件開始。常見做法是先用 PDB2PQR 生成 PQR再用 VMD 或 PyMOL 做可視化最后通過命令行或 Gemstone 跑計算。如果你是做生物信息學(xué)軟件開發(fā)也可以把這套流程寫進(jìn)自動化管線。下面按這條鏈路逐層拆。2. PDB2PQR 與力場選型電荷半徑參數(shù)從哪里來2.1 PQR 與 XML兩種參數(shù)載體的取舍PQR 格式是 PDB 的超集將 occupancy 和 temperature factor 兩列替換為電荷Q與半徑R。例如ATOM 1 I ION 1 0.000 0.000 0.000 1.00 3.00就是一個半徑為 3.0 ?、帶 1 電荷的離子。PQR 足夠簡單任意文本編輯器都能修改缺點是可擴(kuò)展性差。非標(biāo)準(zhǔn)殘基、配體或自定義原子類型需要借助 PDB2PQR 這類工具才能寫入。XML 格式則提供了結(jié)構(gòu)化參數(shù)定義修改原子類型和力場參數(shù)更方便但文件體積大多數(shù)可視化工具對它的支持不如 PQR。在實際項目中我一般只在調(diào)試自定義參數(shù)時用 XML生產(chǎn)計算都用 PQR。2.2 PDB2PQR 的使用從 Web 到命令行PDB2PQR 不只是格式轉(zhuǎn)換器它會補缺失側(cè)鏈原子、加氫、優(yōu)化氫鍵網(wǎng)絡(luò)、處理滴定狀態(tài)甚至直接生成 APBS 輸入文件。最簡單的入口是 Web 服務(wù)輸入四字符 PDB ID如 1FAS、1MYK或上傳本地 PDB。批量處理或集成到腳本時命令行版本更可靠。以 PDB2PQR 1.2.1 兼容語法為例pdb2pqr --ffPARSE \ --chain \ --apbs-input1fas.in \ --with-ph7.0 \ --titration-state-methodpropka \ 1fas.pdb 1fas.pqr這條命令把 1fas.pdb 轉(zhuǎn)為帶 PARSE 參數(shù)的 1fas.pqr同時生成 1fas.in。--chain保留肽鏈邊界--with-ph設(shè)置目標(biāo) pH配合propka方法決定可滴定殘基的質(zhì)子化狀態(tài)--apbs-input讓 PDB2PQR 順便寫好 APBS 輸入文件后續(xù)只需微調(diào)網(wǎng)格和介電參數(shù)。需要留意的是PDB2PQR 不能修復(fù)大段缺失的骨架坐標(biāo)。如果結(jié)構(gòu)里缺失殘基超過 10%先用同源模建補齊再轉(zhuǎn)換。2.3 力場選擇為什么默認(rèn)是 PARSE力場決定了每個原子的電荷和半徑直接影響 PB 方程源項和介電邊界。教程建議“大多數(shù)情況選 PARSE”這個建議至今有效。PARSE 是針對隱式溶劑模型校準(zhǔn)的參數(shù)集在蛋白靜電勢可視化和溶劑化能估算中表現(xiàn)穩(wěn)定。需要核酸參數(shù)、與 CHARMM/AMBER 分子動力學(xué)軌跡對接時才考慮更換力場。力場適用場景注意點PARSE蛋白/核酸靜電勢可視化、溶劑化能初估配體參數(shù)需額外補充AMBER與 AMBER MD 軌跡配套、MM/PBSA力場版本需與 MD 一致CHARMMCHARMM 體系、脂膜、核酸原子類型命名差異大自定義非天然殘基、藥物分子手寫 XML 或力場文件還有一個容易踩的坑PDB2PQR 1.2.0 之前對配體幾乎沒有支持分子對接或酶-抑制劑體系經(jīng)常在配體這一步失敗。新版雖然能處理常見配體但金屬中心、共價配體仍需手動驗證電荷。對這類體系我一般先用 Gaussian 或 ORCA 在 HF/6-31G* 級別擬合 ESP 電荷再寫自定義力場。2.4 輸出命名與額外處理輸出命名配置雖然不影響 APBS 能量計算但會直接影響 VMD/PyMOL 的顯示效果。PDB2PQR 默認(rèn)使用 internal 命名與 IUPAC 標(biāo)準(zhǔn)一致能避免 VMD 中奇怪的成鍵。如果結(jié)構(gòu)里的原子名已經(jīng)被 MD 工具改過保留原命名可能更安全。加氫方式也在這里決定pH 7.0 下組氨酸可帶正電天冬氨酸和谷氨酸去質(zhì)子化。propka方法會比簡單按模型 pKa 更精細(xì)會考慮周圍殘基的靜電影響。生成 PQR 后可以用grep快速檢查總電荷是否合理比如一個普通蛋白在中性 pH 下總電荷在 -20 到 20 之間。3. APBS 輸入文件拆解玻恩離子算例與溶劑化能3.1 輸入文件的三段式結(jié)構(gòu)APBS 輸入文件由read、elec、print三類塊組成。read讀入 PQR/XML 分子elec定義一次靜電計算print把多個elec的能量相減。溶劑化能計算通常需要兩個elec塊溶劑態(tài)溶質(zhì)介電 pdie1.0溶劑介電 sdie78.54和參考態(tài)pdie1.0sdie1.0。APBS 返回的充電自由能中包含分子的自相互作用能這項對網(wǎng)格離散化非常敏感因此兩個塊必須使用完全相同的dime、grid、gcent讓自作用能精確抵消。3.2 網(wǎng)格與邊界參數(shù)怎么調(diào)網(wǎng)格參數(shù)決定計算精度和內(nèi)存。dime是每個方向的格點數(shù)grid是格點間距gcent mol 1表示網(wǎng)格中心在分子 1。APBS 推薦溶質(zhì)到盒子邊界至少留 10 ?。如果網(wǎng)格太小邊界條件會污染能量網(wǎng)格太大內(nèi)存占用又按三次方膨脹。關(guān)鍵參數(shù)如下參數(shù)作用常見取值dime網(wǎng)格維數(shù)97 / 129grid格點間距?0.33 / 0.50pdie溶質(zhì)介電常數(shù)1.0 - 4.0sdie溶劑介電常數(shù)78.54水bcfl邊界條件mdh / mgautochgm電荷離散方式spl2srfm表面定義mol / spl2小分子或單離子用dime 97足夠2 萬原子以上的蛋白通常需要dime 129。bcfl mdh是多重極展開邊界精度高分子很大時可以考慮bcfl mgauto。srfm mol生成分子表面srad 1.4是水探針半徑。3.3 玻恩離子算例極性溶劑化能計算玻恩離子是最簡單的驗證體系一個半徑 3 ?、電荷 1 的球在均勻水介質(zhì)中理論極性溶劑化能約 -230.62 kJ/mol。PQR 文件只有一行原子記錄REMARK Born ion, radius 3 A, charge 1 e ATOM 1 I ION 1 0.000 0.000 0.000 1.00 3.00 ENDAPBS 輸入文件如下read mol pqr born.pqr end elec name solv mg-manual dime 97 97 97 nlev 4 grid 0.33 0.33 0.33 gcent mol 1 mol 1 lpbe bcfl mdh pdie 1.0 sdie 78.54 chgm spl2 srfm mol srad 1.4 swin 0.3 sdens 10.0 temp 298.15 calcenergy total calcforce no end elec name ref mg-manual dime 97 97 97 nlev 4 grid 0.33 0.33 0.33 gcent mol 1 mol 1 lpbe bcfl mdh pdie 1.0 sdie 1.0 chgm spl2 srfm mol srad 1.4 swin 0.3 sdens 10.0 temp 298.15 calcenergy total calcforce no end print energy solv - ref quit運行apbs born.in結(jié)果約 -229.59 kJ/mol與理論值非常接近。lpbe求解線性化 PB 方程chgm spl2用樣條函數(shù)把點電荷分散到格點比chgm spl0更平滑。calcenergy total讓每個塊輸出總能量最后print energy solv - ref做減法。如果想把能量換成 kcal/mol需要自己乘 0.239。3.4 非極性溶劑化Wagoner-Baker 模型參數(shù)非極性溶劑化包含空穴形成和溶質(zhì)-溶劑分散吸引兩部分。APBS 實現(xiàn)了 Wagoner-Baker 模型空穴項由press壓強和gamma表面張力決定分散項由bconc溶劑密度和 Lennard-Jones 參數(shù)決定。將press和bconc設(shè)為 0 可以退化為純 SASA 模型。實際計算里gamma常用 0.00542 kJ/mol/?2但這是經(jīng)驗值不同體系需要重標(biāo)定。在elec塊中追加press 0.0 gamma 0.00542 bconc 0.0333若需要溶劑化力只要把calcforce no改成calcforce compsAPBS 就會輸出每個原子的受力。前提是網(wǎng)格已經(jīng)收斂否則力的噪聲會非常大。4. 結(jié)合能與 pKa自由能循環(huán)在 APBS 中的落法4.1 結(jié)合自由能循環(huán)拆解隱式溶劑模型下結(jié)合自由能可以分成溶劑化貢獻(xiàn)和庫侖貢獻(xiàn)。溶劑化貢獻(xiàn)由復(fù)合物和分離組分的溶劑化能差給出需要分別按第三章的流程計算復(fù)合物、受體、配體三個體系的極性和非極性項。庫侖貢獻(xiàn)則是在統(tǒng)一介電常數(shù)下計算復(fù)合物與受體、配體兩兩庫侖作用的差值。這樣做的前提是庫侖作用的介電常數(shù)要與溶劑化能參考態(tài)一致。如果溶劑化能計算中 pdie1.0那么庫侖能也要除以 1.0如果 pdie4.0則除以 4.0。4.2 用 tools/manip/coulomb 算庫侖貢獻(xiàn)APBS 安裝包附帶的tools/manip/coulomb直接從 PQR 文件計算真空庫侖能介電常數(shù) 1。對復(fù)合物、受體、配體分別運行tools/manip/coulomb complex.pqr complex.coul tools/manip/coulomb receptor.pqr receptor.coul tools/manip/coulomb ligand.pqr ligand.coul輸出文件的最后一行是總庫侖能。結(jié)合庫侖貢獻(xiàn)為dG_coul E_complex - E_receptor - E_ligand注意三個 PQR 必須來自同一個構(gòu)象至少坐標(biāo)系要對齊。常見的錯誤是直接從 PDB 里分別提取受體和配體坐標(biāo)把復(fù)合物中的構(gòu)象變化忽略了這會帶來不小的誤差。對于柔性對接應(yīng)該從復(fù)合物軌跡中提取同一時刻的坐標(biāo)再拆開。4.3 配體參數(shù)缺失最常見的失敗路徑“配體沒有設(shè)置參數(shù)”幾乎每個人都會遇到。PDB2PQR 1.2.0 之后雖然能處理配體但依賴力場庫里已有的殘基定義。只要配體不是標(biāo)準(zhǔn)氨基酸、核酸或常見輔因子就可能報Unknown atom type。常見做法有兩種用 AMBER 的antechamber生成 GAFF 參數(shù)再轉(zhuǎn) PQR或者把配體單獨寫成自定義 XML 參數(shù)文件與蛋白 PQR 合并。后者更可控項目里做藥物篩選時我傾向用這種方式。4.4 pKa 計算模型化合物與溶菌酶案例pKa 計算是溶劑化能的進(jìn)階應(yīng)用。蛋白質(zhì)中一個可滴定基團(tuán)的 pKa 相對模型值發(fā)生偏移本質(zhì)上是帶電態(tài)與中性態(tài)在蛋白環(huán)境和水環(huán)境間的遷移自由能差。APBS 方案是對每個可滴定殘基準(zhǔn)備質(zhì)子化/去質(zhì)子化兩套 PQR分別算蛋白環(huán)境和水環(huán)境的溶劑化能得到遷移自由能再換算出 pKa 偏移。模型化合物 pKa 是參照系常用值如下基團(tuán)模型 pKaArginine13.0Aspartic acid4.0Cysteine8.7C-terminus3.8Glutamic acid4.4Histidine6.3Lysine10.4N-terminus8.0Tyrosine9.6溶菌酶HEWL是教程里的經(jīng)典驗證體系?；钚晕稽c Glu35 與 Asp52 的實驗 pKa 分別偏移到約 6.2 和 3.5用 APBS 能定性重現(xiàn)這個趨勢。要復(fù)現(xiàn)結(jié)構(gòu)必須先加氫網(wǎng)格收斂要足夠好。pKa 對去溶劑化能極度敏感1 kJ/mol 誤差會帶來約 0.2 個 pH 單位的偏移。我習(xí)慣先跑一遍第三章的溶劑化能收斂測試再跑 pKa 掃描。5. 可視化、并行計算與排錯把 APBS 變成可復(fù)現(xiàn)管線5.1 在 VMD 和 PyMOL 里做同樣的勢能圖VMD 的完整操作載入 PQR 后從主菜單打開Extensions → Analysis → APBS Electrostatics在設(shè)置里把 APBS 二進(jìn)制路徑指對。計算完成后用Graphics → Representations新建 Rep繪圖方式選Isosurfacedraw 改為Solid Surface透明顯示等勢面值設(shè)為 1正值藍(lán)色和 -1負(fù)值紅色。如果只想看分子表面電荷分布把繪圖方式改為Surf著色方式選Volume色階范圍設(shè)到 -10 到 10kT/e顏色標(biāo)度在Graphics → Colors → Color Scale里選 RWB。PyMOL 路線更短Plugin → APBS ToolsMain 標(biāo)簽選擇Use another PQR指定 PDB2PQR 生成的文件Configuration 標(biāo)簽里先Set grid再Run APBS。跑完切到 Visualization 標(biāo)簽點Update載入勢能然后分別設(shè)置正負(fù)等勢面閾值。常用閾值是 ±1、±5、±10 kT/e。注意 PyMOL 早期版本在 contour 值上有單位換算小誤差顯示結(jié)果和 APBS 原始數(shù)據(jù)可能有量級差異建議使用較新版本。5.2 讓大分子計算不爆內(nèi)存dime 129對 2 萬原子體系內(nèi)存約幾個 GB。若內(nèi)存吃緊先用mg-auto代替mg-manual讓 APBS 自動生成網(wǎng)格。仍然不夠就放寬grid從 0.33 到 0.5 ?但前提是能量隨網(wǎng)格變化已經(jīng)收斂。并行計算用 MPImpirun -np 8 apbs big.innlev是多層網(wǎng)格層數(shù)一般 4 就夠了。排錯時優(yōu)先看輸出里的 Warning提示grid too small就擴(kuò)大盒子提示charges exceed就檢查 PQR 電荷總和。想快速驗證參數(shù)可以用 Gemstone 在線跑免去本地安裝。最后把以上步驟寫進(jìn) Makefile 或 Snakemake就是一條可復(fù)現(xiàn)的從 PDB 到靜電勢圖的 APBS 管線。記得把缺失原子檢查和 grid 收斂測試寫進(jìn)去換機器、換 APBS 版本時先跑一遍回歸。本文還有配套的精品資源點擊獲取 從PDB到PQR：APBS靜電勢與溶劑化能計算全解析
簡介APBS 教程匯編是面向生物信息學(xué)、生物物理學(xué)及計算化學(xué)研究者的靜電計算入門與進(jìn)階指南。教程以“怎樣做”為主線系統(tǒng)講解如何準(zhǔn)備 PQR/XML 格式結(jié)構(gòu)文件、借助 PDB2PQR 完成 PDB 轉(zhuǎn)換并演示利用 VMD、PyMOL、PMV 觀察靜電勢分布隨后依次覆蓋溶劑化能、結(jié)合能、溶劑化力與 pKa 值計算并給出并行與異步時序等降內(nèi)存方案以及將 APBS 接入 MM/PBSA 分子動力學(xué)流程和通過 Gemstone 在線運行的實踐方法。全文配有玻恩離子等實例、方程與圖像清單便于讀者對照操作。資源為單個 PDF 文件大小約 3.37MB閱讀輕量、便于離線保存。目前已有 107 人學(xué)習(xí)下載適合希望系統(tǒng)掌握 APBS 靜電計算全流程、需要解決大分子體系內(nèi)存瓶頸或想將靜電分析融入 MD 工作流的中高級科研用戶。1. 從 PDB 到 PQRAPBS 靜電計算的第一步APBSAdaptive Poisson-Boltzmann Solver解決的是生物大分子周圍和內(nèi)部的靜電性質(zhì)問題。它的核心是求解連續(xù)介質(zhì)下的泊松-玻爾茲曼方程輸入是帶電荷與半徑的分子結(jié)構(gòu)輸出是空間格點上的靜電勢、溶劑化能和受力。真正攔住新手的不是求解器本身而是輸入數(shù)據(jù)PDB 文件只有坐標(biāo)和元素名沒有電荷和半徑。APBS 需要原子點電荷作為源項同時用范德華半徑構(gòu)造溶質(zhì)與溶劑的介電邊界。因此幾乎所有 APBS 管線都從 PQR 或 XML 參數(shù)文件開始。常見做法是先用 PDB2PQR 生成 PQR再用 VMD 或 PyMOL 做可視化最后通過命令行或 Gemstone 跑計算。如果你是做生物信息學(xué)軟件開發(fā)也可以把這套流程寫進(jìn)自動化管線。下面按這條鏈路逐層拆。2. PDB2PQR 與力場選型電荷半徑參數(shù)從哪里來2.1 PQR 與 XML兩種參數(shù)載體的取舍PQR 格式是 PDB 的超集將 occupancy 和 temperature factor 兩列替換為電荷Q與半徑R。例如ATOM 1 I ION 1 0.000 0.000 0.000 1.00 3.00就是一個半徑為 3.0 ?、帶 1 電荷的離子。PQR 足夠簡單任意文本編輯器都能修改缺點是可擴(kuò)展性差。非標(biāo)準(zhǔn)殘基、配體或自定義原子類型需要借助 PDB2PQR 這類工具才能寫入。XML 格式則提供了結(jié)構(gòu)化參數(shù)定義修改原子類型和力場參數(shù)更方便但文件體積大多數(shù)可視化工具對它的支持不如 PQR。在實際項目中我一般只在調(diào)試自定義參數(shù)時用 XML生產(chǎn)計算都用 PQR。2.2 PDB2PQR 的使用從 Web 到命令行PDB2PQR 不只是格式轉(zhuǎn)換器它會補缺失側(cè)鏈原子、加氫、優(yōu)化氫鍵網(wǎng)絡(luò)、處理滴定狀態(tài)甚至直接生成 APBS 輸入文件。最簡單的入口是 Web 服務(wù)輸入四字符 PDB ID如 1FAS、1MYK或上傳本地 PDB。批量處理或集成到腳本時命令行版本更可靠。以 PDB2PQR 1.2.1 兼容語法為例pdb2pqr --ffPARSE \ --chain \ --apbs-input1fas.in \ --with-ph7.0 \ --titration-state-methodpropka \ 1fas.pdb 1fas.pqr這條命令把 1fas.pdb 轉(zhuǎn)為帶 PARSE 參數(shù)的 1fas.pqr同時生成 1fas.in。--chain保留肽鏈邊界--with-ph設(shè)置目標(biāo) pH配合propka方法決定可滴定殘基的質(zhì)子化狀態(tài)--apbs-input讓 PDB2PQR 順便寫好 APBS 輸入文件后續(xù)只需微調(diào)網(wǎng)格和介電參數(shù)。需要留意的是PDB2PQR 不能修復(fù)大段缺失的骨架坐標(biāo)。如果結(jié)構(gòu)里缺失殘基超過 10%先用同源模建補齊再轉(zhuǎn)換。2.3 力場選擇為什么默認(rèn)是 PARSE力場決定了每個原子的電荷和半徑直接影響 PB 方程源項和介電邊界。教程建議“大多數(shù)情況選 PARSE”這個建議至今有效。PARSE 是針對隱式溶劑模型校準(zhǔn)的參數(shù)集在蛋白靜電勢可視化和溶劑化能估算中表現(xiàn)穩(wěn)定。需要核酸參數(shù)、與 CHARMM/AMBER 分子動力學(xué)軌跡對接時才考慮更換力場。力場適用場景注意點PARSE蛋白/核酸靜電勢可視化、溶劑化能初估配體參數(shù)需額外補充AMBER與 AMBER MD 軌跡配套、MM/PBSA力場版本需與 MD 一致CHARMMCHARMM 體系、脂膜、核酸原子類型命名差異大自定義非天然殘基、藥物分子手寫 XML 或力場文件還有一個容易踩的坑PDB2PQR 1.2.0 之前對配體幾乎沒有支持分子對接或酶-抑制劑體系經(jīng)常在配體這一步失敗。新版雖然能處理常見配體但金屬中心、共價配體仍需手動驗證電荷。對這類體系我一般先用 Gaussian 或 ORCA 在 HF/6-31G* 級別擬合 ESP 電荷再寫自定義力場。2.4 輸出命名與額外處理輸出命名配置雖然不影響 APBS 能量計算但會直接影響 VMD/PyMOL 的顯示效果。PDB2PQR 默認(rèn)使用 internal 命名與 IUPAC 標(biāo)準(zhǔn)一致能避免 VMD 中奇怪的成鍵。如果結(jié)構(gòu)里的原子名已經(jīng)被 MD 工具改過保留原命名可能更安全。加氫方式也在這里決定pH 7.0 下組氨酸可帶正電天冬氨酸和谷氨酸去質(zhì)子化。propka方法會比簡單按模型 pKa 更精細(xì)會考慮周圍殘基的靜電影響。生成 PQR 后可以用grep快速檢查總電荷是否合理比如一個普通蛋白在中性 pH 下總電荷在 -20 到 20 之間。3. APBS 輸入文件拆解玻恩離子算例與溶劑化能3.1 輸入文件的三段式結(jié)構(gòu)APBS 輸入文件由read、elec、print三類塊組成。read讀入 PQR/XML 分子elec定義一次靜電計算print把多個elec的能量相減。溶劑化能計算通常需要兩個elec塊溶劑態(tài)溶質(zhì)介電 pdie1.0溶劑介電 sdie78.54和參考態(tài)pdie1.0sdie1.0。APBS 返回的充電自由能中包含分子的自相互作用能這項對網(wǎng)格離散化非常敏感因此兩個塊必須使用完全相同的dime、grid、gcent讓自作用能精確抵消。3.2 網(wǎng)格與邊界參數(shù)怎么調(diào)網(wǎng)格參數(shù)決定計算精度和內(nèi)存。dime是每個方向的格點數(shù)grid是格點間距gcent mol 1表示網(wǎng)格中心在分子 1。APBS 推薦溶質(zhì)到盒子邊界至少留 10 ?。如果網(wǎng)格太小邊界條件會污染能量網(wǎng)格太大內(nèi)存占用又按三次方膨脹。關(guān)鍵參數(shù)如下參數(shù)作用常見取值dime網(wǎng)格維數(shù)97 / 129grid格點間距?0.33 / 0.50pdie溶質(zhì)介電常數(shù)1.0 - 4.0sdie溶劑介電常數(shù)78.54水bcfl邊界條件mdh / mgautochgm電荷離散方式spl2srfm表面定義mol / spl2小分子或單離子用dime 97足夠2 萬原子以上的蛋白通常需要dime 129。bcfl mdh是多重極展開邊界精度高分子很大時可以考慮bcfl mgauto。srfm mol生成分子表面srad 1.4是水探針半徑。3.3 玻恩離子算例極性溶劑化能計算玻恩離子是最簡單的驗證體系一個半徑 3 ?、電荷 1 的球在均勻水介質(zhì)中理論極性溶劑化能約 -230.62 kJ/mol。PQR 文件只有一行原子記錄REMARK Born ion, radius 3 A, charge 1 e ATOM 1 I ION 1 0.000 0.000 0.000 1.00 3.00 ENDAPBS 輸入文件如下read mol pqr born.pqr end elec name solv mg-manual dime 97 97 97 nlev 4 grid 0.33 0.33 0.33 gcent mol 1 mol 1 lpbe bcfl mdh pdie 1.0 sdie 78.54 chgm spl2 srfm mol srad 1.4 swin 0.3 sdens 10.0 temp 298.15 calcenergy total calcforce no end elec name ref mg-manual dime 97 97 97 nlev 4 grid 0.33 0.33 0.33 gcent mol 1 mol 1 lpbe bcfl mdh pdie 1.0 sdie 1.0 chgm spl2 srfm mol srad 1.4 swin 0.3 sdens 10.0 temp 298.15 calcenergy total calcforce no end print energy solv - ref quit運行apbs born.in結(jié)果約 -229.59 kJ/mol與理論值非常接近。lpbe求解線性化 PB 方程chgm spl2用樣條函數(shù)把點電荷分散到格點比chgm spl0更平滑。calcenergy total讓每個塊輸出總能量最后print energy solv - ref做減法。如果想把能量換成 kcal/mol需要自己乘 0.239。3.4 非極性溶劑化Wagoner-Baker 模型參數(shù)非極性溶劑化包含空穴形成和溶質(zhì)-溶劑分散吸引兩部分。APBS 實現(xiàn)了 Wagoner-Baker 模型空穴項由press壓強和gamma表面張力決定分散項由bconc溶劑密度和 Lennard-Jones 參數(shù)決定。將press和bconc設(shè)為 0 可以退化為純 SASA 模型。實際計算里gamma常用 0.00542 kJ/mol/?2但這是經(jīng)驗值不同體系需要重標(biāo)定。在elec塊中追加press 0.0 gamma 0.00542 bconc 0.0333若需要溶劑化力只要把calcforce no改成calcforce compsAPBS 就會輸出每個原子的受力。前提是網(wǎng)格已經(jīng)收斂否則力的噪聲會非常大。4. 結(jié)合能與 pKa自由能循環(huán)在 APBS 中的落法4.1 結(jié)合自由能循環(huán)拆解隱式溶劑模型下結(jié)合自由能可以分成溶劑化貢獻(xiàn)和庫侖貢獻(xiàn)。溶劑化貢獻(xiàn)由復(fù)合物和分離組分的溶劑化能差給出需要分別按第三章的流程計算復(fù)合物、受體、配體三個體系的極性和非極性項。庫侖貢獻(xiàn)則是在統(tǒng)一介電常數(shù)下計算復(fù)合物與受體、配體兩兩庫侖作用的差值。這樣做的前提是庫侖作用的介電常數(shù)要與溶劑化能參考態(tài)一致。如果溶劑化能計算中 pdie1.0那么庫侖能也要除以 1.0如果 pdie4.0則除以 4.0。4.2 用 tools/manip/coulomb 算庫侖貢獻(xiàn)APBS 安裝包附帶的tools/manip/coulomb直接從 PQR 文件計算真空庫侖能介電常數(shù) 1。對復(fù)合物、受體、配體分別運行tools/manip/coulomb complex.pqr complex.coul tools/manip/coulomb receptor.pqr receptor.coul tools/manip/coulomb ligand.pqr ligand.coul輸出文件的最后一行是總庫侖能。結(jié)合庫侖貢獻(xiàn)為dG_coul E_complex - E_receptor - E_ligand注意三個 PQR 必須來自同一個構(gòu)象至少坐標(biāo)系要對齊。常見的錯誤是直接從 PDB 里分別提取受體和配體坐標(biāo)把復(fù)合物中的構(gòu)象變化忽略了這會帶來不小的誤差。對于柔性對接應(yīng)該從復(fù)合物軌跡中提取同一時刻的坐標(biāo)再拆開。4.3 配體參數(shù)缺失最常見的失敗路徑“配體沒有設(shè)置參數(shù)”幾乎每個人都會遇到。PDB2PQR 1.2.0 之后雖然能處理配體但依賴力場庫里已有的殘基定義。只要配體不是標(biāo)準(zhǔn)氨基酸、核酸或常見輔因子就可能報Unknown atom type。常見做法有兩種用 AMBER 的antechamber生成 GAFF 參數(shù)再轉(zhuǎn) PQR或者把配體單獨寫成自定義 XML 參數(shù)文件與蛋白 PQR 合并。后者更可控項目里做藥物篩選時我傾向用這種方式。4.4 pKa 計算模型化合物與溶菌酶案例pKa 計算是溶劑化能的進(jìn)階應(yīng)用。蛋白質(zhì)中一個可滴定基團(tuán)的 pKa 相對模型值發(fā)生偏移本質(zhì)上是帶電態(tài)與中性態(tài)在蛋白環(huán)境和水環(huán)境間的遷移自由能差。APBS 方案是對每個可滴定殘基準(zhǔn)備質(zhì)子化/去質(zhì)子化兩套 PQR分別算蛋白環(huán)境和水環(huán)境的溶劑化能得到遷移自由能再換算出 pKa 偏移。模型化合物 pKa 是參照系常用值如下基團(tuán)模型 pKaArginine13.0Aspartic acid4.0Cysteine8.7C-terminus3.8Glutamic acid4.4Histidine6.3Lysine10.4N-terminus8.0Tyrosine9.6溶菌酶HEWL是教程里的經(jīng)典驗證體系?；钚晕稽c Glu35 與 Asp52 的實驗 pKa 分別偏移到約 6.2 和 3.5用 APBS 能定性重現(xiàn)這個趨勢。要復(fù)現(xiàn)結(jié)構(gòu)必須先加氫網(wǎng)格收斂要足夠好。pKa 對去溶劑化能極度敏感1 kJ/mol 誤差會帶來約 0.2 個 pH 單位的偏移。我習(xí)慣先跑一遍第三章的溶劑化能收斂測試再跑 pKa 掃描。5. 可視化、并行計算與排錯把 APBS 變成可復(fù)現(xiàn)管線5.1 在 VMD 和 PyMOL 里做同樣的勢能圖VMD 的完整操作載入 PQR 后從主菜單打開Extensions → Analysis → APBS Electrostatics在設(shè)置里把 APBS 二進(jìn)制路徑指對。計算完成后用Graphics → Representations新建 Rep繪圖方式選Isosurfacedraw 改為Solid Surface透明顯示等勢面值設(shè)為 1正值藍(lán)色和 -1負(fù)值紅色。如果只想看分子表面電荷分布把繪圖方式改為Surf著色方式選Volume色階范圍設(shè)到 -10 到 10kT/e顏色標(biāo)度在Graphics → Colors → Color Scale里選 RWB。PyMOL 路線更短Plugin → APBS ToolsMain 標(biāo)簽選擇Use another PQR指定 PDB2PQR 生成的文件Configuration 標(biāo)簽里先Set grid再Run APBS。跑完切到 Visualization 標(biāo)簽點Update載入勢能然后分別設(shè)置正負(fù)等勢面閾值。常用閾值是 ±1、±5、±10 kT/e。注意 PyMOL 早期版本在 contour 值上有單位換算小誤差顯示結(jié)果和 APBS 原始數(shù)據(jù)可能有量級差異建議使用較新版本。5.2 讓大分子計算不爆內(nèi)存dime 129對 2 萬原子體系內(nèi)存約幾個 GB。若內(nèi)存吃緊先用mg-auto代替mg-manual讓 APBS 自動生成網(wǎng)格。仍然不夠就放寬grid從 0.33 到 0.5 ?但前提是能量隨網(wǎng)格變化已經(jīng)收斂。并行計算用 MPImpirun -np 8 apbs big.innlev是多層網(wǎng)格層數(shù)一般 4 就夠了。排錯時優(yōu)先看輸出里的 Warning提示grid too small就擴(kuò)大盒子提示charges exceed就檢查 PQR 電荷總和。想快速驗證參數(shù)可以用 Gemstone 在線跑免去本地安裝。最后把以上步驟寫進(jìn) Makefile 或 Snakemake就是一條可復(fù)現(xiàn)的從 PDB 到靜電勢圖的 APBS 管線。記得把缺失原子檢查和 grid 收斂測試寫進(jìn)去換機器、換 APBS 版本時先跑一遍回歸。本文還有配套的精品資源點擊獲取