💨

RDkit+Psi4で量子計算+可視化

2024/02/15に公開
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