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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
|
/* -------------------------------------------------------------------------- *
* SimTK Core: SimTK Simbody(tm) *
* -------------------------------------------------------------------------- *
* This is part of the SimTK Core biosimulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2006-7 Stanford University and the Authors. *
* Authors: Christopher Bruns *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "SimTKmolmodel.h"
#include "molmodel/internal/Pdb.h"
#include <sstream>
using namespace SimTK;
using namespace std;
#define ASSERT(cond) {SimTK_ASSERT_ALWAYS(cond, "Assertion failed");}
void testInputMatchesOutput() {
// A very small protein
String inputPdb = ""\
"ATOM 00001 N HIS A 19 28.165 29.227 23.618 1.00 91.78 N\n"
"ATOM 00002 CA HIS A 19 27.004 29.173 22.731 1.00 91.74 C\n"
"ATOM 00003 C HIS A 19 26.321 27.818 22.666 1.00 81.05 C\n"
"ATOM 00004 O HIS A 19 25.105 27.739 22.805 1.00 89.56 O\n"
"ATOM 00005 CB HIS A 19 27.248 29.629 21.270 1.00 98.21 C\n"
"ATOM 00006 CG HIS A 19 27.954 30.936 21.082 1.00100.00 C\n"
"ATOM 00007 ND1 HIS A 19 28.852 31.112 20.015 1.00100.00 N\n"
"ATOM 00008 CD2 HIS A 19 27.882 32.111 21.798 1.00100.00 C\n"
"ATOM 00009 CE1 HIS A 19 29.310 32.368 20.116 1.00100.00 C\n"
"ATOM 00010 NE2 HIS A 19 28.753 32.997 21.176 1.00100.00 N\n"
"ATOM 00011 N SER A 20 27.057 26.778 22.303 1.00 58.48 N\n"
"ATOM 00012 CA SER A 20 26.412 25.501 22.190 1.00 45.35 C\n"
"ATOM 00013 C SER A 20 26.089 24.909 23.558 1.00 30.40 C\n"
"ATOM 00014 O SER A 20 26.808 25.067 24.542 1.00 29.60 O\n"
"ATOM 00015 CB SER A 20 27.206 24.513 21.344 1.00 49.69 C\n"
"ATOM 00016 OG SER A 20 26.466 23.310 21.049 1.00 22.13 O\n"
"TER\n"
"END\n";
istringstream inStream(inputPdb);
PdbStructure pdbStructure(inStream, PdbStructure::InputType::PDB);
ostringstream outStream;
pdbStructure.write(outStream);
cout << inputPdb << "\n";
cout << outStream.str();
ASSERT(inputPdb == outStream.str());
}
// SimTK core bug 653
// Avoid distorting input geometry
void testMatchDefaultBreaksPlanarity() {
// A very small RNA
String pdbString = ""\
"ATOM 1 O5' A A 1 0.000 0.000 0.000 1.00 0.00 O\n"
"ATOM 2 C5' A A 1 -0.731 1.221 0.000 1.00 0.00 C\n"
"ATOM 3 H5' A A 1 -0.380 1.863 -0.837 1.00 0.00 H\n"
"ATOM 4 H5'' A A 1 -0.570 1.749 0.965 1.00 0.00 H\n"
"ATOM 5 C4' A A 1 -2.202 0.927 -0.173 1.00 0.00 C\n"
"ATOM 6 H4' A A 1 -2.775 1.880 -0.152 1.00 0.00 H\n"
"ATOM 7 O4' A A 1 -2.416 0.262 -1.447 1.00 0.00 O\n"
"ATOM 8 C3' A A 1 -2.649 -0.071 0.888 1.00 0.00 C\n"
"ATOM 9 H3' A A 1 -1.861 -0.843 1.026 1.00 0.00 H\n"
"ATOM 10 O3' A A 1 -2.867 0.610 2.118 1.00 0.00 O\n"
"ATOM 11 C2' A A 1 -3.982 -0.569 0.339 1.00 0.00 C\n"
"ATOM 12 H2' A A 1 -4.164 -1.607 0.693 1.00 0.00 H\n"
"ATOM 13 C1' A A 1 -3.707 -0.583 -1.164 1.00 0.00 C\n"
"ATOM 14 H1' A A 1 -4.620 -0.263 -1.711 1.00 0.00 H\n"
"ATOM 15 O2' A A 1 -4.991 0.378 0.624 1.00 0.00 O\n"
"ATOM 16 HO2' A A 1 -5.049 0.485 1.576 1.00 0.00 H\n"
"ATOM 17 N9 A A 1 -3.416 -1.965 -1.639 1.00 0.00 N\n"
"ATOM 18 C8 A A 1 -2.272 -2.678 -1.387 1.00 0.00 C\n"
"ATOM 19 H8 A A 1 -1.653 -2.510 -0.518 1.00 0.00 H\n"
"ATOM 20 N7 A A 1 -2.001 -3.578 -2.293 1.00 0.00 N\n"
"ATOM 21 C5 A A 1 -3.039 -3.450 -3.209 1.00 0.00 C\n"
"ATOM 22 C4 A A 1 -3.909 -2.469 -2.818 1.00 0.00 C\n"
"ATOM 23 C6 A A 1 -3.333 -4.123 -4.406 1.00 0.00 C\n"
"ATOM 24 N3 A A 1 -5.033 -2.063 -3.454 1.00 0.00 N\n"
"ATOM 25 C2 A A 1 -5.370 -2.657 -4.589 1.00 0.00 C\n"
"ATOM 26 N1 A A 1 -4.448 -3.699 -5.015 1.00 0.00 N\n"
"ATOM 1 H2 A A 1 -6.262 -2.352 -5.114 1.00 0.00 H\n"
"ATOM 28 N6 A A 1 -2.564 -5.104 -4.896 1.00 0.00 N\n"
"ATOM 29 H61 A A 1 -1.732 -5.389 -4.399 1.00 0.00 H\n"
"ATOM 30 H62 A A 1 -2.817 -5.558 -5.761 1.00 0.00 H\n"
"ATOM 31 P A A 1 1.610 0.000 0.000 1.00 0.00 P\n"
"ATOM 32 OP1 A A 1 2.103 0.698 1.208 1.00 0.00 O\n"
"ATOM 33 OP2 A A 1 2.103 0.698 -1.208 1.00 0.00 O\n"
"ATOM 34 OP3 A A 1 2.103 -1.395 0.000 1.00 0.00 O\n";
istringstream pdbStream(pdbString);
PdbStructure pdbStructure(pdbStream, PdbStructure::InputType::PDB);
// 1) First reproduce undesired behavior
RNA mol1("A");
mol1.setPdbChainId("A");
Compound::AtomTargetLocations atomTargets = mol1.createAtomTargets(pdbStructure);
double myMaxObservedSinePlaneDeviation;
mol1.matchDefaultAtomChirality(atomTargets, myMaxObservedSinePlaneDeviation);
mol1.matchDefaultBondLengths(atomTargets);
mol1.matchDefaultBondAngles(atomTargets);
mol1.matchDefaultDihedralAngles(atomTargets, Compound::DistortPlanarBonds);
mol1.matchDefaultTopLevelTransform(atomTargets);
// Most distorted part in problem report is distance between atoms N1 and C2
const ResidueInfo& res1 = mol1.getResidue( ResidueInfo::Index(0) );
Vec3 atomN1Pos = mol1.calcDefaultAtomFrameInCompoundFrame(res1.getAtomIndex("N1")).p();
Vec3 atomC2Pos = mol1.calcDefaultAtomFrameInCompoundFrame(res1.getAtomIndex("C2")).p();
Real bondLength = (atomN1Pos - atomC2Pos).norm();
ASSERT(bondLength > 0.20); // distorted
// 2) Repair with extra parameter on matchDefaultAtomChirality
RNA mol2("A");
mol2.setPdbChainId("A");
atomTargets = mol2.createAtomTargets(pdbStructure);
mol2.matchDefaultAtomChirality(atomTargets,myMaxObservedSinePlaneDeviation, 0.20);
mol2.matchDefaultBondLengths(atomTargets);
mol2.matchDefaultBondAngles(atomTargets);
mol2.matchDefaultDihedralAngles(atomTargets, Compound::DistortPlanarBonds);
mol2.matchDefaultTopLevelTransform(atomTargets);
// Most distorted part in problem report is distance between atoms N1 and C2
const ResidueInfo& res2 = mol2.getResidue( ResidueInfo::Index(0) );
atomN1Pos = mol2.calcDefaultAtomFrameInCompoundFrame(res2.getAtomIndex("N1")).p();
atomC2Pos = mol2.calcDefaultAtomFrameInCompoundFrame(res2.getAtomIndex("C2")).p();
bondLength = (atomN1Pos - atomC2Pos).norm();
cout << "bondLength = " << bondLength << endl;
mol2.writeDefaultPdb(cout);
ASSERT(bondLength < 0.20); // not distorted
// 3) Repair by not setting torsion of bonded planar atoms
RNA mol3("A");
atomTargets = mol3.createAtomTargets(pdbStructure);
mol3.matchDefaultAtomChirality(atomTargets,myMaxObservedSinePlaneDeviation);
mol3.matchDefaultBondLengths(atomTargets);
mol3.matchDefaultBondAngles(atomTargets);
mol3.matchDefaultDihedralAngles(atomTargets, Compound::KeepPlanarBonds);
mol3.matchDefaultTopLevelTransform(atomTargets);
// Most distorted part in problem report is distance between atoms N1 and C2
const ResidueInfo& res3 = mol3.getResidue( ResidueInfo::Index(0) );
atomN1Pos = mol3.calcDefaultAtomFrameInCompoundFrame(res3.getAtomIndex("N1")).p();
atomC2Pos = mol3.calcDefaultAtomFrameInCompoundFrame(res3.getAtomIndex("C2")).p();
bondLength = (atomN1Pos - atomC2Pos).norm();
ASSERT(bondLength < 0.20); // not distorted
}
int main() {
try {
testInputMatchesOutput();
testMatchDefaultBreaksPlanarity();
}
catch (const std::exception& e) {
printf("EXCEPTION THROWN: %s\n", e.what());
return 1;
}
catch (...) {
printf("UNKNOWN EXCEPTION THROWN\n");
return 1;
}
cout << "Done" << endl;
return 0;
}
|
