化學資訊學小技巧：用Psi4計算分子的吸收光譜

TDDFT (Time-dependent density functional theory) 是一種 time-dependent electronic many-body problem 形式的精確算法，廣泛用於分子的計算激發能。這次介紹 Psi4 的 TDDFT 計算方法與計算吸收光譜的功能。

TDDFT 計算

設定分子結構時，通常使用 XYZ 座標提供每個原子的位置資訊。將當時指定的分子結構儲存在變數中，並指定稍後要執行計算的結構。這次我們用 Prodan 這個螢光色素分子來進行 HOMO、LUMO 計算。

mol = Chem.MolFromSmiles("CCC(=O)C1=CC2=C(C=C1)C=C(C=C2)N(C)C")
mol

接下來，我們可以使用 psi4.set_XXX 設定計算環境，我們使用 2 個線程和 4GB 記憶體來計算，並準備好 xyz 格式的分子。

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

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

再來是先進行結構最佳化然後再進行能量計算。需要注意的地方是需要執行這一行指令 psi4.set_options({‘save_jk’: True})，不然在 TDDFT 計算時有可能會碰到錯誤。這裡將 return_wfn 參數設為 True 來取得波函數資訊，TDDFT 計算時會需要用到波函數。

molecule = psi4.geometry(xyz)
psi4.set_options({'save_jk': True}) 
psi4.optimize("B3LYP/6-31G(d)", molecule=molecule)
scf_e, scf_wfn = psi4.energy("B3LYP/6-31G(d)", molecule=molecule,return_wfn=True)

利用 tdscf_excitations 這個函數，Psi4 能夠以線性響應形式計算 SCF 參考波函數的激發能和基態到激發態躍遷特性。需要輸入波函數 scf_wfn 並計算的激發態數量 states = 8，以及三重激發態 triplets=’also’。8 個激發態的計算結果會儲存在 res 變數裡。

from psi4.driver.procrouting.response.scf_response import tdscf_excitations
#TD-DFT
res = tdscf_excitations(scf_wfn, states=8, triplets='also')

吸收光譜計算

吸收光譜的資訊可以從 TDDFT 計算結果 res 中得到，利用 psi4.driver.p4util.spectrum 函數可以繪製吸收光譜。輸入變數的 poles 和 residues 可利用以下程式碼從 res 中提取。

from psi4.driver.p4util import spectrum
import matplotlib.pyplot as plt

#Spectra
poles = [r["EXCITATION ENERGY"] for r in res]
opa_residues = [np.linalg.norm(r["ELECTRIC DIPOLE TRANSITION MOMENT (LEN)"])**2 for r in res]
opa_spectrum = spectrum(poles=poles, residues=opa_residues, gamma=0.01, out_units="nm")

#plot
plt.plot(opa_spectrum["convolution"]["x"], opa_spectrum["convolution"]["y"])
plt.show()

比較實際的吸收光譜

我們來比較一下計算結果和實際實驗結果的差距，由於最簡單的 TDDFT 計算是假設單分子在真空狀態下進行計算，和實際的光譜會有所差異。從文獻中我們找到了 Prodan 的吸收光譜 (excitation) 和螢光光譜 (emission)，吸收波長的區域是落在 300 nm~ 400 nm 這個區間，和 Psi4 的計算結果相近，雖然 Psi4 計算結果的吸收波長範圍稍微窄了一點。

doi: 10.1016/S0006–3495(98)77905–6

結語

這次示範了如何使用 Psi4 進行TDDFT 計算並模擬出吸收光譜，是一種 DFT 能量之後的衍生應用，我們可以輕鬆地利用 Python 上執行也可以模擬出吸收光譜。