• 巨幅降低Multiwfn結合VMD繪制分子表面靜電勢圖耗時的一個關鍵技巧

    巨幅降低Multiwfn結合VMD繪制分子表面靜電勢圖耗時的一個關鍵技巧

    文/Sobereva@北京科音  2021-Jun-30


    1 前言

    Multiwfn(http://www.shanxitv.org/multiwfn)有十分強大的靜電勢分析功能,按照《使用Multiwfn+VMD快速地繪制靜電勢著色的分子范德華表面圖和分子間穿透圖》(http://www.shanxitv.org/443)的做法還可以非常方便地結合VMD繪制效果很好的分子表面靜電勢圖,此方法已經被大量文章所使用。此文中有一個腳本是ESPiso.bat,它會調用Multiwfn來分別計算電子密度和靜電勢的格點數據并導出為cub格式,之后就可以靠VMD繪制靜電勢著色的電子密度等值面圖。文中提到,由于計算靜電勢格點數據耗時很高,所以強烈建議讓Multiwfn調用Gaussian自帶的cubegen程序來實現這一步,因為cubegen計算靜電勢格點數據比Multiwfn內部代碼速度更快。

    有一日筆者突然想到一個主意:VMD里將靜電勢映射到電子密度等值面上的時候是利用電子密度等值面附近格點上的靜電勢插值來實現的,這些被實際利用到的格點只不過占所有格點的很小比例,如果在Multiwfn計算靜電勢的時候忽略掉其它格點(即不計算距離范德華表面較遠的格點的靜電勢),不就能巨幅節約計算靜電勢格點數據步驟的耗時?筆者遂將這個想法實現進了Multiwfn,發現耗時降低了一個數量級,以前很難算得動的體系現在都能算得很輕松,而且比調用cubegen的時候耗時還低得多!下面就對這個極為重要的節約耗時的方法做一下具體說明。讀者必須用2021-Jun-30及以后更新的Multiwfn版本。


    2 設置方法

    在Multiwfn的settings.ini里有一個參數ESPrhoiso,默認為0,代表不使用這種降低耗時的方法。如果設為某個不為0的值,比如設為0.001,就代表在通過主功能5計算靜電勢格點數據的時候,在0.001 a.u.電子密度等值面附近的格點才被計算靜電勢,而其它格點的靜電勢直接被當做0。為了方便,這個參數也可以通過命令行來設置,比如可以寫Multiwfn nico.wfn -ESPrhoiso 0.001。用命令行參數設置的優先級高于settings.ini里的設置。

    由于這個做法節約耗時效果極佳,而且對所得圖像質量沒有任何不良影響,因此從2021-Jun-30更新的Multiwfn開始,在自帶的ESPiso.bat和ESPiso.sh腳本里已經默認加了-ESPrhoiso 0.001選項。因此在計算靜電勢格點數據的時候你會從Multiwfn的窗口中看到類似以下提示:
    Note: ESP will be calculated only for the grids around isosurface of electron density of 0.001000 a.u.
    Detecting the grids for calculating ESP...
    Number of grids to calculate ESP:       56473
    這里Multiwfn首先計算了電子密度格點數據,然后判斷出共有56473個格點在電子密度0.001 a.u.等值面附近,因此需要之后被計算靜電勢。

    顯然,使用《使用Multiwfn+VMD快速地繪制靜電勢著色的分子范德華表面圖和分子間穿透圖》(http://www.shanxitv.org/443)里的ESPiso.vmd腳本在VMD中繪圖的話,你想繪制哪個電子密度等值面的靜電勢著色圖,ESPrhoiso就應該設多少。比如你設ESPrhoiso為0.001,那么若在VMD里把電子密度等值面手動改為0.002,等值面顏色就會很古怪了。

    使用這種節約耗時的做法時不能讓Multiwfn調用cubegen來算靜電勢格點數據,否則沒有效果。

    下面把Multiwfn的上述節約耗時的做法的細節交代一下,感興趣的讀者可以了解一下,一般用戶不用看。在settings.ini里有個ESPrhonlay參數,下面簡寫為n,其值默認為1。Multiwfn會對每個格點進行檢測,看看圍繞它周圍的n層格點里是否電子密度有的大于ESPrhoiso而有的小于ESPrhoiso,如果是的話,則這個格點就被視為電子密度等值面附近的格點(邊界格點),故需要被計算靜電勢。例如n=1,就需要對當前格點周圍26個格點進行檢驗,如果n=2就需要對周圍兩層厚度的共124個格點進行檢驗。實際發現n=2時被判斷為邊界格點的數目比n=1時大約多一倍,因此靜電勢計算耗時也多一倍左右。n=2的時候是絕對不犧牲所得圖像質量的,而用更便宜的n=1的時候往往會造成VMD里顯示的分子表面靜電勢圖在個別地方顏色有輕微瑕疵(輕微發白,白色對應靜電勢為0的情況),這大抵是VMD對靜電勢插值時利用到了與邊界格點相鄰的沒有計算靜電勢的格點所致,而這樣的格點的靜電勢被當成0,和周圍邊界格點的靜電勢相差太大。為了彌補這個問題,Multiwfn自動會在計算完所有邊界格點的靜電勢后,對所有非邊界格點進行循環,如果發現某個非邊界格點的上下左右前后的格點里存在邊界格點的話,就用它的靜電勢當做當前這個非邊界格點的靜電勢。實測n=1在這樣處理情況下得到的圖像效果和n=2沒有任何區別,而比用n=2便宜一倍。因此ESPrhonlay參數默認是1,這是非常安全的,不損害圖像視覺效果。


    3 實際加速效果測試

    3.1 小分子:多巴胺

    這里對實際耗時進行對比測試。第一個體系是下圖所示的多巴胺,共22個原子,波函數是B3LYP/6-311G**計算得到的,其表面靜電勢圖如下

    此例筆者用Intel i7-10870H(8核)普通家用CPU在Win10 64bit下進行測試,cubegen是G16 B.01自帶的,測試的是Medium quality grid的設置下計算靜電勢格點數據:

    Multiwfn內部代碼算靜電勢:179秒(對均勻分布的共98*82*67=538412個點計算了靜電勢)
    Multiwfn內部代碼算靜電勢 -ESPrhoiso 0.001:10秒(只對等值面附近28374個點計算了靜電勢)
    Multiwfn調用cubegen算靜電勢:33秒(對均勻分布的共98*82*67=538412個點計算了靜電勢)

    可見,如果以默認情況,即用Multiwfn內部代碼計算矩形盒子里所有格點的靜電勢,即便對這么個小分子耗時也是相當高的,花了三分鐘。而使用ESPrhoiso設置后,由于被計算的28374個點只占全部538412個點的5%,耗時低了約20倍!比調用cubegen耗時還低得多!

    如果直接用《使用Multiwfn+VMD快速地繪制靜電勢著色的分子范德華表面圖和分子間穿透圖》一文里的ESPiso.bat進行計算,由于默認計算的是Low quality grid的靜電勢格點數據,在-ESPrhoiso 0.001設置下僅需4秒即可算完,真是超快。


    3.2 較大分子:瑞德西韋

    下面再看更大的體系,瑞德西韋,共77個原子,波函數在B3LYP/6-31G*下產生,其表面靜電勢圖如下。用的幾何結構是筆者在《使用Molclus結合xtb做的動力學模擬對瑞德西韋(Remdesivir)做構象搜索》(http://bbs.keinsci.com/thread-16255-1-1.html)中做構象搜索得到的最穩定結構。

    這里測試的是High quality grid設置下的靜電勢格點數據計算耗時。測試在雙路E5-2696v3共36核機子上完成,用的是Linux版Multiwfn和cubegen。耗時如下
    Multiwfn內部代碼算靜電勢 -ESPrhoiso 0.001:69 秒
    Multiwfn調用cubegen算靜電勢:105 秒

    可見即便對于較大的體系、即便對于很多格點數、即便對于CPU核數較多的情況,利用本文的降低耗時的技巧,耗時也明顯低于調用cubegen。

    筆者也嘗試了用i7-10870H(8核),直接用《使用Multiwfn+VMD快速地繪制靜電勢著色的分子范德華表面圖和分子間穿透圖》一文里的ESPiso.bat對瑞德西韋進行計算,涉及到計算Medium quality grid的電子密度格點數據和Low quality grid的靜電勢格點數據,在-ESPrhoiso 0.001設置下總共僅花了55秒就算完了。也就是說,在普通個人計算機上獲得上面那張靜電勢圖的計算代價還不到一分鐘!這可謂相當理想了。


    4 總結

    本文介紹了一種特別重要的在Multiwfn結合VMD繪制靜電勢著色的分子表面圖的過程中巨幅降低靜電勢計算耗時的方法。對于一般體系,用這種方法比讓Multiwfn調用cubegen來算靜電勢格點數據還要耗時低得多,因此裝有Gaussian的用戶也沒必要再讓Multiwfn調用cubegen了,而對于沒有Gaussian的用戶,此文的做法更是帶來巨大的福音。而且相對于調用cubegen,直接用Multiwfn自己的代碼計算靜電勢還有一個好處就是可以監控計算進度,而調用cubegen的時候則完全被蒙在鼓里。

    值得一提的是,根據筆者對其它體系的測試發現,對于目前的Multiwfn來說,體系越大,本文做法相對于調用cubegen的耗時優勢越小。對于非常大的體系,比如150~200個原子及以上的,調用cubegen則大概率會快于利用本文介紹的方法。所以如果你要算非常大體系,比如《巨大體系的范德華表面靜電勢圖的快速繪制方法》(http://www.shanxitv.org/481)里所示例的體系,如果有cubegen能調用的話還是建議調用cubegen。

    久久精品国产99久久香蕉