Scanning the torsional potential in Psikit(RDKit+Psi4)

Considering the conformational effects of the compound is important in Structure Based Drug Design, this paper discussed about it, in terms Protein-Ligand binding using torsional scan of each ligands(PDB:2JH0,2JH5,2JH6). They calculated torsional energy, and explained the relation between inhibitory activity and torsional energy.

Torsional scanning is the task of the quantum chemistry rather than that of chemoinformatics. But I wanted to conduct quantum chemical calculation as an extention of chemoinformatics way, so I implemented it in Psikit

日本語訳

jupyterでRDKitからのB3LYP/6-31G*でのDFT計算さっくりうごくの素晴らしい。

Psikit as a Electronic-Structure Informatics library

Now pen and I are planning to develop a library for a Electronic-Structure Informatics named Psikit (Psi4 + RDKit).

Which coding style do you like?

The repository is here.

psikit as a wrapper library?

import psikit
import numpy as np
from rdkit import Chem

psikit.core.set_output_file("out.dat", True)
psikit.set_memory('4 GB')
psikit.set_num_threads(4)

mol = Chem.MolFromSmiles("c1ccccc1")
benz = psikit.geometry(mol)
scf_e, scf_wfn = psikit.energy("B3LYP/cc-pVDZ", return_wfn=True)
HOMO = psikit.get_homo(scf_wfn)
LUMO = psikit.get_lumo(scf_wfn)
print(HOMO, LUMO, scf_e)

or psikit as a set of util functions library?

import psi4
from psikit import mol2xyz, get_homo, get_lumo
import numpy as np
from rdkit import Chem
psi4.core.set_output_file("out.dat", True)
mol = Chem.MolFromSmiles("c1ccccc1")
xyz, mol = mol2xyz(mol)
psi4.set_memory('4 GB')
psi4.set_num_threads(4)
benz = psi4.geometry(xyz)
scf_e, scf_wfn = psi4.energy("B3LYP/cc-pVDZ", return_wfn=True)
HOMO = get_homo(scf_wfn)
LUMO = get_lumo(scf_wfn)
print(HOMO, LUMO, scf_e)

RDKitのMolオブジェクトにarrayデータを属性として持たせる方法(強引に)

RDKitのMolオブジェクトにリストとかnumpyのアレイをプロパティとしてセットしたかったので、Gregさんにメソッド用意されていないのか、やり方がないのか聞いてみたら、C++じゃないと出来ないとのこと。

セット系だと

  • SetBoolProp
  • SetDoubleProp
  • SetIntProp
  • SetProp
  • SetUnsignedProp

しか用意されていなくて、困ったなぁと途方にくれていたところ、 SetPropって文字列格納するからとりあえずシリアライズすればいけるんじゃないの ということでやってみた

>>> from rdkit import Chem
>>> import numpy as np
>>> import pickle
>>> mol = Chem.MolFromSmiles("CC")
>>> ar = np.array([[1,2,3],[4,5,6]])
>>> mol.SetProp("arraydata",pickle.dumps(ar))
>>> pickle.loads(mol.GetProp("arraydata"))
array([[1, 2, 3],
       [4, 5, 6]])

というわけで無事にセットできましたとさ。pickle使っているのはnumpyだとファイルかファイルハンドルとるので面倒くさかったから。