文章簡介
作者采用PM6方法計算了分子幾何優化構型和電子性質，計算基數為PubChem Compounds（截止2016年8月29日的數據）的9,160萬個分子的94%分子。除了中性狀態的分子，作者還計算了56.2%的陽離子狀態、49.7%的陰離子狀態和41.3%的自旋翻轉狀態的分子。因此，作者采用PM6方法計算的分子總數達2.21億個。作者對比了260萬個分子的PM6法和B3LYP/6-31G*法的幾何優化構型結果，鍵長和鍵角的標準差分別為0.016 ?和1.7°；并對兩種方法的HOMO能級和LUMO能級進行線性回歸分析，得出HOMO能級關系為EB3LYP(HOMO) = 0.876EPM6(HOMO) + 1.975（eV），計算確定的系數為0.803；LUMO能級關系為EB3LYP(LUMO) = 1.069EPM6(LUMO) – 0.42（eV），計算確定的系數為0.842。作者還生成了4個子數據集，每個子集中的分子質量都低于500，子集一包含 C、H、O、N四個元素，子集二包含C、H、N、O、P、S六個元素；子集三包含C、H、N、O、P、S、F、Cl八個元素；子集四包含C、H、N、O、P、S、F、Cl、Na、K、Mg、Ca十二個元素。
背景介紹
隨著有機薄膜太陽能電池、電發光材料、有機非線性光學材料、分子傳感器和創新藥的開發和發展，新的有機分子的發現越來越重要。太陽能電池的化學性質如激發能、最高占據分子軌道（HOMO）與最低未占據分子軌道（LUMO）的能隙等非常重要，而量子化學從頭計算法不需要昂貴的物理或化學實驗就可以精確預測化學性質。遺憾的是，這種量子化學計算方法的計算費用雖然已經降低很多，但對于研究分子化合物的化學空間來說仍然太慢。根據Lipinski成藥性規則，目前候選藥物預計已達到1060個，數量相當龐大。機器學習被證實能彌補量子化學計算方法缺陷，該方法在過去十年已經廣泛應用于分子和固體研究，尤其是深度學習，已經引起相當大的關注。深度學習基于定量結構-性質關系（QSPR）數據訓練，采用多層人工神經網絡創建預測模型，而不是采取“第一性原理”計算化學性質?？偟膩碚f，為了使用有監督深度學習方法來開發精準模型，我們需要大量的可信賴數據用于訓練。
我們報告了大規模的量子化學數據集——PubChemQC PM6數據集，數據集采用PM6方法對PubChem Compounds（提取日期2016年8月29日）上的分子進行幾何結構最優化、電子結構和其他性質的計算。PubChem Compounds是最全面的化合物數據庫之一，截止2016年8月29日已收錄了大約9,160萬個分子，且每天都在增加。在撰寫本文時，作者已經成功完成計算的分子包括86,213,135個中性分子、51,555,911個陽離子態分子、45,581,750個陰離子態分子及37,839,619個自旋翻轉態分子，合計計算了221,190,415個分子的最優幾何結構和電子結構。PubChemQC PM6覆蓋了超過94%的中性分子。此外，作者為了方便用戶使用，通過限制元素和分子量、去除鹽基等維度創建了4個子數據集。數據集網址：http://pubchemqc.riken.jp/pm6_datasets.html。
建立這些數據集的目的是向研究人員或從業人員提供全面的量子化學數據集，應用領域包括但不限于機器學習、高通量虛擬篩選等。目前互聯網上已有多個量子化學數據集，但它們的范圍和規模僅限于幾百上千個有機小分子。據目前所知，在半經驗量子化學計算方法創建的數據集中，PubChemQC PM6數據集是最大的。如果忽略計算方法及其計算準確性的不同，PubChemQC PM6數據集里的記錄數量大大超過了任何通過量子化學計算獲得的數據集。采用密度泛函理論（DFT）方法創建的數據集比PM6半經驗方法的準確性高，如哈佛清潔能源項目數據庫（Harvard Clean Energy Project Database）和ANI-1數據集，前者收錄了230萬個由26個構建基塊組成的有機光伏候選化合物，后者收錄了57,462個有機小分子的2,000萬個經計算的非平衡構象。然而，和PubChemQC PM6相比，這兩種數據庫的分子種類非常有限，PM6數據集包含了2.21億個分子幾何和電子結構，覆蓋70種元素。作者使用PM6方法建立數據集的原因有3個，第一，計算成本大大低于那些更準確的計算方法；第二，作為分子結構優化的初步推測，計算結果更優越，因為采用PM6方法的分子幾何優化能用于構象分析，能用于減少各種高精密計算方法的計算成本（如B3LYP/6-31G*, CCSD(T)/ccpVDZ, 等等）。值得一提的是，PM6方法為包括蛋白質在內的分子提供了相當優異的幾何結構；第三，計算結果本身就具有價值，其能用于預測分子的電子結構。據文獻報告，PM7方法結合機器學習可能可以高準確性的預測電子結構，值得一提的是PM6和PM7方法在小分子的幾何優化方面有著相似的準確性。為了確認計算出的分子幾何結構的準確性，作者對比了260萬個PM6法和B3LYP/6-31G*法優化的幾何結構。B3LYP/6-31G*法計算的鍵長和鍵角的標準差分別約為0.016? and 1.7°，這些值幾乎和Stewart的原始結果一樣。
創建PubChemQC PM6數據集的方法
——選擇分子的標準
分子的數量極大，即使限制原子種類和組成的原子數量，分子的數量仍然是個天文數字，例如Ruddigkeit等人創建的數據集GDB-17，里面列舉了1660億個由C、N、O、S和鹵素等17種原子組成的有機小分子，這個數據集的問題就在于不能直接判斷哪個分子必不可少。Ramakrishnan等人建立了GDB-17的子集，包括QM7、QM7b、QM8和QM9數據集，里面采用量子化學方法計算了幾何和電子結構。其中QM9數據集是最大的，但里面只有13.4萬個分子。此外，當分子大小稍微大一點點，采用枚舉法的GDB-17數據集可能就不能正常運行，比如C32的異構體已經能列出27,711,253,769個，盡管這些異構體間的區別幾乎可以說是微不足道。又如順反異構，N個雙鍵的順反異構體有2N個，但是有些異構體不能展示出來，不能直接決定一對異構體中應該選擇哪個。有些異構體非常不同，有些異構體卻相似，如順式不飽和脂肪酸有利于維持良好的膽固醇水平，而反式脂肪酸則被認為是有害的；1, 3二氯丙烯作為農業化工品使用時，順反結構沒有明顯的區別，通?；旌嫌糜谵r藥。因此，與其主觀選擇分子，作者決定基于現有的已收錄必要分子的數據庫來開發化學數據庫，并盡可能多的對數據庫中的分子進行量子化學計算。目前已有很多的化學數據庫包括CAS、ChEMBL、ChemSpider、Zinc等等，作者選擇PubChem化合物庫是因為其可及性、大小和多樣性，Zinc、ChEMBL和PubChem可免費試用，其他數據庫有部分限制。CAS是專有數據庫，ChemSpider需要權限才能組裝超過5000個結構的數據庫。PubChem（2016年）收錄了大約1億個分子，而ChEMBL15只有180萬個分子，Zinc15有1.2億個分子。盡管PubChem和Zinc15的數據庫大小相當，但作者仍選擇了PubChem，因為PubChem數據來源多樣，收錄的信息不止來自化學品供應商目錄和期刊出版商，甚至來自于Zinc，ChEMBL和ChemSpider。PubChem是一個開放式化學數據庫，隸屬于NIH（美國國立衛生研究院）。自2004年起，PubChem就不斷更新來自全世界的分子信息，數據采集自數百個大數據源包括大學、藥企、政府機構、化學品供應商、科學文獻和其他管理機構。PubChem包含3個子數據庫：Substances, Compounds,和BioAssays，Substances收集申報數據，Compounds收錄了從Substances數據庫提取的獨特的標準化合物數據。截止2016年8月16日，PubChem Compounds收錄了91,679,247個化合物，每一個化合物都有InChI和SMILES編碼以及PubChem CID號碼，分子式和其他信息。
創建PubChemQC PM6數據集的方法
——分子的展示
一個分子沒有嚴格的定義，但是我們可以在一定的假設條件下定義一個分子。在Born-Oppenheimer近似氣相、非相對論極限和點電荷模型中，我們可以在笛卡爾坐標系中從一組原子（核電荷）中確定一個分子的漢密爾頓量以及系統中的電子數量。然后我們可以通過解薛定諤方程來獲得波函數和量子數。用這種方法，一個分子就由漢密爾頓量、波函數和量子數定義，這種定義的弊端是不易區分兩個相同的分子或者在其他方面不同的分子如同分異構體。另一方面，最方便的展示方法是有一個通用名，而每個新引入的化合物都必須命名。我們可以簡單采用IUPAC命名法作為系統命名法，因為大部分情況下原子、原子間的連接、鍵級和其他立體信息等就足以確定一個分子。但是我們不容易在電腦運行IUPAC名字?；诖俗髡邲Q定采用人類和機器可讀的分子編碼系統如InChI和SMILES，能與IUPAC命名法兼容，采用類似化學分子式的系統方法來編碼分子。PubChem Compounds采用InChI和SMILES兩種方法來編碼化合物，作者在計算過程中廣泛使用了這些編碼，即使用SMILES系統生成初始幾何結構預測，用InChI系統確認最優結果。需要注意的是這些編碼系統也有一些模糊的地方，只能確定一組原子、電荷和原子連接的一些信息以及立體異構。
PubChemQC PM6數據集的詳細信息
——創建數據集
首先，從PubChem站點下載所有分子的SDF文件，解析和提取每一個分子的CID、分子量、InChI和SMILES信息、分子式、電子電荷和自旋量子數，通過異構體的SMILES編碼計算電子電荷，根據系統的電子數量的奇偶性將自旋量子數設置為0或1。然后根據分子量進行降序排列，剔除分子量大于1000g/mol的分子（上限高于Lipinski成藥性原則的500g/mol）后剩余PubChem Compound中的604,330個（占比0.66%）分子。同時剔除2,188,881個（占比2.39%）帶電分子。輸入文件由Open Babel生成，采用分子的異構SMILES編碼，選項為-addh和-gen3d。通過對比原始數據集中的分子式和分子的SMILES編碼生成的結構來確認氫原子數量，因為Open Babel有時會給分子添加過多或過少的氫原子。作者沒有用PubChem3D來展示分子幾何結構，原因在于（1）SMILES和InChI表示法比分子幾何結構小得多，因此更容易處理；（2）和PM6分子幾何優化相比，Open Babel計算初始幾何預測值所需要的計算資源是微不足道的；（3）PubChem3D的分子幾何構型的可及性受限。最終，作者使用Gaussian09來進行每個分子的PM6幾何構型的優化計算。成功之后，采用PM6幾何優化模型對陽離子態、陰離子態和自旋翻轉狀態的分子也進行了幾何構型最優化，得到初始預測值。接下來使用Open Babel對分子的中性狀態進行確認，確認計算好的最優幾何構型的InChI是否與原始InChI一致。采用sed腳本來描述2個InChI表示法，并驗證輸出的特征。作者忽略了游離氫和原子的形式電荷，因為不能預測或確定它們的量子化學計算輸出的意義。作者還忽略了總電荷，總電荷在幾何優化過程中應不變。同樣的，立體異構體、幾何異構體和構象異構體也不在考慮范圍內。每一個CID數據如下：輸入文件，xyz格式的原子坐標和JSON文件。所有計算都在RIKEN HOKUSAI BigWave超級計算機（Intel Xeon Gold 6148 2.4GHz, 1680CPUs, 33600 cores），QUEST集群 (Intel Core2 L7400 1.50 GHz, 700 nodes, 1400 cores), 和 RIKEN RICC 超級計算機 (Intel Xeon 5570 2.93 GHz, 1024 nodes, 8192 cores)運行?？偟挠嬎銜r間為HOKUSAI BigWave 95天、QUEST 346 天、RICC 126天，2016年12月30日開始計算，到2018年6月9日結束，2018年10月12日完成項目總結。
PubChemQC PM6數據集的詳細信息
——數據
表1. PubChemQC PM6的數據量，總共計算了86,213,135個中性分子，其中85,197,307個分子的原始InChI和PM6優化幾何結構計算得出的InChI的化學分子式和主層原子連接相一致。
表格源自JCIM
?
表2. PubChemQC PM6數據集中每個分子可用的數據，包括振動頻率、振動強度、SCF能量。
表格源自JCIM
?
表3. PubChemQC PM6子數據集的數據，每個子集包含的分子質量都低于500，所有分子都是中性狀態和單峰狀態，CHON意思是化合物只包含C, H, N, O原子，前三個子集均去除了鹽基。
表格源自JCIM
表1和2列出了PM6幾何結構優化數據，包括大量化合物和包含詳細數據的文件名字，比如PubChem CID、分子量、InChI、異構體的SMILES以及分子式。比如，“MW less than 1000”列出了PubChem Compound庫中分子量低于1000的分子數量；“charged molecules”列出了帶電分子的數量和包含詳細數據的文件名；“no results”列出了PM6幾何結構優化失敗的分子數量；“InChI (in)valid”列出了在PM6優化幾何構象中原始InChI和經計算過的InChI的化學分子式和主層原子連接(不)一致的分子數量；“cations” “anions” “spin flipped”是計算的陽離子、陰離子和自旋翻轉狀態的分子，這些分子都是來自PM6優化過的幾何中性狀態分子，較中性分子不穩定，因此更難計算?！癵rand total”是中性狀態、陽離子狀態、陰離子狀態和自旋翻轉狀態的分子總和。為了用戶方便，作者創建了4個較小的數據集（表3），選擇分子質量低于500的分子，符合Lipinski原則之一。第一個子集包含了生物體中最常見的元素，第二個子集包含了生物分子的重要元素，第三個子集加入了輕鹵素原子F和Cl，第四個子集包含了人機體的基礎元素，除了氟。
PM6和B3LYP分子計算的對比
——分子的幾何構型
為了看看PM6和B3LYP/6-31G*優化幾何結構計算方法的區別，作者選取鍵長和鍵角的標準差（RMSD）進行比較，這兩個指標能獲得相對準確的實驗結果，因此常用于評估量子化學方法的準確性。另一方面，比較分子的標準差較難，不做分子的標準差對比的原因如下：第一，實驗已表明乙烷在室溫條件下容易發生異構翻轉，烯醇的互變異構化也頻繁發生，因此較難獲得精確的結果；第二，構象異構體的搜索從計算方面來講具有挑戰性，即使是如L-半胱氨酸這樣的小分子，也需要耗費2個月的時間來獲得穩定的構象。因此，作者只對比PM6和B3LYP/6-31G* 這兩種幾何優化方法所預測的2,606,946個分子的鍵長和鍵角。加入對比的分子數量少于原始PubChemQC數據集（包含3,981,230個分子），原因如下：第一，基礎數據集不同，B3LYP/6-31G*方法用的是截止2014年7月PubChem compounds里的數據，而PM6用的是2016年8月的數據，兩個數據集有很大的不同。雖然兩個數據集共用了2,777,085個分子（CIDs），與2014年7月的版本相比，2016年8月的版本去掉了超過一百萬個分子或移到不同的CIDs。第二，作者優化過程中去掉了每個CID對應的變形分子，確切的說，剔除了PM6和B3LYP優化方法計算后SMILES不一致的分子。?
圖1. PM6和B3LYP/6-31G*方法對2,606,946個分子的幾何優化構象的鍵長標準差。
圖片源自JCIM
圖2. PM6和B3LYP/6-31G*方法對2,606,946個分子的幾何優化構象的鍵角標準差。
圖片源自JCIM
圖1描繪了兩種方法分子優化后的鍵長（單位：埃，?）標準差的直方圖。鍵長標準差直方圖的一個箱寬度（bin size）為0.001 ?，眾數是0.016 ?，尖峰明顯，第85百分位數為0.037 ?。結果顯示PM6分子幾何構象略微好于Stewart的方法，PM6計算的由C、H、O、N組成的化合物的平均不知名錯誤為0.025 ?。圖2描繪了兩種方法分子優化后的鍵角（單位：埃，?）標準差的直方圖。鍵角標準差直方圖的箱寬度為0.1°，眾數為1.70°，尖峰明顯，第85百分位數為4.3°，結果主要與原始論文結果一致，100個由C、H、O、N組成的化合物中鍵角計算的平均不知名錯誤為3.1°。作者發現PM6和B3LYP/6-31G*兩種方法計算的鍵長和鍵角的標準差以在0.016 ?和1.70°出現尖峰為特征，這些數值幾乎與PM6的原始論文報告相同，同時作者也在260萬個分子中驗證過除氫鍵外的所有鍵長和鍵角。有趣的是，所有的計算結果和B3LYP幾何結構優化結果一致，鍵角和鍵長與實驗結果相差也只在幾度和0.02 ?范圍內。
HOMO、LUMO的能級和HOMO-LUMO的能隙
HOMO能級，LUMO能級和HOMO-LUMO能隙是很重要的信息，與反應性、光激發和電荷運輸有關，對于開發有機發光二極管、有機光伏設備、有機薄膜晶體管等材料很重要。分子的HOMO-LUMO能隙越大，分子越穩定。當一個分子在HOMO能級釋放電子，另一個分子在LUMO能級接受電子時發生電荷運輸。光激發發生在一個分子吸收光子產生電子-空穴對，電子靠近HOMO能級，而空穴靠近LUMO能級。因此，比較PM6和幾何優化法和B3LYP/6-31G*幾何優化法獲得的分子的HOMO、LUMO的能級和HOMO-LUMO的能隙會很有意思。根據分析，作者預估兩者的HOMO能級，LUMO能級和HOMO-LUMO能隙關系分別如下式（圖3-5）：EB3LYP(HOMO) = 0.876EPM6(HOMO) + 1.975；EB3LYP(LUMO) = 1.069EPM6(LUMO) – 0.42和EB3LYP(GAP) = 0.959EPM6(GAP) ? 3.165。同時三個關系式對應的系數分別為0.803, 0.842和0.779。
圖3. PM6和B3LYP/6-31G*優化的HOMO能級和線性回歸分析結果。
圖片源自JCIM
圖4. PM6和B3LYP/6-31G*優化的LUMO能級和線性回歸分析結果。圖片源自JCIM
圖5. PM6和B3LYP/6031G*優化的HOMO-LUMO能隙和線性回歸分析結果。
圖片源自JCIM
前述確定的系數大部分與Pereira等人的研究成果（Pereira, F.; Xiao, K.; Latino, D. A. R. S.; Wu, C.; Zhang, Q.; Aires-de Sousa, J. Machine Learning Methods to Predict Density Functional Theory B3LYP Energies of HOMO and LUMO Orbitals. J. Chem. Inf. Model. 2017, 57, 11?21.）一致，他們使用JChem軟件以經驗為主地生成單一構象異構體，采用半經驗PM7法優化幾何構型，再采用B3LYP/6-31G*理論計算單點；他們確定的HOMO、LUMO、HOMO-LUMO能隙的系數分別為0.799、0.895、0.656。作者的方法和Pereira等人的方法稍微不同，和他們直接對比是不可能的，第一個是計算方法的不同，Pereira等人用的是PM7，作者用的是PM6；第二是分子幾何，Pereira等人用PM7優化幾何適用于B3LYP/6-31G*計算，作者采用了2種不同的幾何構型——PM7優化幾何和B3LYP/6-31G*優化幾何；第三用于對比的分子數量也不同，Pereira等人用了111,725個分子，而作者使用了2,606,946個分子，前兩點可能對結果影響不大，但是使用相當大數量的分子作為對比使作者的結果更準確。因此可以預測B3LYP/6-31G*法的HOMO能級、LUMO能級和HOMO-LUMO能隙與Pereira等人機器學習的結果相似。
展望
在PubChemQC PM6數據集的幫助下，許多研究方向都可以開展，作者列出了主要的四個方向。（1）更詳盡的分析，基于PubChemQC PM6的分子幾何優化和電子性質計算結果，可以分析更多詳細的分子性質，作者準備發表更多關于最優振動強度、模式、偶極瞬間、結構轉化等分析結果，這對于材料開發應該有用。（2）使用更精細的方法進行計算。雖然PM6是一種比較完善的半經驗量子化學計算方法，但增強計算資源仍有更多精細的方法能更準確的優化幾何結構。事實上作者基于PubChemQC PM6優化的幾何構型，正在試圖結合B3LYP方法開發計算方法，使用目前的結果創建更全面的數據集。眾所周知色散相互作用對于獲得分子的平衡幾何結構很重要，但PM6方法和B3LYP方法均不能恰當的處理色散相互作用。因此能恰當處理色散相互作用（如PM6-D3H4, DFT-D, ωB97X-D, ωB97M-V, vdW-DF-04, VV10, MP2, CCSD等）的計算方法將是另一個非常有價值的研究。因此相信PubChemQC PM6的結果能減少計算費用。（3）氣相外的計算。目前所有的PubChemQC PM6的計算結果都是用于氣相環境中，但是分子的幾何結構在溶液和晶體環境中與在氣相中有時會有很大的不同。因此在氣相環境外的計算和創建相對應的數據集將非常有用。溶液（如水、乙醇等）環境的計算雖然會較昂貴，但PubChemQC PM6能在一定程度上降低計算負擔。另外，如果沒有事先了解晶體的空間群或晶胞內的分子數，晶體環境中的計算會更加有挑戰性，因為會產生組合爆炸。（4）比較PubChemQC PM6的計算結果和實驗確定的結構。量子化學計算有時會偏離實驗數據，因此，系統對比PubChemQC PM6的計算結果和實驗結果有助于驗證和提高量子化學計算理論和實踐。然而，據我們所知，目前還沒有全面的實驗數據庫公開發表。
參考文獻
Nakata Maho, Shimazaki Tomomi, Hashimoto Masatomo and Maeda Toshiyuki, PubChemQC PM6: Data Sets of 221 Million Molecules with Optimized Molecular Geometries and Electronic Properties. J. Chem. Inf. Model. 2020, ASAP.