Created
July 3, 2023 03:32
-
-
Save TVect/9ea7d6b4fbd8c1fef0a09423c2465976 to your computer and use it in GitHub Desktop.
使用 openbabel 做固定部分结构的分子构象搜索. 来源于 https://www.insilico.jp/blog/2021/11/14/partial_conformation_search/
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #include <fstream> | |
| #include <iostream> | |
| #include <string> | |
| #include <openbabel/mol.h> | |
| #include <openbabel/obconversion.h> | |
| #include <openbabel/forcefield.h> | |
| #include <openbabel/parsmart.h> | |
| #include <openbabel/generic.h> | |
| using namespace std; | |
| using namespace OpenBabel; | |
| int main(int argc, char* argv[]) { | |
| // 力場の設定 | |
| string ff = "MMFF94"; | |
| OBForceField* pFF = OBForceField::FindForceField(ff); | |
| if (!pFF) { | |
| cerr << "指定の力場が見つかりませんでした" << endl; | |
| return -1; | |
| } | |
| pFF->SetLogFile(&cout); | |
| pFF->SetLogLevel(OBFF_LOGLVL_NONE); | |
| // ファイルオープン | |
| ifstream ifs("Donepezil.sdf"); | |
| if (ifs.fail()) { | |
| cerr << "ファイルが見つかりませんでした" << endl; | |
| return -1; | |
| } | |
| // 参照分子の読み込み | |
| OBConversion conv; | |
| conv.SetInAndOutFormats("SDF", "SDF"); | |
| OBMol mol; | |
| conv.Read(&mol, &ifs); | |
| ifs.close(); | |
| mol.DeleteHydrogens(); // 水素を外す | |
| // 母核構造のマッチング | |
| OBSmartsPattern smp; | |
| smp.Init("O=C1CCc2c1cccc2"); | |
| smp.Match(mol); | |
| vector<vector<int> > mlist = smp.GetUMapList(); | |
| // 拘束原子の設定 | |
| OBFFConstraints constraints; | |
| for (int i = 0; i < mlist[0].size(); i++) { | |
| cout << mlist[0][i] << std::endl; | |
| constraints.AddAtomConstraint(mlist[0][i]); | |
| } | |
| // 分子をセットアップ（部分電荷などの計算） | |
| if (!pFF->Setup(mol, constraints)) { | |
| cerr << "分子の設定に失敗しました" << endl; | |
| return -1; | |
| } | |
| // ランダムにコンホメーション発生 | |
| pFF->RandomRotorSearch(10, 0); | |
| pFF->GetConformers(mol); | |
| mol.AddHydrogens(); // 水素を付加 | |
| char s[99]; | |
| ofstream ofs("out.sdf"); | |
| for (int i = 0; i < mol.NumConformers(); i++) { | |
| mol.SetConformer(i); | |
| if (!pFF->Setup(mol, constraints)) { | |
| cerr << "分子の設定に失敗しました" << endl; | |
| return -1; | |
| } | |
| pFF->ConjugateGradients(5000); | |
| pFF->GetCoordinates(mol); | |
| sprintf(s, "%f", pFF->Energy()); | |
| string Unit = pFF->GetUnit(); | |
| cout << s << Unit << std::endl; | |
| auto pd = new OBPairData; | |
| pd->SetAttribute("Energy"); | |
| pd->SetValue(s); | |
| mol.SetData(pd); | |
| conv.Write(&mol, &ofs); | |
| mol.DeleteData(pd); | |
| } | |
| ofs.close(); | |
| return 0; | |
| } |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| from openbabel import openbabel, pybel | |
| pFF = openbabel.OBForceField.FindType("MMFF94") | |
| if not pFF: | |
| print("Error: force field not found") | |
| exit() | |
| pFF.SetLogToStdOut() | |
| pFF.SetLogLevel(0) | |
| mol = next(pybel.readfile("SDF", "./Donepezil.sdf")) | |
| print(mol.write("smi")) | |
| ob_mol = mol.OBMol | |
| ob_mol.DeleteHydrogens() | |
| smp = openbabel.OBSmartsPattern() | |
| smp.Init("O=C1CCc2c1cc(*)cc2") | |
| smp.Match(ob_mol) | |
| mlist = smp.GetUMapList() | |
| constraints = openbabel.OBFFConstraints() | |
| for item in mlist[0]: | |
| constraints.AddAtomConstraint(item) | |
| pFF.Setup(ob_mol, constraints) | |
| pFF.RandomRotorSearch(10, 0) | |
| pFF.GetConformers(ob_mol) | |
| ob_mol.AddHydrogens() | |
| output_file = pybel.Outputfile("sdf", "output_pcs.sdf", overwrite=True) | |
| for i in range(10): | |
| ob_mol.SetConformer(i) | |
| pFF.Setup(ob_mol, constraints) | |
| pFF.ConjugateGradients(5000) | |
| pFF.GetCoordinates(ob_mol) | |
| energy = pFF.Energy() | |
| unit = pFF.GetUnit() | |
| print(f"{i} energy: {round(energy, 4)} unit: {unit}") | |
| new_mol = pybel.Molecule(ob_mol) | |
| new_mol.data["energy"] = round(energy, 4) | |
| output_file.write(new_mol) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment