💨
RDkit+Psi4で量子計算+可視化
def psi4_from_rad_smiles(smiles, charge, multiplicity):
#ラジカル系の計算
# RDKitとPsi4を使用して、メチルラジカルの構造最適化とエネルギー計算を行う
from rdkit import Chem
from rdkit.Chem import AllChem
import psi4
# SMILES文字列から分子を作成
smiles = smiles # 例:エタノール
mol = Chem.MolFromSmiles(smiles) # SMILESから分子を生成
# 3D構造を生成
mol = Chem.AddHs(mol) # 水素を追加
AllChem.EmbedMolecule(mol) # 3D座標を生成
#AllChem.MMFFOptimizeMolecule(mol) # 分子の最適化
AllChem.UFFOptimizeMolecule(mol) # 分子の最適化
# XYZ形式での出力
num_atoms = mol.GetNumAtoms()
xyz = f"{num_atoms}\n\n"
for atom in mol.GetAtoms():
pos = mol.GetConformer().GetAtomPosition(atom.GetIdx())
xyz += f"{atom.GetSymbol()} {pos.x} {pos.y} {pos.z}\n"
print(xyz)
# Psi4の初期設定
psi4.core.set_output_file(smiles+'_output.dat', False)
psi4.set_options({
'reference': 'uhf', # UHFを使用して開殻系を扱う
})
# 分子の定義(スピン多重度を設定)
molecule = psi4.geometry(xyz)
molecule.set_multiplicity(multiplicity) # ダブレット状態
molecule.set_molecular_charge(charge) # 電荷を設定
# 構造最適化とエネルギー計算(例:B3LYP/6-31G*)
# 波動関数とエネルギーを取得
energy, wavefunc = psi4.properties('B3LYP/6-31+G*', # 例:B3LYP/6-31G*
molecule=molecule,
return_wfn=True,
properties=['MULLIKEN_CHARGES',
])
MULLIKEN_CHARGES= psi4.oeprop(wavefunc, 'MULLIKEN_CHARGES')
mulli=wavefunc.variable('MULLIKEN_CHARGES')
print(energy, wavefunc, mulli)
## 類似度マップで可視化
from rdkit.Chem.Draw import SimilarityMaps
fig = SimilarityMaps.GetSimilarityMapFromWeights(mol, mulli, colorMap='RdBu')
return energy, wavefunc, mulli
Discussion