1 // MFINDER.CPP
2
3 // Copyright (C) 2006 Tommi Hassinen.
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 "mfinder.h"
23
24 #include "v3d.h"
25
26 #include "local_i18n.h"
27 #include "notice.h"
28
29 #include <cstring>
30 #include <iomanip>
31 #include <iostream>
32 #include <sstream>
33 using namespace std;
34
35 /*################################################################################################*/
36
mfinder(char * filename)37 mfinder::mfinder(char * filename)
38 {
39 ifstream file;
40 file.unsetf(ios::dec | ios::oct | ios::hex);
41
42 model::OpenLibDataFile(file, false, filename);
43
44 // read information about the main chain...
45
46 while (file.peek() != 'M') file.getline(buffer, sizeof(buffer));
47 file.getline(buffer, sizeof(buffer));
48
49 while (file.peek() != 'E')
50 {
51 mf_data_atm newatm; file >> newatm;
52 while (file.peek() != '(') file.get();
53 newatm.tr = new typerule(& file, & cout);
54 file.getline(buffer, sizeof(buffer));
55 main_vector.push_back(newatm);
56 }
57
58 while (file.peek() != 'C') file.getline(buffer, sizeof(buffer));
59 file.getline(buffer, sizeof(buffer));
60
61 while (file.peek() != 'E')
62 {
63 mf_data_atm newatm; file >> newatm;
64 while (file.peek() != '(') file.get();
65 newatm.tr = new typerule(& file, & cout);
66 file.getline(buffer, sizeof(buffer));
67 conn_vector.push_back(newatm);
68 }
69
70 // read descriptions about the chain initiation...
71
72 while (file.peek() != 'H') file.getline(buffer, sizeof(buffer));
73 file.getline(buffer, sizeof(buffer));
74
75 while (file.peek() != 'E')
76 {
77 typerule newrule = typerule(& file, & cout);
78 file.getline(buffer, sizeof(buffer));
79 head_vector.push_back(newrule);
80 }
81
82 // read descriptions about the chain termination...
83
84 while (file.peek() != 'T') file.getline(buffer, sizeof(buffer));
85 file.getline(buffer, sizeof(buffer));
86
87 while (file.peek() != 'E')
88 {
89 typerule newrule = typerule(& file, & cout);
90 file.getline(buffer, sizeof(buffer));
91 tail_vector.push_back(newrule);
92 }
93
94 // read the default modifications for residues...
95
96 // while (file.peek() != 'B') file.getline(buffer, sizeof(buffer)); // default...
97 // file.getline(buffer, sizeof(buffer)); mod[0] = new mb_data_res; mod[0]->ReadModification(file);
98
99 // while (file.peek() != 'H') file.getline(buffer, sizeof(buffer)); // head-mods
100 // file.getline(buffer, sizeof(buffer)); mod[1] = new mb_data_res; mod[1]->ReadModification(file);
101
102 // while (file.peek() != 'T') file.getline(buffer, sizeof(buffer)); // tail-mods
103 // file.getline(buffer, sizeof(buffer)); mod[2] = new mb_data_res; mod[2]->ReadModification(file);
104
105 // read descriptions for the residues...
106
107 while (file.peek() != 'E')
108 {
109 if (file.peek() == 'R')
110 {
111 // mb_data_res newres; file >> newres;
112 // resi_vector.push_back(newres);
113 }
114 else file.getline(buffer, sizeof(buffer));
115 }
116
117 // ready!!!
118
119 file.close();
120 }
121
~mfinder(void)122 mfinder::~mfinder(void)
123 {
124 // delete mod[0];
125 // delete mod[1];
126 // delete mod[2];
127 }
128
129 // here we will identify the molecules (if there is any), write down the molecules
130 // into model::ref_civ and update the /mol/chn/res-numbering...
131
132 // only model::UpdateChains() should call this???
133 // only model::UpdateChains() should call this???
134 // only model::UpdateChains() should call this???
135
Identify(model * mdl)136 void mfinder::Identify(model * mdl)
137 {
138 if (!mdl->IsGroupsClean()) assertion_failed(__FILE__, __LINE__, "!mdl->IsGroupsClean()");
139
140 if (mdl->ref_civ == NULL) assertion_failed(__FILE__, __LINE__, "mdl->ref_civ == NULL");
141
142 cout << "mfinder::Identify() starts..." << endl;
143 cout << "nmol = " << mdl->nmol << endl;
144
145 // here we will find all possible chains and identify them.
146
147 for (i32s n1 = 0;n1 < mdl->nmol;n1++)
148 {
149 iter_al range[2];
150 mdl->GetRange(0, n1, range);
151
152 vector<atom *> head_atoms;
153 vector<atom *> tail_atoms;
154
155 for (iter_al it1 = range[0];it1 != range[1];it1++)
156 {
157 i32s tmp1 = (* it1).el.GetAtomicNumber();
158
159 if (tmp1 == main_vector.front().el.GetAtomicNumber()) // look for heads...
160 {
161 for (i32u n2 = 0;n2 < head_vector.size();n2++)
162 {
163 bool flag = head_vector[n2].Check(mdl, & (* it1), 0);
164 if (!flag) continue; head_atoms.push_back(& (* it1)); break;
165 }
166 }
167
168 if (tmp1 == main_vector.back().el.GetAtomicNumber()) // look for tails...
169 {
170 for (i32u n2 = 0;n2 < tail_vector.size();n2++)
171 {
172 bool flag = tail_vector[n2].Check(mdl, & (* it1), 0);
173 if (!flag) continue; tail_atoms.push_back(& (* it1)); break;
174 }
175 }
176 }
177
178 if (head_atoms.size() > 0 && tail_atoms.size() > 0)
179 {
180 cout << _("found ") << head_atoms.size() << _(" possible heads and ");
181 cout << tail_atoms.size() << _(" possible tails.") << endl;
182 }
183
184 // now we have all possible head/tail atoms. next we have to check
185 // all possible paths between them to find all possible main chains...
186
187 path_vector.resize(0);
188 for (i32u n2 = 0;n2 < head_atoms.size();n2++)
189 {
190 for (i32u n3 = 0;n3 < tail_atoms.size();n3++)
191 {
192 FindPath(mdl, head_atoms[n2], tail_atoms[n3]);
193 }
194 }
195
196 if (path_vector.size())
197 {
198 cout << path_vector.size() << _(" chains:") << endl;
199 }
200
201 for (i32s n2 = 0;n2 < (i32s) path_vector.size();n2++)
202 {
203 for (i32s n3 = 0;n3 < ((i32s) path_vector[n2].size()) - 1;n3++) // tag the main-chain bonds...
204 {
205 iter_cl it1 = path_vector[n2][n3]->GetConnRecsBegin();
206 while ((* it1).atmr != path_vector[n2][n3 + 1]) it1++;
207 (* it1).bndr->flags[0] = true;
208 }
209
210 //seq vector<char> molecule;
211
212 i32u acount = 0; // atom counter
213 i32u rcount = 0; // residue counter
214
215 while (acount < path_vector[n2].size())
216 {
217 // rsize is the residue size (how many atoms it will chop away from the
218 // path_vector). we use it to determine whether we have reached the last
219 // residue, and to increment acount.
220
221 int head_el = NOT_DEFINED;
222 int tail_el = NOT_DEFINED;
223
224 bool is_first = !rcount;
225 // if (is_first && type == chn_info::amino_acid)
226 // {
227 // // an extra test for an ACE-like block group.
228 // // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
229 //
230 // atom * refA = path_vector[n2][acount + 0];
231 // atom * refB = path_vector[n2][acount + 1];
232 // atom * refX = NULL;
233 //
234 // for (iter_cl itX = refA->cr_list.begin();itX != refA->cr_list.end();itX++)
235 // {
236 // if ((* itX).atmr->el.GetAtomicNumber() <= 1) continue;
237 // if ((* itX).atmr == refB) continue;
238 //
239 // refX = (* itX).atmr;
240 // }
241 //
242 // if (refX != NULL)
243 // {
244 // head_el = refX->el.GetAtomicNumber();
245 // is_first = false;
246 // }
247 // }
248
249 bool is_last = (acount + main_vector.size() + conn_vector.size()) >= path_vector[n2].size();
250 // if (is_last && type == chn_info::amino_acid)
251 // {
252 // // an extra test for an NME-like block group.
253 // // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
254 //
255 // atom * refA = path_vector[n2][acount + 2];
256 // atom * refB = path_vector[n2][acount + 1];
257 // atom * refX = NULL;
258 //
259 // for (iter_cl itX = refA->cr_list.begin();itX != refA->cr_list.end();itX++)
260 // {
261 // if ((* itX).bndr->bt.GetValue() != 1) continue;
262 // if ((* itX).atmr == refB) continue;
263 //
264 // refX = (* itX).atmr;
265 // }
266 //
267 // if (refX != NULL)
268 // {
269 // tail_el = refX->el.GetAtomicNumber();
270 // is_last = false;
271 // }
272 // }
273
274 i32u rsize = main_vector.size();
275 if (!is_last) rsize += conn_vector.size();
276
277 // now we start identifying the residues.
278
279 // we make a simple template for each of them, compare those to
280 // the residue descriptions we have, and pick the closest match.
281 // there may not be any missing atoms, and the closest match is
282 // the one with largest number of correct atoms identified.
283
284 i32s tmp1[2] = { NOT_DEFINED, NOT_DEFINED };
285 /* for (i32s n4 = 0;n4 < (i32s) resi_vector.size();n4++)
286 {
287 vector<mb_tdata> tdata;
288 BuildTemplate(tdata, n4, is_first, is_last);
289
290 // do the template mods if blockgroups!!!
291 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
292
293 if (head_el != NOT_DEFINED)
294 {
295 if (is_first)
296 {
297 cout << "bad head_el" << endl;
298 exit(EXIT_FAILURE);
299 }
300
301 // ACE needs no modifications...
302 }
303
304 if (tail_el != NOT_DEFINED)
305 {
306 if (is_last)
307 {
308 cout << "bad tail_el" << endl;
309 exit(EXIT_FAILURE);
310 }
311
312 for (i32s n5 = 0;n5 < (i32s) tdata.size();n5++)
313 {
314 if (tdata[n5].id[0] == 0x11)
315 {
316 tdata[n5].el = element(tail_el);
317 }
318 }
319 }
320
321 // set the main chain stuff...
322
323 for (i32s n5 = 0;n5 < (i32s) main_vector.size();n5++)
324 {
325 tdata[n5].ref = path_vector[n2][acount + n5];
326 }
327
328 // set the connection stuff...
329
330 if (!is_last)
331 {
332 const int mvsz = main_vector.size();
333 for (i32s n5 = 0;n5 < (i32s) conn_vector.size();n5++)
334 {
335 tdata[mvsz + n5].ref = path_vector[n2][acount + mvsz + n5];
336 }
337 }
338
339 bool result = CheckTemplate(tdata, 0);
340 ///////////////////////////////////////////////////////////////////////////////////////////
341 //cout << "checktemplate " << n4 << " " << is_first << " " << is_last << " ; ";
342 //cout << "result = " << result << endl;
343 //int zzz; cin >> zzz;
344 ///////////////////////////////////////////////////////////////////////////////////////////
345 if (result && (tmp1[0] == NOT_DEFINED || ((i32s) tdata.size()) > tmp1[1]))
346 {
347 tmp1[0] = n4;
348 tmp1[1] = (i32s) tdata.size();
349 }
350 } */
351
352 if (tmp1[0] == NOT_DEFINED)
353 {
354 cout << _("WARNING : residue ") << rcount << _(" was of unknown type!!!") << endl;
355
356 //seq molecule.push_back('?'); // handle the molecule...
357
358 for (i32u n3 = 0;n3 < main_vector.size();n3++) // handle the numbers...
359 {
360 path_vector[n2][acount + n3]->builder_res_id = NOT_DEFINED; // WHAT HERE???
361
362 path_vector[n2][acount + n3]->id[1] = mdl->ref_civ->size();
363 path_vector[n2][acount + n3]->id[2] = rcount;
364 }
365 }
366 else
367 {
368 // molecule.push_back(resi_vector[tmp1[0]].symbol1);
369
370 // rebuild the best matching template and use that to
371 // set up the atom ID numbers...
372
373 vector<mf_tdata> tdata;
374 BuildTemplate(tdata, tmp1[0], is_first, is_last);
375
376 for (i32s n4 = 0;n4 < (i32s) main_vector.size();n4++)
377 {
378 tdata[n4].ref = path_vector[n2][acount + n4];
379 }
380
381 if (!is_last)
382 {
383 const int mvsz = main_vector.size();
384 for (i32s n5 = 0;n5 < (i32s) conn_vector.size();n5++)
385 {
386 tdata[mvsz + n5].ref = path_vector[n2][acount + mvsz + n5];
387 }
388 }
389
390 CheckTemplate(tdata, 0);
391 for (i32s n4 = 0;n4 < (i32s) tdata.size();n4++)
392 {
393 // tdata[n4].ref->builder_res_id = (resi_vector[tmp1[0]].id << 8) + tdata[n4].id[0];
394
395 // tdata[n4].ref->id[1] = mdl->ref_civ->size();
396 // tdata[n4].ref->id[2] = rcount;
397 }
398
399 /////////////////////////////////////////////////////////////////////////// exceptions
400 /////////////////////////////////////////////////////////////////////////// exceptions
401
402 }
403
404 rcount++;
405 acount += rsize;
406 }
407
408 // // next we will write down the molecule...
409 //
410 // chn_info newinfo(type, molecule.size());
411 //
412 // newinfo.id_mol = n1;
413 // newinfo.id_chn = mdl->ref_civ->size();
414 //
415 // for (i32u n3 = 0;n3 < molecule.size();n3++) newinfo.molecule[n3] = molecule[n3];
416 //
417 // mdl->ref_civ->push_back(newinfo);
418 //
419 // // and make some output (using FASTA format again)...
420 //
421 // cout << "> chain " << newinfo.id_chn;
422 // cout << ", length " << newinfo.length << ":" << endl;
423 //
424 // for (i32u n3 = 0;n3 < molecule.size();n3++)
425 // {
426 // cout << molecule[n3];
427 //
428 // bool is_break = !((n3 + 1) % 50);
429 // bool is_end = ((n3 + 1) == molecule.size());
430 //
431 // if (is_break || is_end) cout << endl;
432 // }
433 //
434 // // finally un-tag the main-chain bonds...
435 //
436 // for (i32s n3 = 0;n3 < ((i32s) path_vector[n2].size()) - 1;n3++)
437 // {
438 // iter_cl it1 = path_vector[n2][n3]->cr_list.begin();
439 // while ((* it1).atmr != path_vector[n2][n3 + 1]) it1++;
440 // (* it1).bndr->flags[0] = false;
441 // }
442 }
443 // obsolete obsolete obsolete obsolete obsolete obsolete obsolete obsolete obsolete obsolete
444 // obsolete obsolete obsolete obsolete obsolete obsolete obsolete obsolete obsolete obsolete
445 /* now check the HYDROGENS that are ignored by the seqbuilder generally.
446 for hydrogens, set the same id's that the heavy atom has if they are valid ones.
447 for (iter_al it1 = range[0];it1 != range[1];it1++)
448 {
449 if ((* it1).el.GetAtomicNumber() != 1) continue;
450
451 if ((* it1).cr_list.size() != 1)
452 {
453 cout << "WARNING : seqbuild : H atom with abnormal connectivity found." << endl;
454 continue;
455 }
456
457 atom * heavy = (* it1).cr_list.front().atmr;
458 if (heavy->id[1] < 0) continue; // id < 0 means NOT_DEFINED...
459 if (heavy->id[2] < 0) continue; // id < 0 means NOT_DEFINED...
460
461 // ok, set the heavy atom id's to hydrogen as well...
462
463 (* it1).id[1] = heavy->id[1];
464 (* it1).id[2] = heavy->id[2];
465 } */
466 // obsolete obsolete obsolete obsolete obsolete obsolete obsolete obsolete obsolete obsolete
467 // obsolete obsolete obsolete obsolete obsolete obsolete obsolete obsolete obsolete obsolete
468 }
469
470 // traditionally, we sorted the atom list here to make the newly assigned chains in nice order.
471 // nowadays, this is done at model::UpdateChains() since there might be many different seqbuilder
472 // objects (like for amino/nucleic acid variants) and we like to use them all...
473 }
474
BuildTemplate(vector<mf_tdata> & tdata,i32s res,bool is_first,bool is_last)475 void mfinder::BuildTemplate(vector<mf_tdata> & tdata, i32s res, bool is_first, bool is_last)
476 {
477 BuildPartialT(tdata, main_vector);
478 if (!is_last) BuildPartialT(tdata, conn_vector);
479
480 // if (is_last) BuildPartialT(tdata, mod[2]->atm_vector); // tail-mods
481 // else if (is_first) BuildPartialT(tdata, mod[1]->atm_vector); // head-mods
482 // else BuildPartialT(tdata, mod[0]->atm_vector); // default...
483
484 // BuildPartialT(tdata, resi_vector[res].atm_vector);
485 }
486
BuildPartialT(vector<mf_tdata> & tdata,vector<mf_data_atm> & adata)487 void mfinder::BuildPartialT(vector<mf_tdata> & tdata, vector<mf_data_atm> & adata)
488 {
489 for (i32u n1 = 0;n1 < adata.size();n1++)
490 {
491 mf_tdata newdata;
492 newdata.id[0] = adata[n1].id;
493
494 if (adata[n1].prev[0] & 0xFF00) newdata.id[1] = NOT_DEFINED;
495 else newdata.id[1] = adata[n1].prev[0];
496
497 newdata.el = adata[n1].el;
498 newdata.bt = adata[n1].bt;
499 newdata.ref = NULL;
500
501 tdata.push_back(newdata);
502 }
503 }
504
505 // idea is that this path-searching function should be used only inside a molecule.
506 // this way we are sure to find always at least one path. we collect references to
507 // atoms for all paths and collect this way all possible paths.
508
509 // different flags should be used in type checking and in this path-searching...
510
FindPath(model * mdl,atom * ref1,atom * ref2,i32u index)511 void mfinder::FindPath(model * mdl, atom * ref1, atom * ref2, i32u index)
512 {
513 if (index >= main_vector.size() + conn_vector.size()) index = 0;
514
515 if (index < main_vector.size())
516 {
517 if (ref1->el.GetAtomicNumber() != main_vector[index].el.GetAtomicNumber()) return;
518 if (!main_vector[index].tr->Check(mdl, ref1, 0)) return;
519 }
520 else
521 {
522 i32u index2 = index - main_vector.size();
523
524 if (ref1->el.GetAtomicNumber() != conn_vector[index2].el.GetAtomicNumber()) return;
525 if (!conn_vector[index2].tr->Check(mdl, ref1, 0)) return;
526 }
527
528 temporary_vector.push_back(ref1);
529
530 //for (i32u debug = 0;debug < temporary_vector.size();debug++)
531 //{
532 // cout << temporary_vector[debug]->el.GetSymbol();
533 //} cout << endl;
534
535 if (ref1 == ref2) path_vector.push_back(temporary_vector);
536 else
537 {
538 for (iter_cl it1 = ref1->GetConnRecsBegin();it1 != ref1->GetConnRecsEnd();it1++)
539 {
540 if ((* it1).bndr->flags[2]) continue;
541
542 (* it1).bndr->flags[2] = true;
543 FindPath(mdl, (* it1).atmr, ref2, (index + 1));
544 (* it1).bndr->flags[2] = false;
545 }
546 }
547
548 temporary_vector.pop_back();
549 }
550
CheckTemplate(vector<mf_tdata> & tdata,i32s flag)551 bool mfinder::CheckTemplate(vector<mf_tdata> & tdata, i32s flag)
552 {
553 // miten antaa vinkki atomista nolla??? seh�n on "hotspot"???
554 // olisko templaatti "esit�ytetty" sen alun osalta???????????
555 vector<i32s> tmpv1;
556
557 // look for a suitable starting point to start matching in the template.
558 // if there is no such point we assume the template fits. however, the template
559 // could match since it is too small (for exaple, the template for glycine will
560 // match for all residues). therefore we have to test all templates and choose
561 // the biggest matching one!!!
562
563 i32u index = 0;
564 while (index < tdata.size())
565 {
566 if (tdata[index].ref != NULL)
567 {
568 tmpv1.resize(0);
569 for (i32u n1 = 0;n1 < tdata.size();n1++)
570 {
571 bool test1 = (tdata[n1].ref == NULL);
572 bool test2 = (tdata[n1].id[1] == tdata[index].id[0]);
573 if (test1 && test2) tmpv1.push_back(n1);
574 }
575
576 if (tmpv1.size()) break;
577 }
578
579 index++;
580 }
581
582 // if there is no suitable starting point, we have tested all atoms and
583 // none of them is missing -> template fits -> TRUE!!!
584
585 if (index == tdata.size()) return true;
586
587 // find all possible bonds that we can use to continue testing...
588
589 vector<crec> tmpv2; iter_cl it1;
590 for (it1 = tdata[index].ref->GetConnRecsBegin();it1 != tdata[index].ref->GetConnRecsEnd();it1++)
591 {
592 bool test1 = (* it1).bndr->flags[flag];
593 if (!test1) tmpv2.push_back((* it1));
594 }
595
596 // if there are less bonds than we need, the template can't match -> FALSE!!!
597
598 if (tmpv2.size() < tmpv1.size()) return false;
599
600 vector<i32s> tmpv3; tmpv3.resize(tmpv2.size());
601 for (i32u n1 = 0;n1 < tmpv3.size();n1++) tmpv3[n1] = n1;
602
603 // then we will check all bonds against the template in all possible permutations.
604 // if some combination is a match then leave the ID-numbers untouched, clean the
605 // bond flags and leave -> TRUE!!!
606
607 while (true)
608 {
609 bool test = true;
610 for (i32u n1 = 0;n1 < tmpv1.size();n1++)
611 {
612 i32s el = tdata[tmpv1[n1]].el.GetAtomicNumber();
613 if (el != NOT_DEFINED && el != tmpv2[tmpv3[n1]].atmr->el.GetAtomicNumber()) test = false;
614
615 i32s bt = tdata[tmpv1[n1]].bt.GetValue();
616 if (bt != NOT_DEFINED && bt != tmpv2[tmpv3[n1]].bndr->bt.GetValue()) test = false;
617 }
618
619 if (test)
620 {
621 for (i32u n1 = 0;n1 < tmpv1.size();n1++)
622 {
623 tmpv2[tmpv3[n1]].bndr->flags[flag] = true;
624 tdata[tmpv1[n1]].ref = tmpv2[tmpv3[n1]].atmr;
625 }
626
627 bool result = CheckTemplate(tdata, flag);
628
629 for (i32u n1 = 0;n1 < tmpv1.size();n1++)
630 {
631 tmpv2[tmpv3[n1]].bndr->flags[flag] = false;
632 if (!result) tdata[tmpv1[n1]].ref = NULL;
633 }
634
635 if (result) return true;
636 }
637
638 if (!next_permutation(tmpv3.begin(), tmpv3.end())) break;
639 }
640
641 // none of those combinations matched -> the template doesn't fit -> FALSE!!!
642
643 return false;
644 }
645
646 /*################################################################################################*/
647
mf_data_atm(void)648 mf_data_atm::mf_data_atm(void)
649 {
650 tr = NULL;
651 }
652
mf_data_atm(const mf_data_atm & p1)653 mf_data_atm::mf_data_atm(const mf_data_atm & p1)
654 {
655 id = p1.id;
656
657 el = p1.el;
658 bt = p1.bt;
659
660 ic2 = p1.ic2;
661
662 for (i32s n1 = 0;n1 < 3;n1++)
663 {
664 prev[n1] = p1.prev[n1];
665 ic1[n1] = p1.ic1[n1];
666 }
667
668 if (p1.tr != NULL) tr = new typerule(* p1.tr);
669 else tr = NULL;
670 }
671
~mf_data_atm(void)672 mf_data_atm::~mf_data_atm(void)
673 {
674 if (tr != NULL) delete tr;
675 }
676
operator >>(istream & p1,mf_data_atm & p2)677 istream & operator>>(istream & p1, mf_data_atm & p2)
678 {
679 char buffer[256];
680 while (p1.get() != 'M');
681 p1 >> p2.id; while (p1.get() != ':');
682 p1 >> buffer; p2.el = element(buffer);
683 p1 >> p2.prev[0] >> p2.prev[1] >> p2.prev[2];
684 p1 >> p2.ic1[0] >> p2.ic1[1] >> p2.ic2 >> p2.ic1[2];
685
686 p2.ic1[1] = M_PI * p2.ic1[1] / 180.0;
687 p2.ic1[2] = M_PI * p2.ic1[2] / 180.0;
688
689 p1 >> buffer;
690 p2.bt = bondtype(buffer[0]);
691
692 //////////////////////////////////////////////////////////////////////////////////////////
693 if (p2.bt.GetValue() == 0) { cout << "bad bondtype A" << endl; exit(-1); }
694 //////////////////////////////////////////////////////////////////////////////////////////
695
696 return p1;
697 }
698
operator <<(ostream & p1,mf_data_atm & p2)699 ostream & operator<<(ostream & p1, mf_data_atm & p2)
700 {
701 p1 << hex << "0x" << setw(4) << setfill('0') << p2.id << dec;
702 p1 << " " << p2.el.GetSymbol() << " " << p2.bt.GetSymbol1();
703 if (p2.tr != NULL) p1 << " " << (* p2.tr); p1 << " ";
704
705 p1 << p2.ic1[0] << " " << p2.ic1[1] << " " << p2.ic2 << " " << p2.ic1[2] << " ";
706 p1 << hex << "0x" << setw(4) << setfill('0') << p2.prev[0] << dec << " ";
707 p1 << hex << "0x" << setw(4) << setfill('0') << p2.prev[1] << dec << " ";
708 p1 << hex << "0x" << setw(4) << setfill('0') << p2.prev[2] << dec;
709
710 return p1;
711 }
712
713 /*################################################################################################*/
714
mf_data_bnd(void)715 mf_data_bnd::mf_data_bnd(void)
716 {
717 }
718
~mf_data_bnd(void)719 mf_data_bnd::~mf_data_bnd(void)
720 {
721 }
722
operator >>(istream & p1,mf_data_bnd & p2)723 istream & operator>>(istream & p1, mf_data_bnd & p2)
724 {
725 char buffer[256];
726 while (p1.get() != ':');
727 p1 >> p2.atm[0] >> p2.atm[1];
728 p1 >> buffer; p2.bt = bondtype(buffer[0]);
729
730 //////////////////////////////////////////////////////////////////////////////////////////
731 if (p2.bt.GetValue() == 0) { cout << "bad bondtype B" << endl; exit(-1); }
732 //////////////////////////////////////////////////////////////////////////////////////////
733
734 return p1;
735 }
736
737 /*################################################################################################*/
738
739 // eof
740