1 // UTILITY.CPP
2
3 // Copyright (C) 1998 Tommi Hassinen, Geoff Hutchison.
4
5 // This package is free software; you can redistribute it and/or modify
6 // it under the terms of the GNU General Public License as published by
7 // the Free Software Foundation; either version 2 of the License, or
8 // (at your option) any later version.
9
10 // This package is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU General Public License for more details.
14
15 // You should have received a copy of the GNU General Public License
16 // along with this package; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 /*################################################################################################*/
20
21 #include "libghemicalconfig2.h"
22 #include "utility.h"
23
24 #include "v3d.h"
25 #include "notice.h"
26
27 #include "local_i18n.h"
28
29 #include <cstring>
30 #include <iomanip>
31 #include <iostream>
32 #include <sstream>
33 using namespace std;
34
35 /*################################################################################################*/
36
superimpose(model * p1,i32s p2,i32s p3)37 superimpose::superimpose(model * p1, i32s p2, i32s p3) : conjugate_gradient(10, 1.0e-5, 1.0e+3) // 2003-06-08 ok!!!
38 {
39 mdl = p1;
40 index[0] = p2;
41 index[1] = p3;
42
43 for (i32s n1 = 0;n1 < 3;n1++)
44 {
45 rot[n1] = drot[n1] = 0.0;
46 AddVar(& rot[n1], & drot[n1]);
47
48 loc[n1] = dloc[n1] = 0.0;
49 AddVar(& loc[n1], & dloc[n1]);
50 }
51 }
52
~superimpose(void)53 superimpose::~superimpose(void)
54 {
55 }
56
GetRMS(void)57 f64 superimpose::GetRMS(void)
58 {
59 return sqrt(value / (f64) counter);
60 }
61
Compare(const f64 * d1,const f64 * d2,bool gradient,f64 * result)62 void superimpose::Compare(const f64 * d1, const f64 * d2, bool gradient, f64 * result)
63 {
64 const f64 balance = 0.25;
65
66 // translate...
67 // translate...
68 // translate...
69
70 f64 d2a[3];
71 d2a[0] = d2[0] + loc[0] * balance;
72 d2a[1] = d2[1] + loc[1] * balance;
73 d2a[2] = d2[2] + loc[2] * balance;
74
75 // rotate x (y,z)...
76 // rotate x (y,z)...
77 // rotate x (y,z)...
78
79 f64 d2b[3];
80 d2b[0] = d2a[0];
81 d2b[1] = d2a[1] * cos(rot[0]) - d2a[2] * sin(rot[0]);
82 d2b[2] = d2a[1] * sin(rot[0]) + d2a[2] * cos(rot[0]);
83
84 // rotate y (z,x)...
85 // rotate y (z,x)...
86 // rotate y (z,x)...
87
88 f64 d2c[3];
89 d2c[0] = d2b[2] * sin(rot[1]) + d2b[0] * cos(rot[1]);
90 d2c[1] = d2b[1];
91 d2c[2] = d2b[2] * cos(rot[1]) - d2b[0] * sin(rot[1]);
92
93 // rotate z (x,y)...
94 // rotate z (x,y)...
95 // rotate z (x,y)...
96
97 f64 d2d[3];
98 d2d[0] = d2c[0] * cos(rot[2]) - d2c[1] * sin(rot[2]);
99 d2d[1] = d2c[0] * sin(rot[2]) + d2c[1] * cos(rot[2]);
100 d2d[2] = d2c[2];
101
102 // rot-gradients...
103 // rot-gradients...
104 // rot-gradients...
105
106 f64 drz[3];
107 drz[0] = -(d2c[0] * sin(rot[2]) + d2c[1] * cos(rot[2]));
108 drz[1] = d2c[0] * cos(rot[2]) - d2c[1] * sin(rot[2]);
109 drz[2] = 0.0;
110
111 f64 dry[3];
112 dry[0] = (d2b[2] * cos(rot[1]) - d2b[0] * sin(rot[1])) * cos(rot[2]);
113 dry[1] = (d2b[2] * cos(rot[1]) - d2b[0] * sin(rot[1])) * sin(rot[2]);
114 dry[2] = -(d2b[2] * sin(rot[1]) + d2b[0] * cos(rot[1]));
115
116 f64 drx[3];
117 drx[0] = ((d2a[1] * cos(rot[0]) - d2a[2] * sin(rot[0])) * sin(rot[1])) * cos(rot[2]) +
118 (d2a[1] * sin(rot[0]) + d2a[2] * cos(rot[0])) * sin(rot[2]);
119 drx[1] = ((d2a[1] * cos(rot[0]) - d2a[2] * sin(rot[0])) * sin(rot[1])) * sin(rot[2]) -
120 (d2a[1] * sin(rot[0]) + d2a[2] * cos(rot[0])) * cos(rot[2]);
121 drx[2] = (d2a[1] * cos(rot[0]) - d2a[2] * sin(rot[0])) * cos(rot[1]);
122
123 // loc-gradients...
124 // loc-gradients...
125 // loc-gradients...
126
127 f64 dlx[3];
128 dlx[0] = balance * cos(rot[1]) * cos(rot[2]);
129 dlx[1] = balance * cos(rot[1]) * sin(rot[2]);
130 dlx[2] = -balance * sin(rot[1]);
131
132 f64 dly[3];
133 dly[0] = balance * sin(rot[0]) * sin(rot[1]) * cos(rot[2]) - balance * cos(rot[0]) * sin(rot[2]);
134 dly[1] = balance * sin(rot[0]) * sin(rot[1]) * sin(rot[2]) + balance * cos(rot[0]) * cos(rot[2]);
135 dly[2] = balance * sin(rot[0]) * cos(rot[1]);
136
137 f64 dlz[3];
138 dlz[0] = balance * cos(rot[0]) * sin(rot[1]) * cos(rot[2]) + balance * sin(rot[0]) * sin(rot[2]);
139 dlz[1] = balance * cos(rot[0]) * sin(rot[1]) * sin(rot[2]) - balance * sin(rot[0]) * cos(rot[2]);
140 dlz[2] = balance * cos(rot[0]) * cos(rot[1]);
141
142 // f = (a-b)^2
143 // df/db = -2(a-b)*Db
144
145 f64 diff[3];
146 diff[0] = d1[0] - d2d[0];
147 diff[1] = d1[1] - d2d[1];
148 diff[2] = d1[2] - d2d[2];
149
150 f64 tmp1 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2];
151
152 value += tmp1;
153
154 if (gradient)
155 {
156 dloc[0] -= 2.0 * diff[0] * dlx[0];
157 dloc[0] -= 2.0 * diff[1] * dlx[1];
158 dloc[0] -= 2.0 * diff[2] * dlx[2];
159
160 dloc[1] -= 2.0 * diff[0] * dly[0];
161 dloc[1] -= 2.0 * diff[1] * dly[1];
162 dloc[1] -= 2.0 * diff[2] * dly[2];
163
164 dloc[2] -= 2.0 * diff[0] * dlz[0];
165 dloc[2] -= 2.0 * diff[1] * dlz[1];
166 dloc[2] -= 2.0 * diff[2] * dlz[2];
167
168 drot[0] -= 2.0 * diff[0] * drx[0];
169 drot[0] -= 2.0 * diff[1] * drx[1];
170 drot[0] -= 2.0 * diff[2] * drx[2];
171
172 drot[1] -= 2.0 * diff[0] * dry[0];
173 drot[1] -= 2.0 * diff[1] * dry[1];
174 drot[1] -= 2.0 * diff[2] * dry[2];
175
176 drot[2] -= 2.0 * diff[0] * drz[0];
177 drot[2] -= 2.0 * diff[1] * drz[1];
178 drot[2] -= 2.0 * diff[2] * drz[2];
179 }
180
181 if (result != NULL)
182 {
183 result[0] = d2d[0];
184 result[1] = d2d[1];
185 result[2] = d2d[2];
186 }
187
188 counter++;
189 }
190
GetValue(void)191 f64 superimpose::GetValue(void)
192 {
193 value = 0.0; counter = 0;
194
195 for (iter_al it1 = mdl->GetAtomsBegin();it1 != mdl->GetAtomsEnd();it1++)
196 {
197 // skip the atoms with ATOMFLAG_IS_HIDDEN set, because those are not displayed either and might have invalid coordinates!!!
198 if ((* it1).flags & ATOMFLAG_IS_HIDDEN) continue;
199
200 // also skip the solvent atoms ; 20050120
201 if ((* it1).flags & ATOMFLAG_IS_SOLVENT_ATOM) continue;
202
203 const fGL * crd1 = (* it1).GetCRD(index[0]);
204 f64 d1[3] = { crd1[0], crd1[1], crd1[2] };
205
206 const fGL * crd2 = (* it1).GetCRD(index[1]);
207 f64 d2[3] = { crd2[0], crd2[1], crd2[2] };
208
209 Compare(d1, d2, false);
210 }
211
212 return value;
213 }
214
GetGradient(void)215 f64 superimpose::GetGradient(void)
216 {
217 value = 0.0; counter = 0;
218 dloc[0] = dloc[1] = dloc[2] = 0.0;
219 drot[0] = drot[1] = drot[2] = 0.0;
220
221 for (iter_al it1 = mdl->GetAtomsBegin();it1 != mdl->GetAtomsEnd();it1++)
222 {
223 // skip the atoms with ATOMFLAG_IS_HIDDEN set, because those are not displayed either and might have invalid coordinates!!!
224 if ((* it1).flags & ATOMFLAG_IS_HIDDEN) continue;
225
226 // also skip the solvent atoms ; 20050120
227 if ((* it1).flags & ATOMFLAG_IS_SOLVENT_ATOM) continue;
228
229 const fGL * crd1 = (* it1).GetCRD(index[0]);
230 f64 d1[3] = { crd1[0], crd1[1], crd1[2] };
231
232 const fGL * crd2 = (* it1).GetCRD(index[1]);
233 f64 d2[3] = { crd2[0], crd2[1], crd2[2] };
234
235 Compare(d1, d2, true);
236 }
237
238 return value;
239 }
240
Transform(void)241 void superimpose::Transform(void)
242 {
243 value = 0.0; counter = 0;
244
245 for (iter_al it1 = mdl->GetAtomsBegin();it1 != mdl->GetAtomsEnd();it1++)
246 {
247 // transform all atoms, no matter what flags are set...
248 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
249
250 const fGL * crd1 = (* it1).GetCRD(index[0]);
251 f64 d1[3] = { crd1[0], crd1[1], crd1[2] };
252
253 const fGL * crd2 = (* it1).GetCRD(index[1]);
254 f64 d2[3] = { crd2[0], crd2[1], crd2[2] };
255
256 f64 d3[3];
257 Compare(d1, d2, false, d3);
258
259 (* it1).SetCRD(index[1], d3[0], d3[1], d3[2]);
260 }
261 }
262
263 /*################################################################################################*/
264
265 #define HBOND_TRESHOLD -1.0 // h-bond energy treshold kcal/mol (K&S used -0.5 kcal/mol?)
266
267 #define HELIX_CHECK (M_PI * 45.0 / 180.0) // geometry tresholds for helix torsions...
268 #define HELIX4_REF (M_PI * +50.3 / 180.0)
269
DefineSecondaryStructure(model * mdl)270 void DefineSecondaryStructure(model * mdl)
271 {
272 vector<chn_info> & ci_vector = (* mdl->ref_civ);
273
274 for (i32u n1 = 0;n1 < ci_vector.size();n1++)
275 {
276 if (ci_vector[n1].type != chn_info::amino_acid) continue;
277
278 if (ci_vector[n1].ss_state != NULL) delete[] ci_vector[n1].ss_state;
279 ci_vector[n1].ss_state = new char[ci_vector[n1].length + 1];
280
281 f64 * energy = new f64[ci_vector[n1].length];
282
283 ci_vector[n1].ss_state[ci_vector[n1].length] = 0; // make it a null-terminated string for convenience...
284 for (i32s n2 = 0;n2 < ci_vector[n1].length;n2++)
285 {
286 ci_vector[n1].ss_state[n2] = '.';
287 energy[n2] = HBOND_TRESHOLD;
288 }
289
290 // the 3-turns and 5-turns are disabled (so far) ; 20050527
291 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
292
293 /* for (i32s n2 = 0;n2 < ci_vector[n1].length - 3;n2++) // 3-turns...
294 {
295 i32s tmp1[2] = { n1, n2 };
296 i32s tmp2[2] = { n1, n2 + 3 };
297
298 f64 tmp3 = HBondEnergy(mdl, tmp1, tmp2);
299 if (tmp3 < energy[n2])
300 {
301 ci_vector[n1].ss_state[n2] = '3';
302 energy[n2] = tmp3;
303 }
304 } */
305
306 for (i32s n2 = 0;n2 < ci_vector[n1].length - 4;n2++) // 4-turns...
307 {
308 i32s tmp1[2] = { n1, n2 };
309 i32s tmp2[2] = { n1, n2 + 4 };
310
311 f64 tmp3 = HBondEnergy(mdl, tmp1, tmp2);
312 if (tmp3 < energy[n2])
313 {
314 bool tor_check = true;
315 i32s tc1[4] = { n2+0, n2+1, n2+2, n2+3 }; if (!TorsionCheck(mdl, n1, tc1, HELIX4_REF)) tor_check = false;
316 i32s tc2[4] = { n2+1, n2+2, n2+3, n2+4 }; if (!TorsionCheck(mdl, n1, tc2, HELIX4_REF)) tor_check = false;
317
318 if (tor_check)
319 {
320 ci_vector[n1].ss_state[n2] = '4';
321 energy[n2] = tmp3;
322 }
323 }
324 }
325
326 /* for (i32s n2 = 0;n2 < ci_vector[n1].length - 5;n2++) // 5-turns...
327 {
328 i32s tmp1[2] = { n1, n2 };
329 i32s tmp2[2] = { n1, n2 + 5 };
330
331 f64 tmp3 = HBondEnergy(mdl, tmp1, tmp2);
332 if (tmp3 < energy[n2])
333 {
334 ci_vector[n1].ss_state[n2] = '5';
335 energy[n2] = tmp3;
336 }
337 } */
338
339 delete[] energy;
340
341 cout << _("looking for intrachain strands for chain ") << (n1 + 1) << endl;
342
343 for (i32s n2 = 1;n2 < ci_vector[n1].length - 4;n2++)
344 {
345 cout << "." << flush;
346
347 for (i32s n3 = n2 + 3;n3 < ci_vector[n1].length - 1;n3++)
348 {
349 i32s t1a[2] = { n1, n2 - 1 }; i32s t1b[2] = { n1, n3 };
350 f64 t1c = HBondEnergy(mdl, t1a, t1b);
351
352 i32s t2a[2] = { n1, n3 }; i32s t2b[2] = { n1, n2 + 1 };
353 f64 t2c = HBondEnergy(mdl, t2a, t2b);
354
355 i32s t3a[2] = { n1, n3 - 1 }; i32s t3b[2] = { n1, n2 };
356 f64 t3c = HBondEnergy(mdl, t3a, t3b);
357
358 i32s t4a[2] = { n1, n2 }; i32s t4b[2] = { n1, n3 + 1 };
359 f64 t4c = HBondEnergy(mdl, t4a, t4b);
360
361 i32s t5a[2] = { n1, n2 }; i32s t5b[2] = { n1, n3 };
362 f64 t5c = HBondEnergy(mdl, t5a, t5b);
363
364 i32s t6a[2] = { n1, n3 }; i32s t6b[2] = { n1, n2 };
365 f64 t6c = HBondEnergy(mdl, t6a, t6b);
366
367 i32s t7a[2] = { n1, n2 - 1 }; i32s t7b[2] = { n1, n3 + 1 };
368 f64 t7c = HBondEnergy(mdl, t7a, t7b);
369
370 i32s t8a[2] = { n1, n3 - 1 }; i32s t8b[2] = { n1, n2 + 1 };
371 f64 t8c = HBondEnergy(mdl, t8a, t8b);
372
373 bool bridge = false;
374 if (t1c < HBOND_TRESHOLD && t2c < HBOND_TRESHOLD) bridge = true; // parallel
375 if (t3c < HBOND_TRESHOLD && t4c < HBOND_TRESHOLD) bridge = true; // parallel
376 if (t5c < HBOND_TRESHOLD && t6c < HBOND_TRESHOLD) bridge = true; // antiparallel
377 if (t7c < HBOND_TRESHOLD && t8c < HBOND_TRESHOLD) bridge = true; // antiparallel
378
379 if (bridge)
380 {
381 ci_vector[n1].ss_state[n2] = 'S';
382 ci_vector[n1].ss_state[n3] = 'S';
383 }
384 }
385 }
386
387 cout << endl;
388 }
389
390 cout << _("looking for interchain strands");
391
392 for (i32s n1 = 0;n1 < ((i32s) ci_vector.size()) - 1;n1++)
393 {
394 if (ci_vector[n1].type != chn_info::amino_acid) continue;
395
396 for (i32s n2 = n1 + 1;n2 < (i32s) ci_vector.size();n2++)
397 {
398 if (ci_vector[n2].type != chn_info::amino_acid) continue;
399
400 cout << "." << flush;
401
402 for (i32s n3 = 1;n3 < ci_vector[n1].length - 1;n3++)
403 {
404 for (i32s n4 = 1;n4 < ci_vector[n2].length - 1;n4++)
405 {
406 i32s t1a[2] = { n1, n3 - 1 }; i32s t1b[2] = { n2, n4 };
407 f64 t1c = HBondEnergy(mdl, t1a, t1b);
408
409 i32s t2a[2] = { n2, n4 }; i32s t2b[2] = { n1, n3 + 1 };
410 f64 t2c = HBondEnergy(mdl, t2a, t2b);
411
412 i32s t3a[2] = { n2, n4 - 1 }; i32s t3b[2] = { n1, n3 };
413 f64 t3c = HBondEnergy(mdl, t3a, t3b);
414
415 i32s t4a[2] = { n1, n3 }; i32s t4b[2] = { n2, n4 + 1 };
416 f64 t4c = HBondEnergy(mdl, t4a, t4b);
417
418 i32s t5a[2] = { n1, n3 }; i32s t5b[2] = { n2, n4 };
419 f64 t5c = HBondEnergy(mdl, t5a, t5b);
420
421 i32s t6a[2] = { n2, n4 }; i32s t6b[2] = { n1, n3 };
422 f64 t6c = HBondEnergy(mdl, t6a, t6b);
423
424 i32s t7a[2] = { n1, n3 - 1 }; i32s t7b[2] = { n2, n4 + 1 };
425 f64 t7c = HBondEnergy(mdl, t7a, t7b);
426
427 i32s t8a[2] = { n2, n4 - 1 }; i32s t8b[2] = { n1, n3 + 1 };
428 f64 t8c = HBondEnergy(mdl, t8a, t8b);
429
430 bool bridge = false;
431 if (t1c < HBOND_TRESHOLD && t2c < HBOND_TRESHOLD) bridge = true; // see above???
432 if (t3c < HBOND_TRESHOLD && t4c < HBOND_TRESHOLD) bridge = true; // see above???
433 if (t5c < HBOND_TRESHOLD && t6c < HBOND_TRESHOLD) bridge = true; // see above???
434 if (t7c < HBOND_TRESHOLD && t8c < HBOND_TRESHOLD) bridge = true; // see above???
435
436 if (bridge)
437 {
438 ci_vector[n1].ss_state[n3] = 'S';
439 ci_vector[n2].ss_state[n4] = 'S';
440 }
441 }
442 }
443 }
444 }
445
446 cout << endl;
447
448 for (i32u n1 = 0;n1 < ci_vector.size();n1++)
449 {
450 if (ci_vector[n1].type != chn_info::amino_acid) continue;
451
452 cout << _("found chain ") << n1 << " :" << endl;
453
454 i32u length = strlen(ci_vector[n1].ss_state);
455 for (i32u n2 = 0;n2 < length;n2++)
456 {
457 cout << ci_vector[n1].ss_state[n2];
458
459 bool is_break = !((n2 + 1) % 50);
460 bool is_end = ((n2 + 1) == length);
461
462 if (is_break || is_end) cout << endl;
463 }
464
465 cout << endl;
466 }
467
468 cout << _("DefineSecondaryStructure() is ready.") << endl;
469 }
470
471 // here can be some problems with the carboxyl terminus where
472 // we can accidentally pick "wrong" oxygen...
473
474 // for simplified residues, no changes are needed here?!?!?!
475 // in sf cases 0.0 is just returned since no O/N atoms are found.
476
477 // res1 = H-bond acceptor ; res2 = H-bond donor.
478
HBondEnergy(model * mdl,i32s * res1,i32s * res2)479 f64 HBondEnergy(model * mdl, i32s * res1, i32s * res2)
480 {
481 vector<chn_info> & ci_vector = (* mdl->ref_civ);
482
483 const char * seq = ci_vector[res2[0]].GetSequence1();
484 if (seq[res2[1]] == 'P') return 0.0; // handle the PRO case...
485
486 iter_al r1a[2]; mdl->GetRange(1, res1[0], r1a);
487 iter_al r1b[2]; mdl->GetRange(2, r1a, res1[1], r1b);
488
489 iter_al r2a[2]; mdl->GetRange(1, res2[0], r2a);
490 iter_al r2b[2]; mdl->GetRange(2, r2a, res2[1], r2b);
491
492 iter_al it_o = r1b[0];
493 while (it_o != r1b[1] && ((* it_o).builder_res_id & 0xFF) != 0x10) it_o++;
494 if (it_o == r1b[1]) return 0.0;
495
496 iter_al it_n = r2b[0];
497 while (it_n != r2b[1] && ((* it_n).builder_res_id & 0xFF) != 0x00) it_n++;
498 if (it_n == r2b[1]) return 0.0;
499
500 // if the distance between the residues is large enough we can be sure that there
501 // is no h-bond and skip this time-consumig process. according to K&S this treshold
502 // NO-distance is about 6.0 �, but let's be careful and use a bit larger value...
503
504 v3d<fGL> sv = v3d<fGL>((* it_o).GetCRD(0), (* it_n).GetCRD(0));
505 if (sv.len() > 0.75) return 0.0;
506
507 // the above check is not enough to make this fast. problem is that the check is done in
508 // too late stage to be efficient. calculate first distances between all residues and
509 // use those results for sanity checking in earlier stage???
510
511 iter_al it_c = r1b[0];
512 while (it_c != r1b[1] && ((* it_c).builder_res_id & 0xFF) != 0x02) it_c++;
513 if (it_c == r1b[1]) return 0.0;
514
515 atom * ref_h = NULL;
516 list<crec>::iterator it1 = (* it_n).GetConnRecsBegin();
517 /* while (it1 != (* it_n).GetConnRecsEnd() && ref_h == NULL)
518 {
519 if ((* it1).atmr->el.GetAtomicNumber() != 1) it1++;
520 else ref_h = (* it1).atmr;
521 } */
522
523 // the sequence builder currently ignores hydrogen atoms, so those
524 // are not always available. also PRO does not have the proton (see ABOVE).
525
526 // therefore we don't use the real hydrogens here at all, but we "reconstruct"
527 // the amide hydrogen using the main-chain heavy atoms...
528
529 it1 = (* it_n).GetConnRecsEnd(); // pretend that we couldn't find the amide hydrogen...
530 if (it1 == (* it_n).GetConnRecsEnd())
531 {
532 if (!res2[1]) return 0.0; // skip if N-terminal residue...
533 iter_al r2c[2]; mdl->GetRange(2, r2a, res2[1] - 1, r2c);
534
535 iter_al it_c1 = r2b[0];
536 while (it_c1 != r2b[1] && ((* it_c1).builder_res_id & 0xFF) != 0x01) it_c1++;
537 if (it_c1 == r2b[1]) return 0.0;
538
539 iter_al it_c2 = r2c[0];
540 while (it_c2 != r2c[1] && ((* it_c2).builder_res_id & 0xFF) != 0x02) it_c2++;
541 if (it_c2 == r2c[1]) return 0.0;
542
543 v3d<fGL> v1 = v3d<fGL>((* it_n).GetCRD(0), (* it_c1).GetCRD(0));
544 v3d<fGL> v2 = v3d<fGL>((* it_n).GetCRD(0), (* it_c2).GetCRD(0));
545 v3d<fGL> v3 = v1.vpr(v2); v3d<fGL> v4 = v1.vpr(v3);
546
547 const fGL len = 0.1024;
548 const fGL ang = M_PI * 121.0 / 180.0;
549 v1 = v1 * (len * cos(ang) / v1.len());
550 v4 = v4 * (len * sin(ang) / v4.len());
551 v3d<fGL> v5 = v3d<fGL>((* it_n).GetCRD(0)) + (v1 + v4);
552 ref_h = new atom(element(1), v5.data, 1);
553
554 ///////////////////////////////////////////////////////////////////////////
555 //v3d<fGL> tv1 = v3d<fGL>((* it_n).crd[0], v5.comp);
556 //cout << "the H bond length = " << tv1.len() << endl;
557 //cout << "the H bond angle = " << (tv1.ang(v2) * 180.0 / M_PI) << endl;
558 //i32s zzz; cin >> zzz;
559 ///////////////////////////////////////////////////////////////////////////
560 }
561
562 v3d<fGL> on_v = v3d<fGL>((* it_o).GetCRD(0), (* it_n).GetCRD(0));
563 v3d<fGL> ch_v = v3d<fGL>((* it_c).GetCRD(0), ref_h->GetCRD(0));
564
565 v3d<fGL> oh_v = v3d<fGL>((* it_o).GetCRD(0), ref_h->GetCRD(0));
566 v3d<fGL> cn_v = v3d<fGL>((* it_c).GetCRD(0), (* it_n).GetCRD(0));
567
568 if (it1 == (* it_n).GetConnRecsEnd()) delete ref_h;
569
570 f64 tmp1 = 1.0 / on_v.len() + 1.0 / ch_v.len();
571 f64 tmp2 = 1.0 / oh_v.len() + 1.0 / cn_v.len();
572 f64 tmp3 = 0.42 * 0.20 * (tmp1 - tmp2);
573 return 0.1 * 332.0 * tmp3;
574 }
575
576 // for helices, simply checking the above test seem to give false positives.
577 // perform additional main-chain torsion checks here...
578
TorsionCheck(model * mdl,i32s chn,i32s * res,fGL ref)579 bool TorsionCheck(model * mdl, i32s chn, i32s * res, fGL ref)
580 {
581 iter_al rc[2]; mdl->GetRange(1, chn, rc);
582 iter_al rr1[2]; mdl->GetRange(2, rc, res[0], rr1);
583 iter_al rr2[2]; mdl->GetRange(2, rc, res[1], rr2);
584 iter_al rr3[2]; mdl->GetRange(2, rc, res[2], rr3);
585 iter_al rr4[2]; mdl->GetRange(2, rc, res[3], rr4);
586
587 iter_al it_c1 = rr1[0];
588 while (it_c1 != rr1[1] && ((* it_c1).builder_res_id & 0xFF) != 0x02) it_c1++;
589 if (it_c1 == rr1[1]) assertion_failed(__FILE__, __LINE__, "c_alpha #1 not found.");
590
591 iter_al it_c2 = rr2[0];
592 while (it_c2 != rr2[1] && ((* it_c2).builder_res_id & 0xFF) != 0x02) it_c2++;
593 if (it_c2 == rr2[1]) assertion_failed(__FILE__, __LINE__, "c_alpha #2 not found.");
594
595 iter_al it_c3 = rr3[0];
596 while (it_c3 != rr3[1] && ((* it_c3).builder_res_id & 0xFF) != 0x02) it_c3++;
597 if (it_c3 == rr3[1]) assertion_failed(__FILE__, __LINE__, "c_alpha #3 not found.");
598
599 iter_al it_c4 = rr4[0];
600 while (it_c4 != rr4[1] && ((* it_c4).builder_res_id & 0xFF) != 0x02) it_c4++;
601 if (it_c4 == rr4[1]) assertion_failed(__FILE__, __LINE__, "c_alpha #4 not found.");
602
603 v3d<fGL> v1 = v3d<fGL>((* it_c2).GetCRD(0), (* it_c1).GetCRD(0));
604 v3d<fGL> v2 = v3d<fGL>((* it_c2).GetCRD(0), (* it_c3).GetCRD(0));
605 v3d<fGL> v3 = v3d<fGL>((* it_c3).GetCRD(0), (* it_c4).GetCRD(0));
606 fGL tor = v1.tor(v2, v3);
607
608 fGL dev = tor - ref;
609 if (dev > +M_PI) dev = 2.0 * M_PI - dev;
610 else if (dev < -M_PI) dev = 2.0 * M_PI + dev;
611
612 if (fabs(dev) < HELIX_CHECK) return true;
613 else
614 {
615 cout << _("HELIX CHECK FAILED : ") << dev << endl;
616 // int pause; cin >> pause;
617
618 return false;
619 }
620 }
621
622 /*################################################################################################*/
623
624 // eof
625