化學資訊學小技巧：Psi4的HOMO、LUMO計算與軌域可視化

Psi4 是一套開源的量子化學程式，可以高效、高精度地模擬各種分子特性。我們在之前的幾篇文章中也介紹過了一些用法。

[embed]化學資訊學：Psi4 — Python上的量子化學計算 量子化學計算在化學資訊學與計算化學上是一個非常重要的工具。我們會簡單介紹Psi4的使用方法。medium.com

HOMO和LUMO分別指最高占據分子軌域（Highest Occupied Molecular Orbital）和最低未占分子軌域（Lowest Unoccupied Molecular Orbital）是相當重要的物理量之一。

HOMO 和 LUMO 之間的能量差是 HOMO-LUMO gap，其大小可用於預測過渡金屬配合物的強度和穩定性，以及它們在溶液中產生的顏色。HOMO-LUMO gap 通常用於估計藥物引起的毒性風險。這次就來介紹下面是利用 Psi4 計算 HOMO、LUMO 的方法，配合 Avogadro 來進行軌域可視化。

Psi4 的 HOMO、LUMO 計算

這次我們用螢光素（fluorescein）這個分子來進行 HOMO、LUMO 計算。

mol = Chem.MolFromSmiles("OC(=O)C1=CC=CC=C1C1=C2C=CC(=O)C=C2OC2=CC(O)=CC=C12")
mol

然後是設定 Psi4 的計算條件和準備 xyz 格式：

psi4.set_num_threads(2)
psi4.set_memory('4 GB')

AllChem.EmbedMolecule(mol, useExpTorsionAnglePrefs=True,useBasicKnowledge=True)
AllChem.UFFOptimizeMolecule(mol)
atoms = mol.GetAtoms()
string = ""
for i, atom in enumerate(atoms):
    pos = mol.GetConformer().GetAtomPosition(atom.GetIdx())
    string += "{} {} {} {}\n".format(atom.GetSymbol(), pos.x, pos.y, pos.z)

下面是利用 Psi4 計算 HOMO、LUMO 的方法，這裡將 return_wfn 參數設為 True 來取得波函數資訊，然後從波函數資訊裡把 HOMO、LUMO 的資訊出出來：

molecule = psi4.geometry(string)
scf_e, scf_wfn = psi4.energy("B3LYP/6-31G(d)", molecule=molecule,return_wfn=True)
HOMO = scf_wfn.epsilon_a_subset('AO', 'ALL').np[scf_wfn.nalpha()]
LUMO = scf_wfn.epsilon_a_subset('AO', 'ALL').np[scf_wfn.nalpha() + 1]
psi4.fchk(scf_wfn, "molecule.fchk")
print(HOMO, LUMO, scf_e)

-0.09215683485748985 -0.07077156194021743 -1145.427919209812

scf_wfn.nalpha() 是找出可以找出 HOMO 的位置，scf_wfn.nalpha()+1 就自然是 LUMO 了。

這裡使用 psi4.fchk(scf_wfn, “molecule.fchk”) 把波函數資訊也儲存下來，會在下一步使用 Avogadro 進行軌域可視化

用 Avogadro 進行軌域可視化

Avogadro 是一款分子編輯器和視覺化工具，專為計算化學、分子建模、生物資訊學、材料科學和相關領域的跨平台使用而設計，它提供靈活的高品質渲染和強大的插件架構。

用 Avogadro 開啟剛剛儲存的波函數資訊 (molecule.fchk) ，就可以得到軌域的能量和可視化結果。

結語

利用 Psi4 和 Avogadro 可以輕鬆且免費完成基本量子化學計算和可視化，Psi4 還可以進行更多不一樣的計算，唯一要注意的是計算資源，如果計算的分子太大，可以能需要更好的 CPU 才能跑完計算。