1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
|
#include <openbabel/mol.h>
#include <openbabel/obconversion.h>
#include <iostream>
#include <fstream>
#include <iterator>
#include <string>
#include <map>
#include <vector>
using namespace std;
using namespace OpenBabel;
int main(int argc, char *argv[])
{
// turn off slow sync with C-style output (we don't use it anyway).
std::ios::sync_with_stdio(false);
OBConversion conv;
conv.SetOptions("O", conv.OUTOPTIONS);
conv.SetOutFormat("can");
ifstream ifs;
OBMol mol;
if (argc < 2)
{
cerr << "Usage: obtorsion <file>" << endl;
return(-1);
}
map<string, vector<double> > torsion_angles;
for (int i = 1; i < argc; i++) {
cerr << " Reading file " << argv[i] << endl;
OBFormat *inFormat = conv.FormatFromExt(argv[i]);
if(inFormat==nullptr || !conv.SetInFormat(inFormat))
{
cerr << " Cannot read file format for " << argv[i] << endl;
continue; // try next file
}
ifs.open(argv[i]);
if (!ifs)
{
cerr << "Cannot read input file: " << argv[i] << endl;
continue;
}
while(ifs.peek() != EOF && ifs.good())
{
conv.Read(&mol, &ifs);
mol.DeleteHydrogens();
FOR_BONDS_OF_MOL(bond, mol) {
if(bond->IsRotor()) {
OBBitVec atomsToCopy;
OBAtom *atom = bond->GetBeginAtom();
FOR_NBORS_OF_ATOM(a, &*atom) {
atomsToCopy.SetBitOn(a->GetIdx());
}
atom = bond->GetEndAtom();
FOR_NBORS_OF_ATOM(a, &*atom) {
atomsToCopy.SetBitOn(a->GetIdx());
}
OBMol mol_copy;
mol.CopySubstructure(mol_copy, &atomsToCopy);
string smiles = conv.WriteString(&mol_copy, true);
OBAtom* b = bond->GetBeginAtom();
OBAtom* c = bond->GetEndAtom();
OBAtom* a = nullptr;
FOR_NBORS_OF_ATOM(t, &*b) {
a = &*t;
if(a != c)
break;
}
OBAtom* d = nullptr;
FOR_NBORS_OF_ATOM(t, &*c) {
d = &*t;
if(d != b)
break;
}
double angle = mol.GetTorsion(a, b, c, d);
if(torsion_angles.count(smiles) == 0) {
vector<double> t(36, 0);
torsion_angles[smiles] = t;
}
torsion_angles[smiles][static_cast<size_t>(angle+180)/10] += 1;
}
}
}
}
map<string, vector<double> >::iterator it;
for(it=torsion_angles.begin(); it!=torsion_angles.end(); it++) {
vector<double>::iterator j;
int max_index = distance(it->second.begin(), max_element(it->second.begin(), it->second.end()));
cout << it->first << " " << max_index * 10 - 180 << endl;
}
}
|
