1 // MODEL.CPP
2
3 // Copyright (C) 1998 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 "model.h"
23
24 #include "v3d.h"
25
26 #include "seqbuild.h"
27 #include "resonance.h"
28
29 #include "tab_mm_default.h"
30
31 #include "eng1_qm.h"
32 #include "eng1_mm.h"
33 #include "eng1_sf.h"
34
35 #include "eng2_qm_mm.h"
36
37 #include "geomopt.h"
38 #include "moldyn.h"
39 #include "search.h"
40
41 #include "utility.h"
42
43 #include "tab_mm_default.h" // libghemical_init()
44 #include "tab_mm_tripos52.h" // libghemical_init()
45
46 #include "local_i18n.h"
47 #include "notice.h"
48
49 #include <stdio.h>
50
51 #include <cstring>
52 #include <iomanip>
53 #include <sstream>
54 using namespace std;
55
56 #ifdef WIN32
57 #include <time.h>
58 #endif // WIN32
59
60 /*################################################################################################*/
61
62 const char model::libversion[16] = LIBVERSION;
63 char model::libdata_path[256] = "libdata_path_is_not_yet_set_by_libghemical_init";
64
65 // 2007-01-16 ; previously these were just static members of the class,
66 // but now the header files model.h and seqbuild.h are becoming too dependent
67 // on each other ; this way the dependency can be dropped from model.h header.
68
69 sequencebuilder * model::amino_builder = NULL;
70 //singleton_cleaner<sequencebuilder> amino_cleaner(new sequencebuilder(chn_info::amino_acid, AMINO_BUILDER_FILE));
71
72 sequencebuilder * model::nucleic_builder = NULL;
73 //singleton_cleaner<sequencebuilder> nucleic_cleaner(new sequencebuilder(chn_info::nucleic_acid, NUCLEIC_BUILDER_FILE));
74
75 // sf_symbols and sf_types are in the same order.
76
77 const char model::sf_symbols[20 + 1] = "ARNDCQEGHILKMFPSTWYV";
78
79 const i32s model::sf_types[SF_NUM_TYPES] =
80 {
81 0x0000, // A
82 0x0100, 0x0101, 0x0102, // R
83 0x0200, 0x0201, // N
84 0x0300, 0x0301, // D
85 0x0400, 0x0401, // C
86 0x0500, 0x0501, // Q
87 0x0600, 0x0601, // E
88 0x0700, // G
89 0x0800, 0x0801, // H
90 0x0900, 0x0901, // I
91 0x0A00, 0x0A01, // L
92 0x0B00, 0x0B01, 0x0B02, // K
93 0x0C00, 0x0C01, // M
94 0x0D00, 0x0D01, // F
95 0x0E00, // P
96 0x0F00, // S
97 0x1000, // T
98 0x1100, 0x1101, 0x1102, // W
99 0x1200, 0x1201, // Y
100 0x1300 // V
101 };
102
103 const bool model::sf_is_polar[SF_NUM_TYPES] =
104 {
105 true, // A
106 true, false, true, // R
107 true, true, // N
108 true, true, // D
109 true, false, // C
110 true, true, // Q
111 true, true, // E
112 true, // G
113 true, true, // H
114 true, false, // I
115 true, false, // L
116 true, false, true, // K
117 true, false, // M
118 true, false, // F
119 true, // P
120 true, // S
121 true, // T
122 true, true, false, // W
123 true, true, // Y
124 false // V
125 };
126
model(void)127 model::model(void)
128 {
129 current_setup = new setup1_mm(this);
130
131 rs = NULL;
132
133 crd_table_size_glob = 1;
134 cs_vector.push_back(new crd_set());
135 SetCRDSetVisible(0, true);
136
137 is_index_clean = false;
138 is_groups_clean = false;
139 is_groups_sorted = false;
140
141 qm_total_charge = 0;
142 qm_current_orbital = 0;
143
144 use_boundary_potential = false;
145 saved_boundary_potential_rad_solute = 1.0;
146 saved_boundary_potential_rad_solvent = 1.0;
147
148 use_periodic_boundary_conditions = false;
149 saved_periodic_box_HALFdim[0] = 1.0;
150 saved_periodic_box_HALFdim[1] = 1.0;
151 saved_periodic_box_HALFdim[2] = 1.0;
152
153 nmol = NOT_DEFINED;
154 ref_civ = NULL;
155
156 trajfile = NULL;
157 traj_num_atoms = NOT_DEFINED;
158 total_traj_frames = NOT_DEFINED;
159 current_traj_frame = NOT_DEFINED;
160
161 verbosity = 3;
162
163 // initialize the ecomp stuff:
164 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^
165
166 ecomp_enabled = false;
167
168 const int defind = ecomp_AddGroup("default");
169 if (defind != 0) assertion_failed(__FILE__, __LINE__, "failed to initialize ecomp_grp_names.");
170 }
171
~model(void)172 model::~model(void)
173 {
174 if (current_setup != NULL)
175 {
176 delete current_setup;
177 current_setup = NULL;
178 }
179 else assertion_failed(__FILE__, __LINE__, "current_setup == NULL");
180
181 if (rs != NULL) delete rs;
182
183 for (i32u n1 = 0;n1 < cs_vector.size();n1++)
184 {
185 delete cs_vector[n1];
186 cs_vector[n1] = NULL;
187 }
188
189 if (trajfile != NULL) WarningMessage(_("WARNING : trajectory file was not closed!"));
190
191 for (i32u n1 = 0;n1 < ecomp_grp_names.size();n1++)
192 {
193 delete[] ecomp_grp_names[n1];
194 ecomp_grp_names[n1] = NULL;
195 }
196
197 if (ref_civ != NULL)
198 {
199 delete ref_civ;
200 ref_civ = NULL;
201 }
202 }
203
ThreadLock(void)204 void model::ThreadLock(void)
205 {
206 }
207
ThreadUnlock(void)208 void model::ThreadUnlock(void)
209 {
210 }
211
NoThreadsIterate(void)212 void model::NoThreadsIterate(void)
213 {
214 }
215
SetProgress(double,double *)216 bool model::SetProgress(double, double *)
217 {
218 return false;
219 }
220
Message(const char * msg)221 void model::Message(const char * msg)
222 {
223 cerr << msg << endl;
224 }
225
WarningMessage(const char * msg)226 void model::WarningMessage(const char * msg)
227 {
228 cerr << msg << endl;
229 }
230
ErrorMessage(const char * msg)231 void model::ErrorMessage(const char * msg)
232 {
233 cerr << msg << endl;
234 }
235
GetCurrentSetup(void)236 setup * model::GetCurrentSetup(void)
237 {
238 return current_setup;
239 }
240
ReplaceCurrentSetup(setup * su)241 void model::ReplaceCurrentSetup(setup * su)
242 {
243 if (su->GetModel() != this) assertion_failed(__FILE__, __LINE__, "bad setup passed as parameter.");
244 if (current_setup == NULL) assertion_failed(__FILE__, __LINE__, "current_setup == NULL");
245
246 delete current_setup;
247 current_setup = su;
248 }
249
GetRS(void)250 resonance_structures * model::GetRS(void)
251 {
252 return rs;
253 }
254
CreateRS(void)255 void model::CreateRS(void)
256 {
257 if (rs == NULL) rs = new resonance_structures(this);
258 }
259
DestroyRS(void)260 void model::DestroyRS(void)
261 {
262 if (rs != NULL)
263 {
264 delete rs;
265 rs = NULL;
266 }
267 }
268
OpenLibDataFile(ifstream & file,bool is_binary_file,const char * fn)269 void model::OpenLibDataFile(ifstream & file, bool is_binary_file, const char * fn)
270 {
271 ostringstream oss;
272
273 oss << model::libdata_path << DIR_SEPARATOR;
274 oss << model::libversion << DIR_SEPARATOR;
275 oss << fn << ends;
276
277 cout << _("DEBUG ; preparing to open file ") << oss.str() << endl;
278
279 if (!is_binary_file) file.open(oss.str().c_str(), ios::in);
280 else file.open(oss.str().c_str(), ios::in | ios::binary);
281
282 if (file.good()) return;
283
284 file.close();
285
286 cerr << _("ERROR : could not open data file : ") << oss.str().c_str() << endl;
287 cerr << _("The program will now exit. This file must be installed with this program.") << endl;
288 cerr << _("Re-installing the program and all the data files may solve this problem.") << endl;
289
290 exit(EXIT_FAILURE);
291 }
292
DiscardCurrEng(void)293 void model::DiscardCurrEng(void)
294 {
295 // cout << "discard!!!" << endl;
296 GetCurrentSetup()->DiscardCurrentEngine();
297 }
298
SetupPlotting(void)299 void model::SetupPlotting(void)
300 {
301 engine * eng = GetCurrentSetup()->GetCurrentEngine();
302 if (eng != NULL) eng->SetupPlotting();
303 }
304
305 /*##############################################*/
306 /*##############################################*/
307
GetCRDSetCount(void)308 i32u model::GetCRDSetCount(void)
309 {
310 return cs_vector.size();
311 }
312
GetCRDSetVisible(i32u index)313 bool model::GetCRDSetVisible(i32u index)
314 {
315 if (index < cs_vector.size())
316 {
317 return cs_vector[index]->visible;
318 }
319 else assertion_failed(__FILE__, __LINE__, "index overflow");
320 }
321
SetCRDSetVisible(i32u index,bool visible)322 void model::SetCRDSetVisible(i32u index, bool visible)
323 {
324 if (index < cs_vector.size())
325 {
326 cs_vector[index]->visible = visible;
327 }
328 else assertion_failed(__FILE__, __LINE__, "index overflow");
329 }
330
PushCRDSets(i32u p1)331 void model::PushCRDSets(i32u p1)
332 {
333 i32u old_size = cs_vector.size();
334
335 for (i32u n1 = 0;n1 < p1;n1++) cs_vector.push_back(new crd_set());
336
337 i32u new_size = cs_vector.size();
338
339 // determine whether we need to reallocate the crd_tables, and reallocate if needed.
340 // determine whether we need to reallocate the crd_tables, and reallocate if needed.
341 // determine whether we need to reallocate the crd_tables, and reallocate if needed.
342
343 if (new_size > crd_table_size_glob)
344 {
345 // i32u old_cap = crd_table_size_glob; this is actually not needed...
346 i32u new_cap;
347
348 new_cap = new_size; // just allocate the minimum amount (simple).
349
350 crd_table_size_glob = new_cap;
351 fGL * buff = new fGL[old_size * 3];
352
353 for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
354 {
355 for (i32u n1 = 0;n1 < old_size;n1++)
356 {
357 buff[n1 * 3 + 0] = (* it1).crd_table[n1 * 3 + 0];
358 buff[n1 * 3 + 1] = (* it1).crd_table[n1 * 3 + 1];
359 buff[n1 * 3 + 2] = (* it1).crd_table[n1 * 3 + 2];
360 }
361
362 delete[] (* it1).crd_table;
363 (* it1).crd_table = new fGL[new_cap * 3];
364 (* it1).crd_table_size_loc = new_cap;
365
366 for (i32u n1 = 0;n1 < old_size;n1++)
367 {
368 (* it1).crd_table[n1 * 3 + 0] = buff[n1 * 3 + 0];
369 (* it1).crd_table[n1 * 3 + 1] = buff[n1 * 3 + 1];
370 (* it1).crd_table[n1 * 3 + 2] = buff[n1 * 3 + 2];
371 }
372 }
373
374 delete[] buff;
375 }
376
377 // initialize the new memory blocks.
378 // initialize the new memory blocks.
379 // initialize the new memory blocks.
380
381 for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
382 {
383 for (i32u n1 = old_size;n1 < new_size;n1++)
384 {
385 (* it1).crd_table[n1 * 3 + 0] = 0.0;
386 (* it1).crd_table[n1 * 3 + 1] = 0.0;
387 (* it1).crd_table[n1 * 3 + 2] = 0.0;
388 }
389 }
390 }
391
PopCRDSets(i32u p1)392 void model::PopCRDSets(i32u p1)
393 {
394 // do not deallocate the memory; in the future deallocation could be added... 2003-06-12 TH
395 // do not deallocate the memory; in the future deallocation could be added... 2003-06-12 TH
396 // do not deallocate the memory; in the future deallocation could be added... 2003-06-12 TH
397
398 for (i32u n1 = 0;n1 < p1;n1++)
399 {
400 delete cs_vector.back();
401 cs_vector.pop_back();
402 }
403 }
404
CopyCRDSet(i32u p1,i32u p2)405 void model::CopyCRDSet(i32u p1, i32u p2)
406 {
407 if (p1 >= crd_table_size_glob || p2 >= crd_table_size_glob) assertion_failed(__FILE__, __LINE__, "cs overflow");
408
409 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
410 {
411 (* it1).crd_table[p2 * 3 + 0] = (* it1).crd_table[p1 * 3 + 0];
412 (* it1).crd_table[p2 * 3 + 1] = (* it1).crd_table[p1 * 3 + 1];
413 (* it1).crd_table[p2 * 3 + 2] = (* it1).crd_table[p1 * 3 + 2];
414 }
415 }
416
SwapCRDSets(i32u,i32u)417 void model::SwapCRDSets(i32u, i32u)
418 {
419 assertion_failed(__FILE__, __LINE__, "the method is not yet implemented!");
420 }
421
CenterCRDSet(i32u p1,bool all_atoms)422 void model::CenterCRDSet(i32u p1, bool all_atoms)
423 {
424 if (p1 >= crd_table_size_glob) assertion_failed(__FILE__, __LINE__, "cs overflow");
425
426 fGL sum[3] = { 0.0, 0.0, 0.0 };
427 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
428 {
429 if (!all_atoms && (* it1).flags & ATOMFLAG_IS_HIDDEN) continue;
430
431 fGL * crd_table = (* it1).crd_table;
432 for (i32s n1 = 0;n1 < 3;n1++) sum[n1] += crd_table[p1 * 3 + n1];
433 }
434
435 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
436 {
437 if (!all_atoms && (* it1).flags & ATOMFLAG_IS_HIDDEN) continue;
438
439 fGL * crd_table = (* it1).crd_table;
440 for (i32s n1 = 0;n1 < 3;n1++) crd_table[p1 * 3 + n1] -= sum[n1] / (fGL) GetAtomCount();
441 }
442 }
443
OrientCRDSet(i32u p1,bool all_atoms,fGL * array3)444 void model::OrientCRDSet(i32u p1, bool all_atoms, fGL * array3)
445 {
446 if (p1 >= crd_table_size_glob) assertion_failed(__FILE__, __LINE__, "cs overflow");
447
448 const fGL origo[3] = { 0.0, 0.0, 0.0 };
449
450 v3d<fGL> v1; v3d<fGL> v2;
451 v3d<fGL> vA; v3d<fGL> vB;
452 v3d<fGL> vv; v3d<fGL> vP; v3d<fGL> vT;
453
454 fGL maxd; fGL maxp[3];
455
456 // first orient to the X-axis (the more important step).
457 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
458 maxd = 0.0; maxp[0] = maxp[1] = maxp[2] = 0.0;
459 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
460 {
461 if (!all_atoms && (* it1).flags & ATOMFLAG_IS_HIDDEN) continue;
462
463 fGL * crd_table = (* it1).crd_table;
464
465 fGL tmp; fGL d = 0.0;
466 tmp = crd_table[p1 * 3 + 0]; d += tmp * tmp;
467 tmp = crd_table[p1 * 3 + 1]; d += tmp * tmp;
468 tmp = crd_table[p1 * 3 + 2]; d += tmp * tmp;
469 d = sqrt(d);
470
471 if (d > maxd)
472 {
473 maxd = d;
474 maxp[0] = crd_table[p1 * 3 + 0];
475 maxp[1] = crd_table[p1 * 3 + 1];
476 maxp[2] = crd_table[p1 * 3 + 2];
477 }
478 }
479
480 v1 = v3d<fGL>(origo, maxp); v1 = v1 / v1.len(); // initial...
481 vA = v3d<fGL>(+1.0, 0.0, 0.0); // final...
482
483 if (v1.ang(vA) > 0.1 * M_PI / 180.0)
484 {
485 vv = v1.vpr(vA);
486 vv = vv / vv.len();
487
488 v2 = v1.vpr(vv); v2 = v2 / v2.len(); // initial...
489 vB = vA.vpr(vv); vB = vB / vB.len(); // final...
490
491 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
492 {
493 if (!all_atoms && (* it1).flags & ATOMFLAG_IS_HIDDEN) continue;
494
495 fGL * crd_table = (* it1).crd_table;
496
497 vP = v3d<fGL>(& crd_table[p1 * 3]);
498
499 const fGL proj1 = vP.spr(v1);
500 const fGL proj2 = vP.spr(v2);
501 const fGL projV = vP.spr(vv);
502
503 vT = (vA * proj1) + (vB * proj2) + (vv * projV);
504
505 crd_table[p1 * 3 + 0] = vT[0];
506 crd_table[p1 * 3 + 1] = vT[1];
507 crd_table[p1 * 3 + 2] = vT[2];
508
509 // check that this code preserves chirality:
510 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
511 // ok ; 20060310 TH
512 }
513 }
514 else
515 {
516 ostringstream str;
517 str << _("Skipped stage 1 of Orient.") << endl << ends;
518
519 PrintToLog(str.str().c_str());
520 }
521
522 // then orient to the Y-axis (the less important step).
523 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
524 maxd = 0.0; maxp[0] = maxp[1] = maxp[2] = 0.0;
525 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
526 {
527 if (!all_atoms && (* it1).flags & ATOMFLAG_IS_HIDDEN) continue;
528
529 fGL * crd_table = (* it1).crd_table;
530
531 fGL tmp; fGL d = 0.0;
532 // skip the X-component of distance here!!!
533 tmp = crd_table[p1 * 3 + 1]; d += tmp * tmp;
534 tmp = crd_table[p1 * 3 + 2]; d += tmp * tmp;
535 d = sqrt(d);
536
537 if (d > maxd)
538 {
539 maxd = d;
540 maxp[0] = crd_table[p1 * 3 + 0];
541 maxp[1] = crd_table[p1 * 3 + 1];
542 maxp[2] = crd_table[p1 * 3 + 2];
543 }
544 }
545
546 fGL refpt[3] = { maxp[0], 0.0, 0.0 };
547
548 v1 = v3d<fGL>(refpt, maxp); v1 = v1 / v1.len(); // initial...
549 vA = v3d<fGL>(0.0, +1.0, 0.0); // final...
550
551 if (v1.ang(vA) > 0.1 * M_PI / 180.0)
552 {
553 vv = v1.vpr(vA);
554 vv = vv / vv.len();
555
556 v2 = v1.vpr(vv); v2 = v2 / v2.len(); // initial...
557 vB = vA.vpr(vv); vB = vB / vB.len(); // final...
558
559 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
560 {
561 if (!all_atoms && (* it1).flags & ATOMFLAG_IS_HIDDEN) continue;
562
563 fGL * crd_table = (* it1).crd_table;
564 refpt[0] = crd_table[p1 * 3 + 0];
565
566 vP = v3d<fGL>(refpt, & crd_table[p1 * 3]);
567 const fGL proj1 = vP.spr(v1);
568 const fGL proj2 = vP.spr(v2);
569 const fGL projV = vP.spr(vv);
570
571 vT = (vA * proj1) + (vB * proj2) + (vv * projV);
572
573 crd_table[p1 * 3 + 0] = refpt[0] + vT[0];
574 crd_table[p1 * 3 + 1] = refpt[1] + vT[1];
575 crd_table[p1 * 3 + 2] = refpt[2] + vT[2];
576
577 // check that this code preserves chirality:
578 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
579 // ok ; 20060310 TH
580 }
581 }
582 else
583 {
584 ostringstream str;
585 str << _("Skipped stage 2 of Orient.") << endl << ends;
586
587 PrintToLog(str.str().c_str());
588 }
589
590 // finally compute and return the max-dimensions...
591 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
592 array3[0] = array3[1] = array3[2] = 0.0;
593
594 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
595 {
596 if (!all_atoms && (* it1).flags & ATOMFLAG_IS_HIDDEN) continue;
597
598 fGL * crd_table = (* it1).crd_table;
599
600 fGL tmp;
601
602 tmp = fabs(crd_table[p1 * 3 + 0]);
603 if (tmp > array3[0]) array3[0] = tmp;
604
605 tmp = fabs(crd_table[p1 * 3 + 1]);
606 if (tmp > array3[1]) array3[1] = tmp;
607
608 tmp = fabs(crd_table[p1 * 3 + 2]);
609 if (tmp > array3[2]) array3[2] = tmp;
610 }
611 }
612
ReserveCRDSets(i32u)613 void model::ReserveCRDSets(i32u)
614 {
615 assertion_failed(__FILE__, __LINE__, "the method is not yet implemented!");
616 }
617
618 /*##############################################*/
619 /*##############################################*/
620
AddAtom_lg(atom & p1)621 void model::AddAtom_lg(atom & p1)
622 {
623 SystemWasModified();
624
625 //////////////////////////////////////////////////////////////////////////obsolete...
626 // p1.index = (i32s) atom_list.size(); // try to keep the atom::index records up-to-date...
627 //////////////////////////////////////////////////////////////////////////obsolete...
628
629 i32s added_index = (i32s) atom_list.size();
630 atom_list.push_back(p1);
631
632 atom_list.back().index = added_index; // try to keep the atom::index records up-to-date...
633 atom_list.back().mdl = this;
634 }
635
RemoveAtom(iter_al it1)636 void model::RemoveAtom(iter_al it1)
637 {
638 // todo : make sure that the iterator is valid for the list!
639 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
640
641 SystemWasModified();
642
643 // this strange while-loop is needed because removing bonds
644 // will invalidate all those atom::cr_list-iterators...
645
646 while (!(* it1).cr_list.empty())
647 {
648 crec * ref = & (* it1).cr_list.back();
649 iter_bl it2 = find(bond_list.begin(), bond_list.end(), (* ref->bndr));
650
651 if (it2 != bond_list.end()) RemoveBond(it2);
652 else assertion_failed(__FILE__, __LINE__, "find failed in bond_list.");
653 }
654
655 // remove any constraints that are impacted by the atom removal.
656
657 while(true)
658 {
659 iter_CDl it2 = FindAtomConstraint((* it1));
660 if (it2 == GetConstD_end()) break;
661 else RemoveConstraint(it2);
662 }
663
664 (* it1).mdl = NULL;
665
666 i32s removed_index = (* it1).index;
667 atom_list.erase(it1);
668
669 // try to keep the atom::index records up-to-date...
670
671 for (it1 = atom_list.begin();it1 != atom_list.end();it1++)
672 {
673 if ((* it1).index < removed_index) continue;
674 else (* it1).index--;
675 }
676 }
677
AddBond(bond & p1)678 void model::AddBond(bond & p1)
679 {
680 if (p1.atmr[0] == p1.atmr[1]) assertion_failed(__FILE__, __LINE__, "tried to add an invalid bond.");
681 if (p1.atmr[0]->mdl != p1.atmr[1]->mdl || p1.atmr[0]->mdl == NULL)assertion_failed(__FILE__, __LINE__, "tried to add an invalid bond.");
682
683 SystemWasModified();
684
685 bond_list.push_back(p1);
686
687 crec info1 = crec(p1.atmr[1], & bond_list.back());
688 p1.atmr[0]->cr_list.push_back(info1);
689 /* if (p1.atmr[0]->cr_list.size() > 8) // this is for debugging purposes only...
690 {
691 assertion_failed(__FILE__, __LINE__, "too many bonds!");
692 } */
693
694 crec info2 = crec(p1.atmr[0], & bond_list.back());
695 p1.atmr[1]->cr_list.push_back(info2);
696 /* if (p1.atmr[1]->cr_list.size() > 8) // this is for debugging purposes only...
697 {
698 assertion_failed(__FILE__, __LINE__, "too many bonds!");
699 } */
700 }
701
RemoveBond(iter_bl it1)702 void model::RemoveBond(iter_bl it1)
703 {
704 // todo : make sure that the iterator is valid for the list!
705 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
706
707 SystemWasModified();
708
709 crec tmpinfo = crec(NULL, & (* it1));
710 iter_cl it2;
711
712 it2 = find((* it1).atmr[0]->cr_list.begin(), (* it1).atmr[0]->cr_list.end(), tmpinfo);
713 if (it2 != (* it1).atmr[0]->cr_list.end()) (* it1).atmr[0]->cr_list.erase(it2);
714 else assertion_failed(__FILE__, __LINE__, "find failed in cr_list");
715
716 it2 = find((* it1).atmr[1]->cr_list.begin(), (* it1).atmr[1]->cr_list.end(), tmpinfo);
717 if (it2 != (* it1).atmr[1]->cr_list.end()) (* it1).atmr[1]->cr_list.erase(it2);
718 else assertion_failed(__FILE__, __LINE__, "find failed in cr_list");
719
720 bond_list.erase(it1);
721 }
722
AddConstraint(constraint_dst & p1)723 void model::AddConstraint(constraint_dst & p1)
724 {
725 bool same_atoms = (p1.atmr[0] == p1.atmr[1]);
726 bool null_pointer = (!p1.atmr[0] || !p1.atmr[1]);
727 bool bad_model = (p1.atmr[0]->mdl != this || p1.atmr[1]->mdl != this);
728
729 if (same_atoms || null_pointer || bad_model) assertion_failed(__FILE__, __LINE__, "bad constraint");
730
731 // first check whether the constraint already exists...
732 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
733
734 iter_CDl it1 = find(GetConstD_begin(), GetConstD_end(), p1);
735 if (it1 != GetConstD_end())
736 {
737 SystemWasModified(); // this is needed ; equivalent to adding a bond.
738
739 (* it1).SetType(p1.GetType());
740 (* it1).SetMinDist(p1.GetMinDist());
741 (* it1).SetMinFC(p1.GetMinFC());
742 (* it1).SetMaxDist(p1.GetMaxDist());
743 (* it1).SetMaxFC(p1.GetMaxFC());
744
745 return;
746 }
747
748 // ...and if it doesn't then add it.
749 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
750
751 SystemWasModified(); // this is needed ; equivalent to adding a bond.
752
753 const_D_list.push_back(p1);
754 }
755
RemoveConstraint(iter_CDl it1)756 void model::RemoveConstraint(iter_CDl it1)
757 {
758 // todo : make sure that the iterator is valid for the list!
759 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
760
761 SystemWasModified(); // this is needed ; equivalent to a bond removal.
762
763 const_D_list.erase(it1);
764
765 return;
766 }
767
768 // this should be called ALWAYS if ANY modification is done to the system.
769 // automatically called by Add/Remove/Atom/Bond. GUI should change if change in element etc...
770
SystemWasModified(void)771 void model::SystemWasModified(void)
772 {
773 // if atoms/bonds added/removed, all engine-objects have to be discarded.
774 // also if an element is changed, the engine objects (at least) must be discarded.
775
776 DiscardCurrEng();
777
778 // in the setup object, all atom/bond tables must be discarded for the same reasons as above...
779
780 current_setup->DiscardSetupInfo();
781
782 // the resonance_structures object must be invalidated, if there is one...
783
784 if (rs != NULL)
785 {
786 delete rs;
787 rs = NULL;
788 }
789
790 // any change that might lead to change in detected sequences must invalidate the SF engine objects...
791
792 if (dynamic_cast<setup1_sf *>(current_setup) != NULL)
793 {
794 // delete current_setup;
795 // current_setup = new setup1_mm(this);
796 ReplaceCurrentSetup(new setup1_mm(this));
797 }
798
799 // any addition/removal of atoms/bonds will invalidate information about molecules.
800
801 InvalidateGroups();
802 }
803
FindAtomByIndex(i32s index)804 iter_al model::FindAtomByIndex(i32s index)
805 {
806 iter_al end_it = GetAtomsEnd();
807 if (index < 0) return end_it;
808
809 i32s counter = 0;
810 iter_al iter = GetAtomsBegin();
811 while (counter != index)
812 {
813 counter++; iter++;
814 if (iter == end_it) return end_it;
815 }
816
817 return iter;
818 }
819
FindAtomConstraint(atom & p1)820 iter_CDl model::FindAtomConstraint(atom & p1)
821 {
822 for (iter_CDl it1 = GetConstD_begin();it1 != GetConstD_end();it1++)
823 {
824 if((* it1).atmr[0] == & p1 || (* it1).atmr[1] == & p1)
825 {
826 return it1;
827 }
828 }
829
830 return GetConstD_end();
831 }
832
ClearModel(void)833 void model::ClearModel(void)
834 {
835 while (!bond_list.empty())
836 {
837 iter_bl it1 = bond_list.begin();
838 RemoveBond(it1);
839 }
840
841 while (!atom_list.empty())
842 {
843 iter_al it1 = atom_list.begin();
844 RemoveAtom(it1);
845 }
846
847 // do this last as ideally removing the atoms
848 // should get all the constraints as well....
849 while (!const_D_list.empty())
850 {
851 Message("DEBUG_WARNING : constr_D_list was not empty!");
852 cout << "DEBUG_WARNING : constr_D_list was not empty!" << endl;
853
854 iter_CDl it1 = const_D_list.begin();
855 RemoveConstraint(it1);
856 }
857 }
858
859 /*##############################################*/
860 /*##############################################*/
861
GetRange(i32s ind,i32s value,iter_al * result)862 void model::GetRange(i32s ind, i32s value, iter_al * result)
863 {
864 iter_al range[2] = { atom_list.begin(), atom_list.end() };
865 GetRange(ind, range, value, result); // call GetRange() using full atom_list!!!
866 }
867
GetRange(i32s ind,iter_al * range,i32s value,iter_al * result)868 void model::GetRange(i32s ind, iter_al * range, i32s value, iter_al * result)
869 {
870 if (!is_groups_sorted) assertion_failed(__FILE__, __LINE__, "!is_groups_sorted");
871
872 result[0] = range[0]; while (result[0] != range[1] && (* result[0]).id[ind] != value) result[0]++;
873 result[1] = result[0]; while (result[1] != range[1] && (* result[1]).id[ind] == value) result[1]++;
874 }
875
GetRange(i32s molecule,iter_bl * result)876 void model::GetRange(i32s molecule, iter_bl * result)
877 {
878 if (!is_groups_sorted) assertion_failed(__FILE__, __LINE__, "!is_groups_sorted");
879
880 // assume that the molecule numbers of both atoms in a bond object are identical!!!
881 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
882
883 result[0] = bond_list.begin(); while (result[0] != bond_list.end() && (* result[0]).atmr[0]->id[0] != molecule) result[0]++;
884 result[1] = result[0]; while (result[1] != bond_list.end() && (* result[1]).atmr[0]->id[0] == molecule) result[1]++;
885 }
886
FindPath(atom * ref1,atom * ref2,i32s max,i32s flag,i32s dist)887 i32s model::FindPath(atom * ref1, atom * ref2, i32s max, i32s flag, i32s dist)
888 {
889 if (ref1 == ref2) return dist;
890 if (dist == max) return NOT_FOUND;
891
892 i32s tmp1 = NOT_FOUND; iter_cl it1;
893 for (it1 = ref1->cr_list.begin();it1 != ref1->cr_list.end();it1++)
894 {
895 if ((* it1).bndr->flags[flag]) continue;
896
897 (* it1).bndr->flags[flag] = true;
898 i32s tmp2 = FindPath((* it1).atmr, ref2, max, flag, dist + 1);
899 (* it1).bndr->flags[flag] = false;
900
901 if (tmp2 < tmp1) tmp1 = tmp2;
902 }
903
904 return tmp1;
905 }
906
FindPathV(atom * ref1,atom * ref2,i32s max,i32s flag,i32s dist)907 vector<bond *> * model::FindPathV(atom * ref1, atom * ref2, i32s max, i32s flag, i32s dist)
908 {
909 if (ref1 == ref2) return new vector<bond *>;
910 if (dist == max) return NULL;
911
912 vector<bond *> * tmp1 = NULL; iter_cl it1;
913 for (it1 = ref1->cr_list.begin();it1 != ref1->cr_list.end();it1++)
914 {
915 if ((* it1).bndr->flags[flag]) continue;
916
917 (* it1).bndr->flags[flag] = true;
918 vector<bond *> * tmp2 = FindPathV((* it1).atmr, ref2, max, flag, dist + 1);
919 (* it1).bndr->flags[flag] = false;
920
921 if (tmp2 != NULL)
922 {
923 tmp2->push_back((* it1).bndr);
924 if (tmp1 == NULL) tmp1 = tmp2;
925 else
926 {
927 if (tmp2->size() < tmp1->size())
928 {
929 delete tmp1;
930 tmp1 = tmp2;
931 }
932 }
933 }
934 }
935
936 return tmp1;
937 }
938
FindRing(atom * ref1,atom * ref2,signed char * str,i32s size,i32s flag,i32s dist)939 bool model::FindRing(atom * ref1, atom * ref2, signed char * str, i32s size, i32s flag, i32s dist)
940 {
941 static vector<signed char> ring_vector;
942
943 if (!dist && str != NULL) ring_vector.resize(0);
944 else if (dist && ref1 == ref2)
945 {
946 if (dist != size) return false;
947
948 if (str != NULL)
949 {
950 for (i32u n1 = 0;n1 < strlen((const char *) str);n1++)
951 {
952 if (!(n1 % 2) && str[n1] == '?') continue;
953 if ((n1 % 2) && str[n1] == NOT_DEFINED) continue;
954 if (str[n1] != ring_vector[n1]) return false;
955 }
956 }
957
958 return true;
959 }
960
961 if (dist == size) return false; iter_cl it1;
962 for (it1 = ref1->cr_list.begin();it1 != ref1->cr_list.end();it1++)
963 {
964 if ((* it1).bndr->flags[flag]) continue;
965
966 if (str != NULL)
967 {
968 ring_vector.push_back((* it1).bndr->bt.GetSymbol2());
969 ring_vector.push_back((signed char) (* it1).atmr->el.GetAtomicNumber());
970 }
971
972 (* it1).bndr->flags[flag] = true;
973 bool result = FindRing((* it1).atmr, ref2, str, size, flag, dist + 1);
974 (* it1).bndr->flags[flag] = false;
975
976 if (str != NULL)
977 {
978 ring_vector.pop_back();
979 ring_vector.pop_back();
980 }
981
982 if (result) return true;
983 }
984
985 return false;
986 }
987
InvalidateGroups(void)988 void model::InvalidateGroups(void)
989 {
990 is_index_clean = false;
991 is_groups_clean = false;
992 is_groups_sorted = false;
993
994 nmol = NOT_DEFINED;
995 if (ref_civ != NULL)
996 {
997 delete ref_civ;
998 ref_civ = NULL;
999 }
1000 }
1001
UpdateIndex(void)1002 void model::UpdateIndex(void)
1003 {
1004 i32s counter = 0;
1005 iter_al it1 = atom_list.begin();
1006 while (it1 != atom_list.end())
1007 {
1008 (* it1++).index = counter++;
1009 }
1010
1011 is_index_clean = true;
1012 }
1013
UpdateGroups(void)1014 void model::UpdateGroups(void)
1015 {
1016 InvalidateGroups();
1017 UpdateIndex(); // looks foolish, but this is a quick operation...
1018
1019 nmol = 0;
1020
1021 iter_al it1;
1022 for (it1 = atom_list.begin();it1 != atom_list.end();it1++)
1023 {
1024 (* it1).id[0] = (* it1).id[1] = (* it1).id[2] = (* it1).id[3] = NOT_DEFINED;
1025 }
1026
1027 while (true)
1028 {
1029 for (it1 = atom_list.begin();it1 != atom_list.end() && (* it1).id[0] != NOT_DEFINED;it1++);
1030 if (it1 != atom_list.end()) GatherAtoms(& (* it1), nmol++); else break;
1031 }
1032
1033 UpdateIndex(); // this is so that we can keep is_index_clean true!
1034
1035 is_groups_clean = true;
1036 }
1037
SortGroups(void)1038 void model::SortGroups(void)
1039 {
1040 if (!is_groups_clean) assertion_failed(__FILE__, __LINE__, "!is_groups_clean");
1041
1042 // sorting the atom list to get contiguous molecules/chains/residues is an efficient technique, but may be
1043 // confusing in the cases where atom indexing must not change (like TSS); therefore always inform the user!!!
1044
1045 ostringstream str;
1046 str << _("Calling model::SortGroups() so the atom indexing may change!") << endl << ends;
1047 if (verbosity >= 3) PrintToLog(str.str().c_str());
1048
1049 atom_list.sort(); // this should be the ONLY place where atom_list is sorted!!!
1050 UpdateIndex(); // this is so that we can keep is_index_clean true!
1051
1052 bond_list.sort(); // this should be the ONLY place where bond_list is sorted!!!
1053
1054 // sorting the atom_list may cause that model::atom_list and setup::atmtab are ordered differently; synchronize!!!
1055 // sorting the atom_list may cause that model::atom_list and setup::atmtab are ordered differently; synchronize!!!
1056 // sorting the atom_list may cause that model::atom_list and setup::atmtab are ordered differently; synchronize!!!
1057 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1058 //GetCurrentSetup()->UpdateSetupInfo(); // not needed? called in every setup::CreateEngineByIndex() implementation...
1059
1060 is_groups_sorted = true;
1061 }
1062
UpdateChains(bool skip_nucleic)1063 void model::UpdateChains(bool skip_nucleic)
1064 {
1065 if (!is_groups_clean) UpdateGroups();
1066 if (!is_groups_sorted) SortGroups(); // at sequencebuilder::Identify() we will need model::GetRange().
1067
1068 if (ref_civ != NULL) delete ref_civ;
1069 ref_civ = new vector<chn_info>;
1070
1071 model::amino_builder->Identify(this); // do this 1st!!!
1072 if (!skip_nucleic) model::nucleic_builder->Identify(this); // do this 2nd!!!
1073
1074 // now the protein chains are defined first,
1075 // and DNA/RNA chains later (if not denied)...
1076
1077 SortGroups(); // ok, then sort the atom_list using identified chains/residues...
1078 }
1079
GatherAtoms(atom * ref,i32s id)1080 void model::GatherAtoms(atom * ref, i32s id)
1081 {
1082 if (ref->id[0] != NOT_DEFINED) return;
1083 ref->id[0] = id; iter_cl it1 = ref->cr_list.begin();
1084 while (it1 != ref->cr_list.end()) GatherAtoms((* it1++).atmr, id);
1085 }
1086
1087 /*##############################################*/
1088 /*##############################################*/
1089
cp_FindAtom(iter_al * res_rng,i32s id)1090 atom * model::cp_FindAtom(iter_al * res_rng, i32s id) // see default_tables::e_UT_FindAtom()...
1091 {
1092 iter_al it1 = res_rng[0];
1093 while (it1 != res_rng[1] && ((* it1).builder_res_id & 0xFF) != id) it1++;
1094 if (it1 == res_rng[1]) return NULL; else return & (* it1);
1095 }
1096
CheckProtonation(void)1097 void model::CheckProtonation(void)
1098 {
1099 if (ref_civ == NULL) UpdateChains();
1100
1101 // todo ; use the newer "sequence3" information as well???
1102 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1103 // or at least check whether it is consistent or not...
1104
1105 vector<chn_info> & ci_vector = (* ref_civ);
1106 for (i32u n1 = 0;n1 < ci_vector.size();n1++)
1107 {
1108 iter_al chnR[2]; GetRange(1, n1, chnR);
1109
1110 // amino...
1111 // amino...
1112 // amino...
1113
1114 if (ci_vector[n1].GetType() == chn_info::amino_acid)
1115 {
1116 const char * sequence1 = ci_vector[n1].GetSequence1();
1117 const char * pstate = ci_vector[n1].GetProtonationStates();
1118
1119 if (pstate != NULL)
1120 {
1121 ostringstream str;
1122 str << _("CheckProtonation() : pstate array found for chain ") << n1 << "." << endl << ends;
1123 PrintToLog(str.str().c_str());
1124 }
1125 else
1126 {
1127 ostringstream str;
1128 str << _("CheckProtonation() : no pstate array found for chain ") << n1 << _("; USING DEFAULTS!") << endl << ends;
1129 PrintToLog(str.str().c_str());
1130 }
1131
1132 for (i32s n2 = 0;n2 < ci_vector[n1].GetLength();n2++)
1133 {
1134 iter_al resR[2]; GetRange(2, chnR, n2, resR);
1135
1136 atom * atmr_00 = cp_FindAtom(resR, 0x00); // N
1137 atom * atmr_01 = cp_FindAtom(resR, 0x01); // alpha-C
1138 atom * atmr_02 = cp_FindAtom(resR, 0x02); // carbonyl-C
1139 atom * atmr_10 = cp_FindAtom(resR, 0x10); // carbonyl-O
1140
1141 if (!atmr_00 || !atmr_01 || !atmr_02 || !atmr_10) continue;
1142
1143 int sidechain_charge = 0;
1144 if (pstate != NULL)
1145 {
1146 sidechain_charge = (pstate[n2] & PSTATE_CHARGE_mask);
1147 if (pstate[n2] & PSTATE_SIGN_NEGATIVE) sidechain_charge = -sidechain_charge;
1148 else if (!(pstate[n2] & PSTATE_SIGN_POSITIVE))
1149 {
1150 assertion_failed(__FILE__, __LINE__, "no sign information found from pstate!");
1151 }
1152 }
1153
1154 // cout << "sidechain_charge = " << sidechain_charge << endl; // debug...
1155
1156 if (n2 == 0) // N-terminal residue???
1157 {
1158 bool charged = true;
1159 if (pstate != NULL) charged = ((pstate[n2] & PSTATE_CHARGED_TERMINAL) ? true : false);
1160
1161 ostringstream str;
1162 str << _("CheckProtonation() : setting N-terminal ") << (charged ? _("charged.") : _("neutral.")) << endl << ends;
1163 PrintToLog(str.str().c_str());
1164
1165 if (charged)
1166 {
1167 atmr_00->formal_charge = +1;
1168
1169 sidechain_charge -= 1;
1170 }
1171 else
1172 {
1173 atmr_00->formal_charge = +0;
1174 }
1175 }
1176
1177 if (n2 == ci_vector[n1].GetLength() - 1) // C-terminal residue???
1178 {
1179 // atom * atmr_11 = cp_FindAtom(resR, 0x11);
1180 atom * atmr_11 = NULL; // search OXT directly since seq-builder misses it...
1181 for (iter_cl it1 = atmr_02->cr_list.begin();it1 != atmr_02->cr_list.end();it1++)
1182 {
1183 if ((* it1).atmr->el.GetAtomicNumber() != 8) continue;
1184 if (((* it1).atmr->builder_res_id & 0xFF) == 0x10) continue;
1185
1186 atmr_11 = (* it1).atmr;
1187 }
1188
1189 bool charged = true;
1190 if (pstate != NULL) charged = ((pstate[n2] & PSTATE_CHARGED_TERMINAL) ? true : false);
1191
1192 ostringstream str;
1193 str << _("CheckProtonation() : setting C-terminal ") << (charged ? _("charged.") : _("neutral.")) << endl << ends;
1194 PrintToLog(str.str().c_str());
1195
1196 bond tmpb1(atmr_02, atmr_10, bondtype('S'));
1197 iter_bl tmpi1 = find(bond_list.begin(), bond_list.end(), tmpb1);
1198
1199 bond tmpb2(atmr_02, atmr_11, bondtype('S'));
1200 iter_bl tmpi2 = find(bond_list.begin(), bond_list.end(), tmpb2);
1201
1202 if (charged)
1203 {
1204 atmr_10->formal_charge = +0;
1205 atmr_11->formal_charge = -1;
1206
1207 if (tmpi1 != bond_list.end()) (* tmpi1).bt = bondtype('D');
1208 if (tmpi2 != bond_list.end()) (* tmpi2).bt = bondtype('S');
1209
1210 sidechain_charge += 1;
1211 }
1212 else
1213 {
1214 atmr_10->formal_charge = +0;
1215 atmr_11->formal_charge = +0;
1216
1217 if (tmpi1 != bond_list.end()) (* tmpi1).bt = bondtype('D');
1218 if (tmpi2 != bond_list.end()) (* tmpi2).bt = bondtype('S');
1219 }
1220 }
1221
1222 // start handling the side chains!!!
1223 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1224
1225 /* if (sequence1[n2] == 'A')
1226 {
1227 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1228 } */
1229
1230 /* if (sequence1[n2] == 'R')
1231 {
1232 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1233 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1234 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1235 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1236 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1237 atom * atmr_25 = cp_FindAtom(resR, 0x25);
1238 atom * atmr_26 = cp_FindAtom(resR, 0x26);
1239 } */
1240
1241 /* if (sequence1[n2] == 'N')
1242 {
1243 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1244 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1245 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1246 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1247 } */
1248
1249 if (sequence1[n2] == 'D')
1250 {
1251 /* atom * atmr_20 = cp_FindAtom(resR, 0x20); */
1252 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1253 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1254 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1255
1256 bool charged = true;
1257 if (pstate != NULL)
1258 {
1259 switch (sidechain_charge)
1260 {
1261 case 0: charged = false; break;
1262 case -1: charged = true; break;
1263 default:
1264 ostringstream msg;
1265 msg << "bad charge : " << sidechain_charge << " is out of valid range for residue D (" << n1 << "/" << n2 << ")." << ends;
1266 assertion_failed(__FILE__, __LINE__, msg.str().c_str());
1267 }
1268 }
1269
1270 ostringstream str;
1271 str << _("CheckProtonation() : setting residue ") << n2 << " D " << (charged ? _("charged.") : _("neutral.")) << endl << ends;
1272 PrintToLog(str.str().c_str());
1273
1274 bond tmpb1(atmr_21, atmr_22, bondtype('S'));
1275 iter_bl tmpi1 = find(bond_list.begin(), bond_list.end(), tmpb1);
1276
1277 bond tmpb2(atmr_21, atmr_23, bondtype('S'));
1278 iter_bl tmpi2 = find(bond_list.begin(), bond_list.end(), tmpb2);
1279
1280 if (charged)
1281 {
1282 atmr_22->formal_charge = +0;
1283 atmr_23->formal_charge = -1;
1284
1285 if (tmpi1 != bond_list.end()) (* tmpi1).bt = bondtype('D');
1286 if (tmpi2 != bond_list.end()) (* tmpi2).bt = bondtype('S');
1287 }
1288 else
1289 {
1290 atmr_22->formal_charge = +0;
1291 atmr_23->formal_charge = +0;
1292
1293 if (tmpi1 != bond_list.end()) (* tmpi1).bt = bondtype('D');
1294 if (tmpi2 != bond_list.end()) (* tmpi2).bt = bondtype('S');
1295 }
1296 }
1297
1298 if (sequence1[n2] == 'C')
1299 {
1300 /* atom * atmr_20 = cp_FindAtom(resR, 0x20); */
1301 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1302
1303 bool charged = false;
1304 if (pstate != NULL)
1305 {
1306 switch (sidechain_charge)
1307 {
1308 case 0: charged = false; break;
1309 case -1: charged = true; break;
1310 default:
1311 ostringstream msg;
1312 msg << "bad charge : " << sidechain_charge << " is out of valid range for residue C (" << n1 << "/" << n2 << ")." << ends;
1313 assertion_failed(__FILE__, __LINE__, msg.str().c_str());
1314 }
1315 }
1316
1317 ostringstream str;
1318 str << _("CheckProtonation() : setting residue ") << n2 << " C " << (charged ? _("charged.") : _("neutral.")) << endl << ends;
1319 PrintToLog(str.str().c_str());
1320
1321 if (charged)
1322 {
1323 atmr_21->formal_charge = -1;
1324 }
1325 else
1326 {
1327 atmr_21->formal_charge = +0;
1328 }
1329 }
1330
1331 /* if (sequence1[n2] == 'Q')
1332 {
1333 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1334 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1335 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1336 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1337 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1338 } */
1339
1340 if (sequence1[n2] == 'E')
1341 {
1342 /* atom * atmr_20 = cp_FindAtom(resR, 0x20);
1343 atom * atmr_21 = cp_FindAtom(resR, 0x21); */
1344 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1345 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1346 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1347
1348 bool charged = true;
1349 if (pstate != NULL)
1350 {
1351 switch (sidechain_charge)
1352 {
1353 case 0: charged = false; break;
1354 case -1: charged = true; break;
1355 default:
1356 ostringstream msg;
1357 msg << "bad charge : " << sidechain_charge << " is out of valid range for residue E (" << n1 << "/" << n2 << ")." << ends;
1358 assertion_failed(__FILE__, __LINE__, msg.str().c_str());
1359 }
1360 }
1361
1362 ostringstream str;
1363 str << _("CheckProtonation() : setting residue ") << n2 << " E " << (charged ? _("charged.") : _("neutral.")) << endl << ends;
1364 PrintToLog(str.str().c_str());
1365
1366 bond tmpb1(atmr_22, atmr_23, bondtype('S'));
1367 iter_bl tmpi1 = find(bond_list.begin(), bond_list.end(), tmpb1);
1368
1369 bond tmpb2(atmr_22, atmr_24, bondtype('S'));
1370 iter_bl tmpi2 = find(bond_list.begin(), bond_list.end(), tmpb2);
1371
1372 if (charged)
1373 {
1374 atmr_23->formal_charge = +0;
1375 atmr_24->formal_charge = -1;
1376
1377 if (tmpi1 != bond_list.end()) (* tmpi1).bt = bondtype('D');
1378 if (tmpi2 != bond_list.end()) (* tmpi2).bt = bondtype('S');
1379 }
1380 else
1381 {
1382 atmr_23->formal_charge = +0;
1383 atmr_24->formal_charge = +0;
1384
1385 if (tmpi1 != bond_list.end()) (* tmpi1).bt = bondtype('D');
1386 if (tmpi2 != bond_list.end()) (* tmpi2).bt = bondtype('S');
1387 }
1388 }
1389
1390 /* if (sequence1[n2] == 'G')
1391 {
1392 } */
1393
1394 if (sequence1[n2] == 'H')
1395 {
1396 /* atom * atmr_20 = cp_FindAtom(resR, 0x20);
1397 atom * atmr_21 = cp_FindAtom(resR, 0x21); */
1398 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1399 /* atom * atmr_23 = cp_FindAtom(resR, 0x23); */
1400 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1401 /* atom * atmr_25 = cp_FindAtom(resR, 0x25); */
1402
1403 bool charged = false;
1404 if (pstate != NULL)
1405 {
1406 switch (sidechain_charge)
1407 {
1408 case 0: charged = false; break;
1409 case +1: charged = true; break;
1410 default:
1411 ostringstream msg;
1412 msg << "bad charge : " << sidechain_charge << " is out of valid range for residue H (" << n1 << "/" << n2 << ")." << ends;
1413 assertion_failed(__FILE__, __LINE__, msg.str().c_str());
1414 }
1415 }
1416
1417 // if neutral, should it be epsilon (HIE) or delta (HID)???
1418
1419 ostringstream str;
1420 str << _("CheckProtonation() : setting residue ") << n2 << " H " << (charged ? _("charged.") : _("neutral(HIE).")) << endl << ends;
1421 PrintToLog(str.str().c_str());
1422
1423 if (charged)
1424 {
1425 atmr_22->formal_charge = +0;//fixme
1426 atmr_24->formal_charge = +0;//fixme
1427
1428 // todo ; check the bonds!!!
1429 // todo ; check the bonds!!!
1430 }
1431 else
1432 {
1433 atmr_22->formal_charge = +0;//fixme
1434 atmr_24->formal_charge = +0;//fixme
1435
1436 // todo ; check the bonds!!!
1437 // todo ; check the bonds!!!
1438 }
1439 }
1440
1441 /* if (sequence1[n2] == 'I')
1442 {
1443 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1444 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1445 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1446 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1447 } */
1448
1449 /* if (sequence1[n2] == 'L')
1450 {
1451 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1452 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1453 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1454 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1455 } */
1456
1457 if (sequence1[n2] == 'K')
1458 {
1459 /* atom * atmr_20 = cp_FindAtom(resR, 0x20);
1460 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1461 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1462 atom * atmr_23 = cp_FindAtom(resR, 0x23); */
1463 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1464
1465 bool charged = true;
1466 if (pstate != NULL)
1467 {
1468 switch (sidechain_charge)
1469 {
1470 case 0: charged = false; break;
1471 case +1: charged = true; break;
1472 default:
1473 ostringstream msg;
1474 msg << "bad charge : " << sidechain_charge << " is out of valid range for residue K (" << n1 << "/" << n2 << ")." << ends;
1475 assertion_failed(__FILE__, __LINE__, msg.str().c_str());
1476 }
1477 }
1478
1479 ostringstream str;
1480 str << _("CheckProtonation() : setting residue ") << n2 << " K " << (charged ? _("charged.") : _("neutral.")) << endl << ends;
1481 PrintToLog(str.str().c_str());
1482
1483 if (charged)
1484 {
1485 atmr_24->formal_charge = +1;
1486 }
1487 else
1488 {
1489 atmr_24->formal_charge = +0;
1490 }
1491 }
1492
1493 /* if (sequence1[n2] == 'M')
1494 {
1495 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1496 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1497 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1498 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1499 } */
1500
1501 /* if (sequence1[n2] == 'F')
1502 {
1503 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1504 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1505 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1506 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1507 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1508 atom * atmr_25 = cp_FindAtom(resR, 0x25);
1509 atom * atmr_26 = cp_FindAtom(resR, 0x26);
1510 } */
1511
1512 /* if (sequence1[n2] == 'P')
1513 {
1514 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1515 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1516 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1517 } */
1518
1519 /* if (sequence1[n2] == 'S')
1520 {
1521 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1522 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1523 } */
1524
1525 /* if (sequence1[n2] == 'T')
1526 {
1527 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1528 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1529 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1530 } */
1531
1532 /* if (sequence1[n2] == 'W')
1533 {
1534 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1535 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1536 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1537 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1538 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1539 atom * atmr_25 = cp_FindAtom(resR, 0x25);
1540 atom * atmr_26 = cp_FindAtom(resR, 0x26);
1541 atom * atmr_27 = cp_FindAtom(resR, 0x27);
1542 atom * atmr_28 = cp_FindAtom(resR, 0x28);
1543 atom * atmr_29 = cp_FindAtom(resR, 0x29);
1544 } */
1545
1546 /* if (sequence1[n2] == 'Y')
1547 {
1548 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1549 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1550 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1551 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1552 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1553 atom * atmr_25 = cp_FindAtom(resR, 0x25);
1554 atom * atmr_26 = cp_FindAtom(resR, 0x26);
1555 atom * atmr_27 = cp_FindAtom(resR, 0x27);
1556 } */
1557
1558 /* if (sequence1[n2] == 'V')
1559 {
1560 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1561 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1562 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1563 } */
1564 }
1565 }
1566
1567 // nucleic...
1568 // nucleic...
1569 // nucleic...
1570
1571 if (ci_vector[n1].GetType() == chn_info::nucleic_acid)
1572 {
1573 // const char * sequence1 = ci_vector[n1].GetSequence1();
1574
1575 for (i32s n2 = 0;n2 < ci_vector[n1].GetLength();n2++)
1576 {
1577 iter_al resR[2]; GetRange(2, chnR, n2, resR);
1578
1579 /* atom * atmr_01 = cp_FindAtom(resR, 0x01);
1580 atom * atmr_02 = cp_FindAtom(resR, 0x02);
1581 atom * atmr_03 = cp_FindAtom(resR, 0x03);
1582 atom * atmr_04 = cp_FindAtom(resR, 0x04);
1583 atom * atmr_05 = cp_FindAtom(resR, 0x05);
1584 atom * atmr_06 = cp_FindAtom(resR, 0x06);
1585 atom * atmr_07 = cp_FindAtom(resR, 0x07);
1586 atom * atmr_08 = cp_FindAtom(resR, 0x08);
1587
1588 if (!atmr_01 || !atmr_02 || !atmr_03 || !atmr_04) continue;
1589 if (!atmr_05 || !atmr_06 || !atmr_07 || !atmr_08) continue; */
1590
1591 //fixme atom * atmr_10 = cp_FindAtom(resR, 0x10);
1592 //fixme atom * atmr_11 = cp_FindAtom(resR, 0x11);
1593 //fixme atom * atmr_12 = cp_FindAtom(resR, 0x12); // 5'-terminal...
1594 //fixme if (atmr_10 != NULL) { av.push_back(atmr_10); vv.push_back(1.0); }
1595 //fixme if (atmr_11 != NULL) { av.push_back(atmr_11); vv.push_back(1.0); }
1596 //fixme if (atmr_12 != NULL) { av.push_back(atmr_12); vv.push_back(1.0); }
1597
1598 /* if (n2 == 0) // 5'-terminal residue???
1599 {
1600 }
1601
1602 if (n2 == ci_vector[n1].GetLength() - 1) // 3'-terminal residue???
1603 {
1604 } */
1605
1606 // start handling the side chains!!!
1607 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1608
1609 /* if (sequence1[n2] == 'A')
1610 {
1611 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1612 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1613 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1614 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1615 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1616 atom * atmr_25 = cp_FindAtom(resR, 0x25);
1617 atom * atmr_26 = cp_FindAtom(resR, 0x26);
1618 atom * atmr_27 = cp_FindAtom(resR, 0x27);
1619 atom * atmr_28 = cp_FindAtom(resR, 0x28);
1620 atom * atmr_40 = cp_FindAtom(resR, 0x40);
1621 } */
1622
1623 /* if (sequence1[n2] == 'C')
1624 {
1625 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1626 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1627 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1628 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1629 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1630 atom * atmr_25 = cp_FindAtom(resR, 0x25);
1631 atom * atmr_40 = cp_FindAtom(resR, 0x40);
1632 atom * atmr_41 = cp_FindAtom(resR, 0x41);
1633 } */
1634
1635 /* if (sequence1[n2] == 'G')
1636 {
1637 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1638 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1639 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1640 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1641 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1642 atom * atmr_25 = cp_FindAtom(resR, 0x25);
1643 atom * atmr_26 = cp_FindAtom(resR, 0x26);
1644 atom * atmr_27 = cp_FindAtom(resR, 0x27);
1645 atom * atmr_28 = cp_FindAtom(resR, 0x28);
1646 atom * atmr_40 = cp_FindAtom(resR, 0x40);
1647 atom * atmr_41 = cp_FindAtom(resR, 0x41);
1648 } */
1649
1650 /* if (sequence1[n2] == 'T')
1651 {
1652 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1653 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1654 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1655 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1656 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1657 atom * atmr_25 = cp_FindAtom(resR, 0x25);
1658 atom * atmr_40 = cp_FindAtom(resR, 0x40);
1659 atom * atmr_41 = cp_FindAtom(resR, 0x41);
1660 atom * atmr_42 = cp_FindAtom(resR, 0x42);
1661 } */
1662
1663 /* if (sequence1[n2] == 'a')
1664 {
1665 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1666 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1667 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1668 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1669 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1670 atom * atmr_25 = cp_FindAtom(resR, 0x25);
1671 atom * atmr_26 = cp_FindAtom(resR, 0x26);
1672 atom * atmr_27 = cp_FindAtom(resR, 0x27);
1673 atom * atmr_28 = cp_FindAtom(resR, 0x28);
1674 atom * atmr_40 = cp_FindAtom(resR, 0x40);
1675 } */
1676
1677 /* if (sequence1[n2] == 'c')
1678 {
1679 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1680 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1681 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1682 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1683 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1684 atom * atmr_25 = cp_FindAtom(resR, 0x25);
1685 atom * atmr_40 = cp_FindAtom(resR, 0x40);
1686 atom * atmr_41 = cp_FindAtom(resR, 0x41);
1687 } */
1688
1689 /* if (sequence1[n2] == 'g')
1690 {
1691 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1692 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1693 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1694 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1695 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1696 atom * atmr_25 = cp_FindAtom(resR, 0x25);
1697 atom * atmr_26 = cp_FindAtom(resR, 0x26);
1698 atom * atmr_27 = cp_FindAtom(resR, 0x27);
1699 atom * atmr_28 = cp_FindAtom(resR, 0x28);
1700 atom * atmr_40 = cp_FindAtom(resR, 0x40);
1701 atom * atmr_41 = cp_FindAtom(resR, 0x41);
1702 } */
1703
1704 /* if (sequence1[n2] == 'u')
1705 {
1706 atom * atmr_20 = cp_FindAtom(resR, 0x20);
1707 atom * atmr_21 = cp_FindAtom(resR, 0x21);
1708 atom * atmr_22 = cp_FindAtom(resR, 0x22);
1709 atom * atmr_23 = cp_FindAtom(resR, 0x23);
1710 atom * atmr_24 = cp_FindAtom(resR, 0x24);
1711 atom * atmr_25 = cp_FindAtom(resR, 0x25);
1712 atom * atmr_40 = cp_FindAtom(resR, 0x40);
1713 atom * atmr_41 = cp_FindAtom(resR, 0x41);
1714 } */
1715 }
1716 }
1717 }
1718 }
1719
AddHydrogens(void)1720 void model::AddHydrogens(void)
1721 {
1722 srand(time(NULL));
1723 if (ref_civ != NULL)
1724 {
1725 ostringstream str;
1726 str << _("Sequence information found; calling CheckProtonation().") << endl;
1727 str << _("WARNING ; formal_charge may be changed for some atoms.") << endl << ends;
1728 PrintToLog(str.str().c_str());
1729
1730 CheckProtonation();
1731 }
1732 else
1733 {
1734 ostringstream str;
1735 str << _("Using default rules in AddHydrogens().") << endl << ends;
1736 PrintToLog(str.str().c_str());
1737 }
1738
1739 // first check to see if the user selected anything.
1740 // if so, only add hydrogens to the selected atoms.
1741 iter_al it1 = atom_list.begin();
1742 bool some_atoms_selected = false; // is anything selected by the user?
1743 while (it1 != atom_list.end())
1744 {
1745 if ((* it1++).flags & ATOMFLAG_USER_SELECTED)
1746 {
1747 some_atoms_selected = true;
1748 break;
1749 }
1750 }
1751
1752 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
1753 {
1754 // if the user selected atoms, only do this if the atom is selected.
1755 if (some_atoms_selected && (!(* it1).flags & ATOMFLAG_USER_SELECTED))
1756 {
1757 continue;
1758 }
1759
1760 AddHydrogens(& (* it1));
1761 }
1762 }
1763
AddHydrogens(atom * atmr)1764 void model::AddHydrogens(atom * atmr)
1765 {
1766 bool skip = true;
1767
1768 fGL valence = NOT_DEFINED;
1769 int an = atmr->el.GetAtomicNumber();
1770 switch (an)
1771 {
1772 case 6: valence = 4; skip = false; break;
1773 case 7: valence = 3; skip = false; break;
1774 case 8: valence = 2; skip = false; break;
1775 case 16: valence = 2; skip = false; break;
1776 }
1777
1778 if (skip) return;
1779
1780 valence += atmr->formal_charge;
1781
1782 // codes:
1783 // ^^^^^^
1784 // 0 = tetrahedral
1785 // 1 = planar
1786 // 2 = linear
1787
1788 i32s geometry = 0; vector<atom *> other_atoms;
1789 for (iter_cl itX = atmr->cr_list.begin();itX != atmr->cr_list.end();itX++)
1790 {
1791 other_atoms.push_back((* itX).atmr);
1792
1793 switch ((* itX).bndr->bt.GetValue())
1794 {
1795 case 1: valence -= 1.0; break;
1796 case 2: valence -= 2.0; if (!geometry) geometry++; break;
1797 case 3: valence -= 3.0; geometry = 2; break;
1798 case 0: valence -= 1.5; if (!geometry) geometry++; break;
1799 }
1800 }
1801
1802 // a special check that tries to detect planar nitrogen atoms...
1803
1804 if (atmr->el.GetAtomicNumber() == 7 && geometry == 0)
1805 {
1806 bool has_sp2_neighbor = false;
1807 for (i32u n1 = 0;n1 < other_atoms.size();n1++)
1808 {
1809 i32s num_S_bonds = 0;
1810 i32s num_D_bonds = 0;
1811 i32s num_T_bonds = 0;
1812
1813 for (iter_cl itX = other_atoms[n1]->cr_list.begin();itX != other_atoms[n1]->cr_list.end();itX++)
1814 {
1815 switch ((* itX).bndr->bt.GetValue())
1816 {
1817 case 1: num_S_bonds++; break;
1818 case 2: num_D_bonds++; break;
1819 case 3: num_T_bonds++; break;
1820 }
1821 }
1822
1823 const int oa_an = other_atoms[n1]->el.GetAtomicNumber();
1824
1825 if (oa_an == 6 && num_D_bonds) has_sp2_neighbor = true;
1826 if (oa_an == 7 && num_D_bonds) has_sp2_neighbor = true;
1827 }
1828
1829 if (has_sp2_neighbor) geometry = 1;
1830 }
1831
1832 // add the atoms/bonds...
1833 // ^^^^^^^^^^^^^^^^^^^^^^
1834
1835 fGL _len = 0.11;
1836
1837 fGL _ang;
1838 if (geometry == 2) _ang = 180.0 * M_PI / 180.0;
1839 else if (geometry == 1) _ang = 120.0 * M_PI / 180.0;
1840 else _ang = 109.5 * M_PI / 180.0;
1841
1842 const fGL * crdx = atmr->GetCRD(0);
1843 while (valence > 0.0)
1844 {
1845 fGL rx = ((fGL) rand() / (fGL) RAND_MAX) - 0.5;
1846 fGL ry = ((fGL) rand() / (fGL) RAND_MAX) - 0.5;
1847 fGL rz = ((fGL) rand() / (fGL) RAND_MAX) - 0.5;
1848 v3d<fGL> v1(rx, ry, rz); v1 = v1 / v1.len();
1849
1850 if (other_atoms.size() == 1)
1851 {
1852 // at this stage we have a good chance of improving
1853 // the resulting geometry if we can find a better v1.
1854
1855 atom * other = other_atoms.front();
1856
1857 vector<atom *> o2h_v;
1858 for (iter_cl itX = other->cr_list.begin();itX != other->cr_list.end();itX++)
1859 {
1860 if ((* itX).atmr == atmr) continue;
1861 if ((* itX).atmr->el.GetAtomicNumber() == 1) continue;
1862
1863 o2h_v.push_back((* itX).atmr);
1864 }
1865
1866 if (o2h_v.size() == 1)
1867 {
1868 const fGL * crd0 = other->GetCRD(0);
1869 const fGL * crdA = o2h_v.front()->GetCRD(0);
1870 v3d<fGL> v2(crd0, crdx); v2 = v2 / v2.len();
1871 v3d<fGL> v3(crd0, crdA); v3 = v3 / v3.len();
1872 v1 = v2.vpr(v3); v1 = v1 / v1.len();
1873 }
1874
1875 if (o2h_v.size() == 2)
1876 {
1877 const fGL * crd0 = other->GetCRD(0);
1878 const fGL * crdA = o2h_v[0]->GetCRD(0);
1879 const fGL * crdB = o2h_v[1]->GetCRD(0);
1880 v3d<fGL> v2(crd0, crdA); v2 = v2 / v2.len();
1881 v3d<fGL> v3(crd0, crdB); v3 = v3 / v3.len();
1882 v1 = v2.vpr(v3); v1 = v1 / v1.len();
1883 }
1884
1885 // continue...
1886
1887 const fGL * crd0 = other_atoms[0]->GetCRD(0);
1888 v3d<fGL> v2(crdx, crd0); v2 = v2 / v2.len();
1889 v3d<fGL> v3 = v1.vpr(v2); v3 = v3 / v3.len();
1890 v1 = (v2 * cos(_ang)) + (v3 * sin(_ang));
1891 }
1892 else if (other_atoms.size() == 2)
1893 {
1894 const fGL * crd0 = other_atoms[0]->GetCRD(0);
1895 const fGL * crd1 = other_atoms[1]->GetCRD(0);
1896 v3d<fGL> v2(crdx, crd0); v2 = v2 / v2.len();
1897 v3d<fGL> v3(crdx, crd1); v3 = v3 / v3.len();
1898 if (!geometry)
1899 {
1900 v3d<fGL> v4 = v2 + v3; v4 = v4 / -v4.len();
1901 v3d<fGL> v5 = v2.vpr(v3); v5 = v5 / v5.len();
1902 v1 = (v4 * 0.5) + (v5 * 0.5); v1 = v1 / v1.len();
1903 }
1904 else
1905 {
1906 v1 = v2 + v3; v1 = v1 / -v1.len();
1907 }
1908 }
1909 else if (other_atoms.size() > 2)
1910 {
1911 const fGL * crd0 = other_atoms[0]->GetCRD(0);
1912 const fGL * crd1 = other_atoms[1]->GetCRD(0);
1913 const fGL * crd2 = other_atoms[2]->GetCRD(0);
1914 v3d<fGL> v2(crdx, crd0); v2 = v2 / v2.len();
1915 v3d<fGL> v3(crdx, crd1); v3 = v3 / v3.len();
1916 v3d<fGL> v4(crdx, crd2); v4 = v4 / v4.len();
1917 v1 = v2 + v3 + v4; v1 = v1 / -v1.len();
1918 }
1919
1920 fGL crdh[3];
1921 crdh[0] = crdx[0] + _len * v1.data[0];
1922 crdh[1] = crdx[1] + _len * v1.data[1];
1923 crdh[2] = crdx[2] + _len * v1.data[2];
1924
1925 atom newatom(element(1), crdh, cs_vector.size());
1926 AddAtom_lg(newatom); other_atoms.push_back(& atom_list.back());
1927
1928 bond newbond(atmr, & atom_list.back(), bondtype('S'));
1929 AddBond(newbond);
1930
1931 valence -= 1.0;
1932 }
1933 }
1934
RemoveHydrogens(void)1935 void model::RemoveHydrogens(void)
1936 {
1937 iter_bl itb1 = bond_list.begin();
1938 while (itb1 != bond_list.end())
1939 {
1940 bool flag = false;
1941 if ((* itb1).atmr[0]->el.GetAtomicNumber() == 1) flag = true;
1942 if ((* itb1).atmr[1]->el.GetAtomicNumber() == 1) flag = true;
1943
1944 if (flag)
1945 {
1946 RemoveBond(itb1); // now this iterator is invalidated?!?!?!
1947 itb1 = bond_list.begin();
1948 }
1949 else itb1++;
1950 }
1951
1952 iter_al ita1 = atom_list.begin();
1953 while (ita1 != atom_list.end())
1954 {
1955 bool flag = false;
1956 if ((* ita1).el.GetAtomicNumber() == 1) flag = true;
1957
1958 if (flag)
1959 {
1960 RemoveAtom(ita1); // now this iterator is invalidated?!?!?!
1961 ita1 = atom_list.begin();
1962 }
1963 else ita1++;
1964 }
1965 }
1966
SolvateBox(fGL dimx,fGL dimy,fGL dimz,fGL density,model * solvent,const char * export_gromacs)1967 void model::SolvateBox(fGL dimx, fGL dimy, fGL dimz, fGL density, model * solvent, const char * export_gromacs)
1968 {
1969 use_periodic_boundary_conditions = true;
1970 saved_periodic_box_HALFdim[0] = dimx;
1971 saved_periodic_box_HALFdim[1] = dimy;
1972 saved_periodic_box_HALFdim[2] = dimz;
1973
1974 SystemWasModified();
1975
1976 if (density <= 0.0) return; // zero density -> infinite distance.
1977 const fGL distance = S_Initialize(density, & solvent);
1978
1979 bool system_was_empty = (!atom_list.size() && !bond_list.size());
1980
1981 srand(time(NULL));
1982
1983 // make the box...
1984 // ^^^^^^^^^^^^^^^
1985 // 20060314 ; do not use a lattice that is oriented in a same way as the
1986 // box is -> this will lead into inaccurate densities since many molecules
1987 // are either included or excluded if whole rows are included/excluded...
1988
1989 v3d<fGL> xv; v3d<fGL> yv; v3d<fGL> zv;
1990
1991 zv = v3d<fGL>((fGL) rand() / (fGL) RAND_MAX - 0.5, (fGL) rand() / (fGL) RAND_MAX - 0.5, (fGL) rand() / (fGL) RAND_MAX - 0.5);
1992 xv = v3d<fGL>((fGL) rand() / (fGL) RAND_MAX - 0.5, (fGL) rand() / (fGL) RAND_MAX - 0.5, (fGL) rand() / (fGL) RAND_MAX - 0.5);
1993 xv = xv / xv.len();
1994
1995 yv = zv.vpr(xv);
1996 yv = yv / yv.len();
1997
1998 zv = xv.vpr(yv);
1999 zv = zv / zv.len();
2000
2001 const fGL maxdist = sqrt(dimx * dimx + dimy * dimy + dimz * dimz);
2002
2003 const i32s lim = (i32s) floor(maxdist / distance) + 2; // +1 is the minimum.
2004
2005 i32s solvent_molecules_added = 0;
2006 for (i32s lx = (-lim + 1);lx < lim;lx++)
2007 {
2008 for (i32s ly = (-lim + 1);ly < lim;ly++)
2009 {
2010 for (i32s lz = (-lim + 1);lz < lim;lz++)
2011 {
2012 fGL crdS1[3];
2013 crdS1[0] = distance * lx;
2014 crdS1[1] = distance * ly;
2015 crdS1[2] = distance * lz;
2016
2017 if (lz % 2) // twist the cubic lattice to an octahedral one ; OPTIONAL!
2018 {
2019 crdS1[0] += 0.5 * distance;
2020 crdS1[1] += 0.5 * distance;
2021 }
2022
2023 v3d<fGL> pv = (xv * crdS1[0]) + (yv * crdS1[1]) + (zv * crdS1[2]);
2024 const fGL crdS2[3] = { pv[0], pv[1], pv[2] };
2025
2026 if (crdS2[0] < -dimx || crdS2[0] > +dimx) continue; // skip if the molecule is too far...
2027 if (crdS2[1] < -dimy || crdS2[1] > +dimy) continue; // skip if the molecule is too far...
2028 if (crdS2[2] < -dimz || crdS2[2] > +dimz) continue; // skip if the molecule is too far...
2029
2030 bool skip = false; // skip if there is overlap with solute atoms ; FOR SMALL SOLVENTS ONLY?!?
2031 for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
2032 {
2033 if ((* it1).flags & ATOMFLAG_IS_SOLVENT_ATOM) continue;
2034
2035 const fGL * crdX = (* it1).GetCRD(0);
2036
2037 fGL dx = crdS2[0] - crdX[0];
2038 fGL dy = crdS2[1] - crdX[1];
2039 fGL dz = crdS2[2] - crdX[2];
2040
2041 fGL d2s = sqrt(dx * dx + dy * dy + dz * dz);
2042 if (d2s < 0.175) skip = true;
2043
2044 if (skip) break;
2045 } if (skip) continue;
2046
2047 solvent_molecules_added++;
2048
2049 f64 rot[3];
2050 rot[0] = 2.0 * M_PI * (f64) rand() / (f64) RAND_MAX;
2051 rot[1] = 2.0 * M_PI * (f64) rand() / (f64) RAND_MAX;
2052 rot[2] = 2.0 * M_PI * (f64) rand() / (f64) RAND_MAX;
2053
2054 vector<atom *> av1;
2055 vector<atom *> av2;
2056
2057 for (iter_al it1 = solvent->GetAtomsBegin();it1 != solvent->GetAtomsEnd();it1++)
2058 {
2059 const fGL * orig = (* it1).GetCRD(0);
2060
2061 fGL d2b[3]; // rotate x (y,z)...
2062 d2b[0] = orig[0];
2063 d2b[1] = orig[1] * cos(rot[0]) - orig[2] * sin(rot[0]);
2064 d2b[2] = orig[1] * sin(rot[0]) + orig[2] * cos(rot[0]);
2065
2066 fGL d2c[3]; // rotate y (z,x)...
2067 d2c[0] = d2b[2] * sin(rot[1]) + d2b[0] * cos(rot[1]);
2068 d2c[1] = d2b[1];
2069 d2c[2] = d2b[2] * cos(rot[1]) - d2b[0] * sin(rot[1]);
2070
2071 fGL d2d[3]; // rotate z (x,y)...
2072 d2d[0] = d2c[0] * cos(rot[2]) - d2c[1] * sin(rot[2]);
2073 d2d[1] = d2c[0] * sin(rot[2]) + d2c[1] * cos(rot[2]);
2074 d2d[2] = d2c[2];
2075
2076 d2d[0] += crdS2[0];
2077 d2d[1] += crdS2[1];
2078 d2d[2] += crdS2[2];
2079
2080 atom newA((* it1).el, d2d, GetCRDSetCount());
2081 newA.flags |= ATOMFLAG_IS_SOLVENT_ATOM;
2082 // newA.flags |= ATOMFLAG_MEASURE_ND_RDF; // what about this???
2083
2084 AddAtom_lg(newA);
2085 av1.push_back(& (* it1));
2086 av2.push_back(& GetLastAtom());
2087 }
2088
2089 for (iter_bl it1 = solvent->GetBondsBegin();it1 != solvent->GetBondsEnd();it1++)
2090 {
2091 i32u ind1 = 0;
2092 while (ind1 < av1.size())
2093 {
2094 if ((* it1).atmr[0] == av1[ind1]) break;
2095 else ind1++;
2096 }
2097
2098 i32u ind2 = 0;
2099 while (ind2 < av1.size())
2100 {
2101 if ((* it1).atmr[1] == av1[ind2]) break;
2102 else ind2++;
2103 }
2104
2105 if (ind1 == av1.size() || ind2 == av1.size())
2106 {
2107 assertion_failed(__FILE__, __LINE__, "index search failed!");
2108 }
2109
2110 bond newB(av2[ind1], av2[ind2], (* it1).bt);
2111 AddBond(newB);
2112 }
2113 }
2114 }
2115 }
2116
2117 cout << _("added ") << solvent_molecules_added << _(" solvent molecules.") << endl;
2118
2119 ///////////////////////////////////////////////////////////////////////////
2120
2121 // 20060811 ; check the dimensions so that the density would be as
2122 // accurate as possible (since we can't add a fraction of a molecule).
2123
2124 double sum_of_masses = 0.0; // kg/mol ; all atoms.
2125 for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
2126 {
2127 // in the above code there are no skip-checks either...
2128 sum_of_masses += (* it1).mass * 1.6605402e-27 * 6.0221367e+23;
2129 }
2130
2131 double currV;
2132 double currD;
2133
2134 int cycle = 0;
2135 while (true)
2136 {
2137 currV = dimx * dimy * dimz * 8.0 * 6.0221367e-4; // m^3 / mol
2138 currD = 0.001 * sum_of_masses / currV; // kg/dm^3
2139
2140 double ratio = pow(currD / density, 0.33333);
2141
2142 //cout << "currD = " << currD << " ; ";
2143 //cout << "ratio = " << ratio << endl;
2144
2145 if (cycle == 1) break;
2146
2147 dimx *= ratio;
2148 dimy *= ratio;
2149 dimz *= ratio;
2150 cycle++;
2151 }
2152
2153 ostringstream str;
2154 str << _("Density is ") << currD << " kg/dm^3." << endl;
2155 str << _("Adjusted dimensions are: ") << dimx << " nm, ";
2156 str << dimy << " nm, " << dimz << " nm." << endl << ends;
2157 PrintToLog(str.str().c_str());
2158
2159 use_periodic_boundary_conditions = true;
2160 saved_periodic_box_HALFdim[0] = dimx;
2161 saved_periodic_box_HALFdim[1] = dimy;
2162 saved_periodic_box_HALFdim[2] = dimz;
2163
2164 SystemWasModified();
2165
2166 ///////////////////////////////////////////////////////////////////////////
2167
2168 if (export_gromacs != NULL)
2169 {
2170 if (!system_was_empty)
2171 {
2172 Message(_("Sorry, the export option is available for pure solvents only!"));
2173 delete solvent; return;
2174 }
2175
2176 if (dynamic_cast<setup1_mm *>(GetCurrentSetup()) == NULL)
2177 {
2178 Message(_("Sorry, the export option is available for MM setups only!"));
2179 delete solvent; return;
2180 }
2181
2182 solvent->ReplaceCurrentSetup(new setup1_mm(solvent));
2183
2184 i32s my_index = GetCurrentSetup()->GetCurrEngIndex();
2185 solvent->current_setup->SetCurrEngIndex(my_index);
2186
2187 solvent->current_setup->CreateCurrentEngine();
2188 engine * eng1 = solvent->current_setup->GetCurrentEngine();
2189 eng1_mm_default_bt * eng2 = dynamic_cast<eng1_mm_default_bt *>(eng1);
2190 if (eng1 == NULL || eng2== NULL)
2191 {
2192 Message(_("Export failed!"));
2193 delete solvent; return;
2194 }
2195
2196 CopyCRD(solvent, eng2, 0);
2197 eng2->Compute(0);
2198
2199 atom ** atmtab = solvent->current_setup->GetMMAtoms();
2200
2201 // write the top file
2202 // ^^^^^^^^^^^^^^^^^^
2203
2204 ostringstream tfns;
2205 tfns << export_gromacs << ".top" << ends;
2206
2207 ofstream top_f;
2208 top_f.open(tfns.str().c_str(), ios::out);
2209
2210 top_f << "; this is a top file exported by libghemical " << LIBVERSION << ", not an original gromacs file!" << endl;
2211 top_f << endl;
2212
2213 top_f << "[ defaults ]" << endl;
2214 top_f << "; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ" << endl;
2215 top_f << " 1 2 yes " << eng1_mm::fudgeLJ << "\t" << eng1_mm::fudgeQQ << endl;
2216 top_f << endl;
2217
2218 top_f << "[ atomtypes ]" << endl;
2219 top_f << "; name mass charge ptype sigma epsilon" << endl;
2220
2221 for (i32s ai = 0;ai < solvent->current_setup->GetMMAtomCount();ai++)
2222 {
2223 const default_at * at; f64 r1 = 0.150; f64 e1 = 0.175;
2224 at = default_tables::GetInstance()->GetAtomType(atmtab[ai]->atmtp);
2225 if (at != NULL) { r1 = at->vdw_R; e1 = at->vdw_E; }
2226
2227 top_f << "\t";
2228 top_f << atmtab[ai]->el.GetSymbol() << (ai + 1) << "\t";
2229 top_f << atmtab[ai]->el.GetAtomicMass() << "\t";
2230 top_f << atmtab[ai]->charge << "\t";
2231 top_f << "A" << "\t";
2232 top_f << (2.0 * r1) << "\t";
2233 top_f << e1 << endl;
2234 }
2235
2236 top_f << endl;
2237
2238 top_f << "[ moleculetype ]" << endl;
2239 top_f << "; name nrexcl" << endl;
2240 top_f << " SOL 3" << endl;
2241 top_f << endl;
2242
2243 top_f << "[ atoms ]" << endl;
2244 top_f << "; nr type resnr residu atom cgnr charge" << endl;
2245
2246 for (i32s ind = 0;ind < solvent->current_setup->GetMMAtomCount();ind++)
2247 {
2248 const default_at * at; f64 r1 = 0.150; f64 e1 = 0.175;
2249 at = default_tables::GetInstance()->GetAtomType(atmtab[ind]->atmtp);
2250 if (at != NULL) { r1 = at->vdw_R; e1 = at->vdw_E; }
2251
2252 top_f << "\t";
2253 top_f << (ind + 1) << "\t";
2254 top_f << atmtab[ind]->el.GetSymbol() << (ind + 1) << "\t"; // what this really should be?
2255 top_f << "1\t";
2256 top_f << "SOL\t";
2257 top_f << atmtab[ind]->el.GetSymbol() << (ind + 1) << "\t"; // what this really should be?
2258 top_f << "1\t";
2259 top_f << atmtab[ind]->charge << endl;
2260 }
2261
2262 top_f << endl;
2263
2264 if (eng2->bt1_vector.size() > 0)
2265 {
2266 top_f << "[ bonds ]" << endl;
2267 top_f << "; ai aj funct c0 c1" << endl;
2268
2269 for (i32u ind = 0;ind < eng2->bt1_vector.size();ind++)
2270 {
2271 top_f << "\t";
2272 top_f << (eng2->bt1_vector[ind].atmi[0] + 1) << "\t";
2273 top_f << (eng2->bt1_vector[ind].atmi[1] + 1) << "\t";
2274 top_f << "1\t";
2275 top_f << eng2->bt1_vector[ind].opt << "\t";
2276 top_f << eng2->bt1_vector[ind].fc << endl;
2277 }
2278
2279 top_f << endl;
2280 }
2281
2282 if (eng2->bt2_vector.size() > 0)
2283 {
2284 top_f << "[ angles ]" << endl;
2285 top_f << "; ai aj ak funct c0 c1" << endl;
2286
2287 for (i32u ind = 0;ind < eng2->bt2_vector.size();ind++)
2288 {
2289 top_f << "\t";
2290 top_f << (eng2->bt2_vector[ind].atmi[0] + 1) << "\t";
2291 top_f << (eng2->bt2_vector[ind].atmi[1] + 1) << "\t";
2292 top_f << (eng2->bt2_vector[ind].atmi[2] + 1) << "\t";
2293 top_f << "1\t";
2294 top_f << (180.0 * eng2->bt2_vector[ind].opt / M_PI) << "\t"; // convert rad->deg!!!
2295 top_f << eng2->bt2_vector[ind].fc << endl;
2296 }
2297
2298 top_f << endl;
2299 }
2300
2301 if (eng2->bt3_vector.size() > 0)
2302 {
2303 top_f << "[ dihedrals ]" << endl;
2304 top_f << "; ai aj ak al funct phi cp mult" << endl;
2305
2306 for (i32u ind = 0;ind < eng2->bt3_vector.size();ind++)
2307 {
2308 // can set only a single term? then select the most important term!
2309 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2310
2311 int mult = 1; f64 maxv = eng2->bt3_vector[ind].fc1;
2312 if (fabs(eng2->bt3_vector[ind].fc2) > fabs(maxv)) { mult = 2; maxv = eng2->bt3_vector[ind].fc2; }
2313 if (fabs(eng2->bt3_vector[ind].fc3) > fabs(maxv)) { mult = 3; maxv = eng2->bt3_vector[ind].fc3; }
2314
2315 top_f << "\t";
2316 top_f << (eng2->bt3_vector[ind].atmi[0] + 1) << "\t";
2317 top_f << (eng2->bt3_vector[ind].atmi[1] + 1) << "\t";
2318 top_f << (eng2->bt3_vector[ind].atmi[2] + 1) << "\t";
2319 top_f << (eng2->bt3_vector[ind].atmi[3] + 1) << "\t";
2320 top_f << "1\t";
2321 top_f << (maxv < 0.0 ? 180.0 : 0.0) << "\t";
2322 top_f << fabs(maxv) << "\t";
2323 top_f << mult << endl;
2324 }
2325
2326 top_f << endl;
2327 }
2328
2329 if (eng2->bt4_vector.size() > 0)
2330 {
2331 top_f << "[ dihedrals ]" << endl;
2332 top_f << "; ai aj ak al funct q0 cq" << endl;
2333
2334 for (i32u ind = 0;ind < eng2->bt4_vector.size();ind++)
2335 {
2336 top_f << "\t";
2337 top_f << (eng2->bt4_vector[ind].atmi[1] + 1) << "\t"; // plane
2338 top_f << (eng2->bt4_vector[ind].atmi[0] + 1) << "\t"; // plane
2339 top_f << (eng2->bt4_vector[ind].atmi[2] + 1) << "\t"; // plane
2340 top_f << (eng2->bt4_vector[ind].atmi[3] + 1) << "\t"; // ref
2341 top_f << "2\t";
2342 top_f << (180.0 * eng2->bt4_vector[ind].opt / M_PI) << "\t"; // convert rad->deg!!!
2343 top_f << eng2->bt4_vector[ind].fc << endl;
2344 }
2345
2346 top_f << endl;
2347 }
2348
2349 top_f << "[ system ]" << endl;
2350 top_f << "exported_from_libghemical" << endl;
2351 top_f << endl;
2352
2353 top_f << "[ molecules ]" << endl;
2354 top_f << "SOL\t" << (atom_list.size() / solvent->current_setup->GetMMAtomCount()) << endl;
2355
2356 top_f.close();
2357
2358 // write the gro file
2359 // ^^^^^^^^^^^^^^^^^^
2360
2361 ostringstream gfns;
2362 gfns << export_gromacs << ".gro" << ends;
2363
2364 ofstream gro_f;
2365 gro_f.open(gfns.str().c_str(), ios::out);
2366
2367 gro_f << "; this is a gro file exported by libghemical " << LIBVERSION << ", not an original gromacs file!" << endl;
2368 gro_f << atom_list.size() << endl;
2369
2370 i32u counter = 0; i32u rescnt = 0;
2371 for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
2372 {
2373 i32u atmcnt = (counter % solvent->current_setup->GetMMAtomCount());
2374 if (!atmcnt) rescnt++;
2375
2376 ostringstream ans; // atomname
2377 ans << (* it1).el.GetSymbol() << (atmcnt + 1) << ends;
2378
2379 gro_f << setw(5) << rescnt;
2380 gro_f << "SOL ";
2381
2382 int nb = 5 - strlen(ans.str().c_str()); if (nb < 0) nb = 0;
2383 for (int blank = 0;blank < nb;blank++) gro_f << " ";
2384 gro_f << ans.str().c_str();
2385
2386 gro_f << setw(5) << (counter + 1);
2387
2388 const fGL * crd = (* it1).GetCRD(0);
2389 gro_f << setw(8) << setprecision(3) << crd[0];
2390 gro_f << setw(8) << setprecision(3) << crd[1];
2391 gro_f << setw(8) << setprecision(3) << crd[2];
2392 gro_f << " 0.000";
2393 gro_f << " 0.000";
2394 gro_f << " 0.000";
2395 gro_f << endl;
2396
2397 counter++;
2398 }
2399
2400 gro_f << (dimx * 2.0) << " " << (dimy * 2.0) << " " << (dimz * 2.0) << endl;
2401
2402 gro_f.close();
2403 }
2404
2405 delete solvent;
2406
2407 // for (iter_al itX = GetAtomsBegin();itX != GetAtomsEnd();itX++)
2408 // { cout << (* itX).index << " " << ((* itX).flags & ATOMFLAG_IS_SOLVENT_ATOM) << endl; }
2409 }
2410
SolvateSphere(fGL rad1,fGL rad2,fGL density,model * solvent)2411 void model::SolvateSphere(fGL rad1, fGL rad2, fGL density, model * solvent)
2412 {
2413 use_boundary_potential = true;
2414 saved_boundary_potential_rad_solute = rad1;
2415 saved_boundary_potential_rad_solvent = rad2;
2416
2417 SystemWasModified();
2418
2419 if (density <= 0.0) return; // zero density -> infinite distance.
2420 const fGL distance = S_Initialize(density, & solvent);
2421
2422 srand(time(NULL));
2423
2424 // make the sphere of size rad2...
2425 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2426
2427 const i32s lim = (i32s) floor(rad2 / distance) + 2; // +1 is the minimum.
2428
2429 i32s solvent_molecules_added = 0;
2430 for (i32s lx = (-lim + 1);lx < lim;lx++)
2431 {
2432 for (i32s ly = (-lim + 1);ly < lim;ly++)
2433 {
2434 for (i32s lz = (-lim + 1);lz < lim;lz++)
2435 {
2436 fGL crdS[3];
2437 crdS[0] = distance * lx;
2438 crdS[1] = distance * ly;
2439 crdS[2] = distance * lz;
2440
2441 if (lz % 2) // twist the cubic lattice to an octahedral one ; OPTIONAL!
2442 {
2443 crdS[0] += 0.5 * distance;
2444 crdS[1] += 0.5 * distance;
2445 }
2446
2447 fGL dist = sqrt(crdS[0] * crdS[0] + crdS[1] * crdS[1] + crdS[2] * crdS[2]);
2448 if (dist > rad2) continue; // skip if the molecule is too far...
2449
2450 bool skip = false; // skip if there is overlap with solute atoms ; FOR SMALL SOLVENTS ONLY?!?
2451 for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
2452 {
2453 if ((* it1).flags & ATOMFLAG_IS_SOLVENT_ATOM) continue;
2454
2455 const fGL * crdX = (* it1).GetCRD(0);
2456
2457 fGL dx = crdS[0] - crdX[0];
2458 fGL dy = crdS[1] - crdX[1];
2459 fGL dz = crdS[2] - crdX[2];
2460
2461 fGL d2s = sqrt(dx * dx + dy * dy + dz * dz);
2462 if (d2s < 0.175) skip = true;
2463
2464 if (skip) break;
2465 } if (skip) continue;
2466
2467 solvent_molecules_added++;
2468
2469 f64 rot[3];
2470 rot[0] = 2.0 * M_PI * (f64) rand() / (f64) RAND_MAX;
2471 rot[1] = 2.0 * M_PI * (f64) rand() / (f64) RAND_MAX;
2472 rot[2] = 2.0 * M_PI * (f64) rand() / (f64) RAND_MAX;
2473
2474 vector<atom *> av1;
2475 vector<atom *> av2;
2476
2477 for (iter_al it1 = solvent->GetAtomsBegin();it1 != solvent->GetAtomsEnd();it1++)
2478 {
2479 const fGL * orig = (* it1).GetCRD(0);
2480
2481 fGL d2b[3]; // rotate x (y,z)...
2482 d2b[0] = orig[0];
2483 d2b[1] = orig[1] * cos(rot[0]) - orig[2] * sin(rot[0]);
2484 d2b[2] = orig[1] * sin(rot[0]) + orig[2] * cos(rot[0]);
2485
2486 fGL d2c[3]; // rotate y (z,x)...
2487 d2c[0] = d2b[2] * sin(rot[1]) + d2b[0] * cos(rot[1]);
2488 d2c[1] = d2b[1];
2489 d2c[2] = d2b[2] * cos(rot[1]) - d2b[0] * sin(rot[1]);
2490
2491 fGL d2d[3]; // rotate z (x,y)...
2492 d2d[0] = d2c[0] * cos(rot[2]) - d2c[1] * sin(rot[2]);
2493 d2d[1] = d2c[0] * sin(rot[2]) + d2c[1] * cos(rot[2]);
2494 d2d[2] = d2c[2];
2495
2496 d2d[0] += crdS[0];
2497 d2d[1] += crdS[1];
2498 d2d[2] += crdS[2];
2499
2500 atom newA((* it1).el, d2d, GetCRDSetCount());
2501 newA.flags |= ATOMFLAG_IS_SOLVENT_ATOM;
2502 // newA.flags |= ATOMFLAG_MEASURE_ND_RDF; // what about this???
2503
2504 AddAtom_lg(newA);
2505 av1.push_back(& (* it1));
2506 av2.push_back(& GetLastAtom());
2507 }
2508
2509 for (iter_bl it1 = solvent->GetBondsBegin();it1 != solvent->GetBondsEnd();it1++)
2510 {
2511 i32u ind1 = 0;
2512 while (ind1 < av1.size())
2513 {
2514 if ((* it1).atmr[0] == av1[ind1]) break;
2515 else ind1++;
2516 }
2517
2518 i32u ind2 = 0;
2519 while (ind2 < av1.size())
2520 {
2521 if ((* it1).atmr[1] == av1[ind2]) break;
2522 else ind2++;
2523 }
2524
2525 if (ind1 == av1.size() || ind2 == av1.size())
2526 {
2527 assertion_failed(__FILE__, __LINE__, "index search failed!");
2528 }
2529
2530 bond newB(av2[ind1], av2[ind2], (* it1).bt);
2531 AddBond(newB);
2532 }
2533 }
2534 }
2535 }
2536
2537 cout << _("added ") << solvent_molecules_added << _(" solvent molecules.") << endl;
2538
2539 delete solvent;
2540
2541 // for (iter_al itX = GetAtomsBegin();itX != GetAtomsEnd();itX++)
2542 // { cout << (* itX).index << " " << ((* itX).flags & ATOMFLAG_IS_SOLVENT_ATOM) << endl; }
2543 }
2544
S_Initialize(fGL density,model ** ref2solv)2545 fGL model::S_Initialize(fGL density, model ** ref2solv)
2546 {
2547 // here we set H2O to default solvent if there is no other molecule set,
2548 // and calculate the distances between the solvent molecules using density.
2549
2550 if (!(* ref2solv))
2551 {
2552 // build a water molecule!!!
2553 // ^^^^^^^^^^^^^^^^^^^^^^^^^
2554
2555 (* ref2solv) = new model();
2556 const fGL angle = 109.5 * M_PI / 180.0;
2557
2558 fGL crdO[3] = { 0.0, 0.0, 0.0 };
2559 atom newO(element(8), crdO, (* ref2solv)->GetCRDSetCount());
2560 (* ref2solv)->AddAtom_lg(newO); atom * ref_O = & (* ref2solv)->GetLastAtom();
2561
2562 fGL crdH1[3] = { 0.095, 0.0, 0.0 };
2563 atom newH1(element(1), crdH1, (* ref2solv)->GetCRDSetCount());
2564 (* ref2solv)->AddAtom_lg(newH1); atom * ref_H1 = & (* ref2solv)->GetLastAtom();
2565
2566 fGL crdH2[3] = { cos(angle)*0.095, sin(angle)*0.095, 0.0 };
2567 atom newH2(element(1), crdH2, (* ref2solv)->GetCRDSetCount());
2568 (* ref2solv)->AddAtom_lg(newH2); atom * ref_H2 = & (* ref2solv)->GetLastAtom();
2569
2570 bond newb1(ref_O, ref_H1, bondtype('S')); (* ref2solv)->AddBond(newb1);
2571 bond newb2(ref_O, ref_H2, bondtype('S')); (* ref2solv)->AddBond(newb2);
2572 }
2573
2574 f64 molarmass = 0.0;
2575 for (iter_al it1 = (* ref2solv)->GetAtomsBegin();it1 != (* ref2solv)->GetAtomsEnd();it1++)
2576 {
2577 molarmass += (* it1).el.GetAtomicMass();
2578 }
2579
2580 if (molarmass < 0.1)
2581 {
2582 ostringstream str;
2583 str << _("Could not calculate molar mass!") << endl;
2584 str << _("Failed to read the solvent file.") << ends;
2585 Message(str.str().c_str());
2586
2587 return NOT_DEFINED;
2588 }
2589
2590 // density = mass / volume -> mass = density * volume
2591 // calculate the mass (in grams) for 1 dm^3 of solvent.
2592
2593 f64 m1 = 1000.0 * density * 1.0;
2594
2595 // quantity = mass / molar_mass ; convert the mass into quantity
2596 // and further into number of particles using Avogadro's constant.
2597
2598 f64 n1 = m1 / molarmass;
2599 f64 n2 = n1 * 6.022e+23;
2600
2601 // now assume that all these particles are placed into a cubic lattice
2602 // of volume 1 dm^3 with even spacings ; calculate the distance (in nm)
2603 // between the particles in this lattice.
2604
2605 f64 distance = pow(1.0e+24 / n2, 1.0 / 3.0);
2606 return distance;
2607 }
2608
DoEnergy(void)2609 void model::DoEnergy(void)
2610 {
2611 engine * eng = GetCurrentSetup()->GetCurrentEngine();
2612 if (eng == NULL) GetCurrentSetup()->CreateCurrentEngine();
2613 eng = GetCurrentSetup()->GetCurrentEngine();
2614 if (eng == NULL) return;
2615
2616 ostringstream str1;
2617 str1 << _("Calculating Energy ");
2618 str1 << _("(setup = ") << GetCurrentSetup()->GetClassName_lg();
2619 str1 << _(", engine = ") << GetCurrentSetup()->GetEngineName(GetCurrentSetup()->GetCurrEngIndex());
2620 str1 << ")." << endl << ends;
2621
2622 PrintToLog(str1.str().c_str());
2623
2624 CopyCRD(this, eng, 0);
2625
2626 eng->Compute(0); // normal
2627 // eng->Check(1); // debug
2628
2629 if (dynamic_cast<eng1_sf *>(eng) != NULL) CopyCRD(eng, this, 0); // for ribbons...
2630
2631 ostringstream str2;
2632 str2.setf(ios::fixed); str2.precision(8);
2633 str2 << _("Energy = ") << eng->energy << " kJ/mol" << endl << ends;
2634
2635 PrintToLog(str2.str().c_str());
2636
2637 // we will not delete current_eng here, so that we can draw plots using it...
2638 // we will not delete current_eng here, so that we can draw plots using it...
2639 // we will not delete current_eng here, so that we can draw plots using it...
2640
2641 SetupPlotting();
2642 }
2643
DoGeomOpt(geomopt_param & param,bool updt)2644 void model::DoGeomOpt(geomopt_param & param, bool updt)
2645 {
2646 // make this thread-safe since this can be called from project::process_job_WhatEver() at the app side...
2647 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2648 // this means: 1) make sure that PrintToLog() is safe, 2) only call UpdateAllGraphicsViews(false) because
2649 // OGL-context is owned by the GUI-thread ; calling with "false" will pass the update to the GUI-thread.
2650
2651 ThreadLock();
2652
2653 engine * eng = GetCurrentSetup()->GetCurrentEngine();
2654 if (eng == NULL) GetCurrentSetup()->CreateCurrentEngine();
2655 eng = GetCurrentSetup()->GetCurrentEngine();
2656 if (eng == NULL)
2657 {
2658 ThreadUnlock();
2659 return;
2660 }
2661
2662 // eng_bp != NULL if we will use a system with boundary potential...
2663 engine_bp * eng_bp = dynamic_cast<engine_bp *>(eng);
2664
2665 // eng_pbc != NULL if we will use a system with periodic boundary conditions...
2666 engine_pbc * eng_pbc = dynamic_cast<engine_pbc *>(eng);
2667
2668 ostringstream str1;
2669 str1 << _("Starting Geometry Optimization ");
2670 str1 << _("(setup = ") << GetCurrentSetup()->GetClassName_lg();
2671 str1 << _(", engine = ") << GetCurrentSetup()->GetEngineName(GetCurrentSetup()->GetCurrEngIndex());
2672 str1 << ")." << endl << ends;
2673
2674 PrintToLog(str1.str().c_str());
2675
2676 CopyCRD(this, eng, 0);
2677
2678 geomopt * opt = new geomopt(eng, 100, 0.025, 10.0); // optimal settings?!?!?
2679
2680 f64 last_energy = 0.0; // this is for output and delta_e test...
2681
2682 ostringstream str2;
2683 str2 << _("Cycle \tEnergy \tGradient \tStep \t\tDeltaE") << endl << ends;
2684
2685 PrintToLog(str2.str().c_str());
2686
2687 ThreadUnlock();
2688
2689 i32s n1 = 0; // n1 counts the number of steps...
2690 bool cancel = false;
2691 while (!cancel)
2692 {
2693 if (!(n1 % 10)) eng->RequestNeighborListUpdate();
2694 if (!(n1 % 10) && eng_pbc != NULL) eng_pbc->CheckLocations();
2695
2696 opt->TakeCGStep(conjugate_gradient::Newton2An);
2697
2698 // problem: the gradient information is in fact not precise in this stage. the current gradient
2699 // is the one that was last calculated in the search, and it is not necessarily the best one.
2700 // to update the gradient, we need to Compute(1) here. JUST SLOWS GEOMOPT DOWN -> DISABLE!
2701 ///////////////////////////////////////////////////////////////////////////////////////////////
2702 //eng->Compute(1); // this is not vital, but will update the gradient vector length...
2703 ///////////////////////////////////////////////////////////////////////////////////////////////
2704
2705 ThreadLock();
2706
2707 if (!(n1 % 5))
2708 {
2709 double progress = 0.0;
2710 if (param.enable_nsteps) progress = (double) n1 / (double) (param.treshold_nsteps);
2711
2712 double graphdata = opt->optval; // this is an array of size 1...
2713
2714 cancel = SetProgress(progress, & graphdata);
2715
2716 if (n1 != 0)
2717 {
2718 ostringstream strX;
2719 strX << n1 << "\t";
2720 strX << opt->optval << "\t";
2721 strX << eng->GetGradientVectorLength() << "\t";
2722 strX << opt->optstp << "\t";
2723 strX << (last_energy - opt->optval) << "\t";
2724 strX << endl << ends;
2725
2726 PrintToLog(strX.str().c_str());
2727
2728 /* char buffer[1024];
2729 sprintf(buffer, "%4d %12.5f kJ/mol %10.4e %10.4e %10.4e \n", n1,
2730 opt->optval, eng->GetGradientVectorLength(), opt->optstp, last_energy - opt->optval);
2731 PrintToLog(buffer); */
2732 }
2733 else
2734 {
2735 ostringstream strX;
2736 strX << n1 << "\t";
2737 strX << opt->optval << "\t";
2738 strX << eng->GetGradientVectorLength() << "\t";
2739 strX << opt->optstp << "\t";
2740 strX << "**********" << "\t";
2741 strX << endl << ends;
2742
2743 PrintToLog(strX.str().c_str());
2744
2745 /* char buffer[1024];
2746 sprintf(buffer, "%4d %12.5f kJ/mol %10.4e %10.4e ********** \n", n1,
2747 opt->optval, eng->GetGradientVectorLength(), opt->optstp);
2748 PrintToLog(buffer); */
2749 }
2750
2751 NoThreadsIterate(); // this will update the GUI when threads disabled...
2752 }
2753
2754 bool terminate = false;
2755
2756 if (param.enable_nsteps) // the nsteps test...
2757 {
2758 if (n1 >= param.treshold_nsteps)
2759 {
2760 terminate = true;
2761
2762 ostringstream strX;
2763 strX << _("The nsteps termination test was passed.") << endl << ends;
2764 PrintToLog(strX.str().c_str());
2765 }
2766 }
2767
2768 if (param.enable_grad) // the grad test...
2769 {
2770 if (eng->GetGradientVectorLength() < param.treshold_grad)
2771 {
2772 terminate = true;
2773
2774 ostringstream strX;
2775 strX << _("The grad termination test was passed.") << endl << ends;
2776 PrintToLog(strX.str().c_str());
2777 }
2778 }
2779
2780 if (param.enable_delta_e) // the delta_e test...
2781 {
2782 bool flag = false; const f64 treshold_step = 1.0e-12; // can we keep this as a constant???
2783 if (n1 != 0 && (last_energy - opt->optval) != 0.0 && fabs(last_energy - opt->optval) < param.treshold_delta_e) flag = true;
2784 if ((opt->optstp != 0.0) && (opt->optstp < treshold_step)) flag = true;
2785
2786 if (flag)
2787 {
2788 terminate = true;
2789
2790 ostringstream strX;
2791 strX << _("The deltaE termination test was passed.") << endl << ends;
2792 PrintToLog(strX.str().c_str());
2793 }
2794 }
2795
2796 last_energy = opt->optval;
2797
2798 if (!(n1 % 10) || terminate)
2799 {
2800 CopyCRD(eng, this, 0);
2801 CenterCRDSet(0, false);
2802
2803 UpdateAllGraphicsViews(updt);
2804 }
2805
2806 ThreadUnlock();
2807
2808 if (terminate) break; // exit the loop here!!!
2809
2810 n1++; // update the number of steps...
2811 }
2812
2813 delete opt;
2814
2815 // we will not delete current_eng here, so that we can draw plots using it...
2816 // we will not delete current_eng here, so that we can draw plots using it...
2817 // we will not delete current_eng here, so that we can draw plots using it...
2818
2819 // above, CopyCRD was done eng->mdl and then CenterCRDSet() was done for mdl.
2820 // this might cause that old coordinates remain in eng object, possibly affecting plots.
2821 // here we sync the coordinates and other plotting data in the eng object.
2822
2823 ThreadLock();
2824 CopyCRD(this, eng, 0);
2825 SetupPlotting();
2826 ThreadUnlock();
2827 }
2828
DoMolDyn(moldyn_param & param,bool updt)2829 void model::DoMolDyn(moldyn_param & param, bool updt)
2830 {
2831 // make this thread-safe since this can be called from project::process_job_WhatEver() at the app side...
2832 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2833 // this means: 1) make sure that PrintToLog() is safe, 2) only call UpdateAllGraphicsViews(false) because
2834 // OGL-context is owned by the GUI-thread ; calling with "false" will pass the update to the GUI-thread.
2835
2836 ThreadLock();
2837
2838 engine * eng = GetCurrentSetup()->GetCurrentEngine();
2839 if (eng == NULL) GetCurrentSetup()->CreateCurrentEngine();
2840 eng = GetCurrentSetup()->GetCurrentEngine();
2841 if (eng == NULL)
2842 {
2843 ThreadUnlock();
2844 return;
2845 }
2846
2847 // eng_bp != NULL if we will use a system with boundary potential...
2848 engine_bp * eng_bp = dynamic_cast<engine_bp *>(eng);
2849
2850 // eng_pbc != NULL if we will use a system with periodic boundary conditions...
2851 engine_pbc * eng_pbc = dynamic_cast<engine_pbc *>(eng);
2852
2853 // THIS IS OPTIONAL!!! FOR BOUNDARY POTENTIAL STUFF ONLY!!!
2854 //if (eng_bp != NULL) eng_bp->nd_eval = new number_density_evaluator(eng_bp, false, 20);
2855
2856 // THIS IS OPTIONAL!!! FOR RADIAL DENSITY FUNCTION STUFF ONLY!!!
2857 //if (eng_bp != NULL) eng_bp->rdf_eval = new radial_density_function_evaluator(eng_bp, 80, 0.15, 0.95); // 0.95 - 0.15 = 0.80
2858 //if (eng_bp != NULL) eng_bp->rdf_eval = new radial_density_function_evaluator(eng_bp, 50, 0.15, 0.65, 0.00, 0.75); // 0.65 - 0.15 = 0.50
2859
2860 ostringstream str1;
2861 str1 << _("Starting Molecular Dynamics ");
2862 str1 << _("(setup = ") << GetCurrentSetup()->GetClassName_lg();
2863 str1 << _(", engine = ") << GetCurrentSetup()->GetEngineName(GetCurrentSetup()->GetCurrEngIndex()) << ")." << endl;
2864 str1 << _("MD steps ") << param.nsteps_h << "/" << param.nsteps_e << "/" << param.nsteps_s << "/" << param.nsteps_c << " ";
2865 str1 << "ts= " << param.timestep << " ";
2866 str1 << "T= " << param.target_T << " (rt " << param.T_rtime_hc << " " << param.T_rtime_es << ") ";
2867 str1 << "P= " << param.target_P << " (rt " << param.P_rtime << " beta " << param.P_beta << ") ";
2868 str1 << "cT= " << param.constant_T << " cP= " << param.constant_P << " ";
2869 str1 << endl << ends;
2870
2871 PrintToLog(str1.str().c_str());
2872
2873 CopyCRD(this, eng, 0);
2874
2875 moldyn * dyn = new moldyn(eng, param.timestep);
2876
2877 char logfilename[256];
2878 strcpy(logfilename, param.filename);
2879 i32s last_dot = NOT_DEFINED;
2880 for (i32u fn = 0;fn < strlen(logfilename);fn++)
2881 {
2882 if (logfilename[fn] == '.') last_dot = (i32s) fn;
2883 } if (last_dot < 0) last_dot = strlen(logfilename);
2884 logfilename[last_dot + 0] = '.'; logfilename[last_dot + 1] = 'l';
2885 logfilename[last_dot + 2] = 'o'; logfilename[last_dot + 3] = 'g';
2886 logfilename[last_dot + 4] = 0;
2887
2888 ofstream logfile;
2889 logfile.open(logfilename, ios::out);
2890
2891 ofstream ofile; // the trajectory file...
2892 ofile.open(param.filename, ios::out | ios::binary);
2893
2894 const int frame_save_frq = 100;
2895 const int total_frames = param.nsteps_s / frame_save_frq;
2896
2897 WriteTrajectoryHeader(ofile, total_frames);
2898
2899 ThreadUnlock();
2900
2901 double sum_pressure = 0.0;
2902 int sum_counter = 0;
2903
2904 if (param.constant_P)
2905 {
2906 dyn->target_pressure = param.target_P;
2907 dyn->pressure_rtime = param.P_rtime;
2908 dyn->isoth_compr = param.P_beta;
2909 }
2910
2911 const i32s tot_nsteps = param.nsteps_h + param.nsteps_e + param.nsteps_s + param.nsteps_c;
2912
2913 for (i32s n1 = 0;n1 < tot_nsteps;n1++)
2914 {
2915 if (!(n1 % 10))
2916 {
2917 if (eng_pbc != NULL) eng_pbc->CheckLocations();
2918 eng->RequestNeighborListUpdate();
2919 }
2920
2921 // check if entering heating stage.
2922 if (n1 == 0)
2923 {
2924 dyn->temperature_rtime = param.T_rtime_hc;
2925 }
2926
2927 // check if entering equilibration stage.
2928 if (n1 == param.nsteps_h)
2929 {
2930 dyn->target_temperature = param.target_T;
2931 dyn->temperature_rtime = param.T_rtime_es;
2932 }
2933
2934 // check if entering cooling stage.
2935 if (n1 == param.nsteps_h + param.nsteps_e + param.nsteps_s)
2936 {
2937 dyn->temperature_rtime = param.T_rtime_hc;
2938 }
2939
2940 // check if T adjustment is needed at heating stage.
2941 if (n1 < param.nsteps_h && !(n1 % 50))
2942 {
2943 int tmp1 = n1;
2944 int tmp2 = param.nsteps_h;
2945
2946 dyn->target_temperature = param.target_T * ((f64) tmp1 / (f64) tmp2);
2947 cout << _("setting T = ") << dyn->target_temperature << endl;
2948 }
2949
2950 // check if T adjustment is needed at cooling stage.
2951 if (n1 >= param.nsteps_h + param.nsteps_e + param.nsteps_s && !(n1 % 50))
2952 {
2953 int tmp1 = n1 - (param.nsteps_h + param.nsteps_e + param.nsteps_s);
2954 int tmp2 = param.nsteps_c;
2955
2956 dyn->target_temperature = param.target_T * ((f64) (tmp2 - tmp1) / (f64) tmp2);
2957 cout << _("setting T = ") << dyn->target_temperature << endl;
2958 }
2959
2960 bool enable_Tc = false;
2961 if (n1 < param.nsteps_h + param.nsteps_e) enable_Tc = true; // heating/equilibration
2962 if (n1 >= param.nsteps_h + param.nsteps_e + param.nsteps_s) enable_Tc = true; // cooling
2963 if (param.constant_T) enable_Tc = true;
2964
2965 bool enable_Pc = param.constant_P;
2966 if (n1 < param.nsteps_h) enable_Pc = false; // heating
2967 if (n1 >= param.nsteps_h + param.nsteps_e + param.nsteps_s) enable_Pc = false; // cooling
2968
2969 dyn->TakeMDStep(enable_Tc, enable_Pc);
2970
2971 if (enable_Pc)
2972 {
2973 sum_pressure += dyn->saved_pressure;
2974 sum_counter++;
2975
2976 if (sum_counter >= 500)
2977 {
2978 double aP = sum_pressure / (f64) sum_counter;
2979
2980 ThreadLock();
2981
2982 ostringstream str2;
2983 str2 << _("pressure ") << aP << " bar ; ";
2984 str2 << _("density ") << dyn->saved_density << " kg/dm^3 ";
2985 str2 << endl << ends;
2986
2987 PrintToLog(str2.str().c_str());
2988 logfile << str2.str().c_str();
2989
2990 ThreadUnlock();
2991
2992 sum_pressure = 0.0;
2993 sum_counter = 0;
2994 }
2995 }
2996
2997 // check if log/progressbar output should be made.
2998 if (!(n1 % 100))
2999 {
3000 ThreadLock();
3001
3002 double progress = (double) (n1 + 1) / (double) (tot_nsteps);
3003 bool cancel = SetProgress(progress, NULL);
3004
3005 if (cancel)
3006 {
3007 ThreadUnlock();
3008 break;
3009 }
3010
3011 ostringstream str2a;
3012 str2a << _("step ") << n1 << " T = " << dyn->ConvEKinTemp(dyn->GetEKin()) << " K ";
3013 str2a << _("Epot = ") << dyn->GetEPot() << _(" kJ/mol Etot = ") << (dyn->GetEKin() + dyn->GetEPot()) << " kJ/mol ";
3014 str2a << endl << ends;
3015
3016 PrintToLog(str2a.str().c_str());
3017
3018 ostringstream str2b;
3019 str2b << _("step ") << n1 << " T = " << dyn->ConvEKinTemp(dyn->GetEKin()) << " ";
3020 str2b << _("Epot = ") << dyn->GetEPot() << _(" Etot = ") << (dyn->GetEKin() + dyn->GetEPot()) << " ;; ";
3021 str2b << endl << ends;
3022
3023 logfile << str2b.str().c_str();
3024
3025 ThreadUnlock();
3026
3027 NoThreadsIterate(); // this will update the GUI when threads disabled...
3028 }
3029
3030 if (!(n1 % 1000))
3031 {
3032 if (eng_bp != NULL && eng_bp->nd_eval != NULL)
3033 {
3034 ThreadLock();
3035
3036 ostringstream str3;
3037 eng_bp->nd_eval->PrintResults(str3); str3 << ends;
3038
3039 PrintToLog(str3.str().c_str());
3040 logfile << str3.str().c_str();
3041
3042 ThreadUnlock();
3043 }
3044
3045 if (eng_bp != NULL && eng_bp->rdf_eval != NULL)
3046 {
3047 ThreadLock();
3048
3049 ostringstream str3;
3050 eng_bp->rdf_eval->PrintResults(str3); str3 << ends;
3051
3052 PrintToLog(str3.str().c_str());
3053 logfile << str3.str().c_str();
3054
3055 ThreadUnlock();
3056 }
3057 }
3058
3059 // check if trajectory output should be made.
3060 if (!(n1 < param.nsteps_h + param.nsteps_e) && !(n1 % frame_save_frq))
3061 {
3062 CopyCRD(eng, this, 0);
3063 WriteTrajectoryFrame(ofile, dyn);
3064 }
3065
3066 // check if graphics update should be made.
3067 if (!(n1 % 100))
3068 {
3069 ThreadLock();
3070
3071 CopyCRD(eng, this, 0);
3072 UpdateAllGraphicsViews(updt);
3073
3074 ThreadUnlock();
3075
3076 // the update frequency here matches with log/progressbar frequency,
3077 // so that there is no need to call NoThreadsIterate() again...
3078 }
3079 }
3080
3081 ofile.close();
3082 logfile.close();
3083
3084 delete dyn;
3085
3086 // we will not delete current_eng here, so that we can draw plots using it...
3087 // we will not delete current_eng here, so that we can draw plots using it...
3088 // we will not delete current_eng here, so that we can draw plots using it...
3089
3090 // above, CopyCRD was done eng->mdl and then CenterCRDSet() was done for mdl.
3091 // this might cause that old coordinates remain in eng object, possibly affecting plots.
3092 // here we sync the coordinates and other plotting data in the eng object.
3093
3094 ThreadLock();
3095 CopyCRD(this, eng, 0);
3096 SetupPlotting();
3097 ThreadUnlock();
3098 }
3099
DoRandomSearch(i32s cycles,i32s optsteps,bool updt)3100 void model::DoRandomSearch(i32s cycles, i32s optsteps, bool updt)
3101 {
3102 // make this thread-safe since this can be called from project::process_job_WhatEver() at the app side...
3103 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
3104 // this means: 1) make sure that PrintToLog() is safe, 2) only call UpdateAllGraphicsViews(false) because
3105 // OGL-context is owned by the GUI-thread ; calling with "false" will pass the update to the GUI-thread.
3106
3107 ThreadLock();
3108
3109 // eng is not really needed here, but just check that it's available...
3110 // see also project::DoEnergyPlot1D and project::DoEnergyPlot2D
3111 engine * eng = GetCurrentSetup()->GetCurrentEngine();
3112 if (eng == NULL) GetCurrentSetup()->CreateCurrentEngine();
3113 eng = GetCurrentSetup()->GetCurrentEngine();
3114 if (eng == NULL)
3115 {
3116 ThreadUnlock();
3117 return;
3118 }
3119
3120 if (cs_vector.size() < 2)
3121 {
3122 PushCRDSets(1);
3123 SetCRDSetVisible(1, false);
3124 }
3125
3126 random_search rs(this, 0, 0, 1, cycles, optsteps);
3127
3128 ThreadUnlock();
3129
3130 bool cancel = false;
3131 while (!cancel)
3132 {
3133 int zzz = rs.TakeStep();
3134 //cout << zzz << endl;
3135
3136 if (rs.last_step != NOT_DEFINED)
3137 {
3138 ThreadLock();
3139
3140 stringstream str1;
3141 str1 << _("step ") << rs.last_step << "/" << cycles << _(" energy = ") << rs.last_E << " kJ/mol" << endl << ends;
3142 PrintToLog(str1.str().c_str());
3143
3144 double progress = (double) rs.last_step / (double) cycles;
3145 cancel = SetProgress(progress, NULL);
3146
3147 ThreadUnlock();
3148
3149 NoThreadsIterate(); // this will update the GUI when threads disabled...
3150 }
3151
3152 static int updatefrq = 0;
3153 if (!(updatefrq % 10))
3154 {
3155 ThreadLock();
3156
3157 UpdateAllGraphicsViews(updt);
3158
3159 ThreadUnlock();
3160
3161 NoThreadsIterate(); // this will update the GUI when threads disabled...
3162 }
3163
3164 if (zzz < 0) break;
3165 }
3166
3167 ThreadLock();
3168
3169 CopyCRDSet(1, 0); // get the best structure!!!
3170 PopCRDSets(1); // remove the extra crd-sets.
3171
3172 DiscardCurrEng(); // the coordinates are changed -> disable plotting.
3173 //SetupPlotting(); // is this an alternative??? need to check this...
3174
3175 UpdateAllGraphicsViews(updt);
3176
3177 stringstream str1;
3178 str1 << _("lowest energy found = ") << rs.GetMinEnergy() << " kJ/mol" << endl << ends;
3179 PrintToLog(str1.str().c_str());
3180
3181 ostringstream oss2;
3182 oss2 << _("RANDOM SEARCH is ready");
3183 if (cancel) oss2 << _(" (cancelled)");
3184 oss2 << "." << endl << ends;
3185
3186 PrintToLog(oss2.str().c_str());
3187
3188 ThreadUnlock();
3189 }
3190
DoSystematicSearch(i32s divisions,i32s optsteps,bool)3191 void model::DoSystematicSearch(i32s divisions, i32s optsteps, bool)
3192 {
3193 // eng is not really needed here, but just check that it's available...
3194 // see also project::DoEnergyPlot1D and project::DoEnergyPlot2D
3195 engine * eng = GetCurrentSetup()->GetCurrentEngine();
3196 if (eng == NULL) GetCurrentSetup()->CreateCurrentEngine();
3197 eng = GetCurrentSetup()->GetCurrentEngine();
3198 if (eng == NULL) return;
3199
3200 if (cs_vector.size() < 2)
3201 {
3202 PushCRDSets(1);
3203 SetCRDSetVisible(1, false);
3204 }
3205
3206 systematic_search ss(this, 0, 0, 1, divisions, optsteps);
3207
3208 while (true)
3209 {
3210 int zzz = ss.TakeStep();
3211 // cout << zzz << endl;
3212
3213 UpdateAllGraphicsViews(true);
3214 if (zzz < 0) break;
3215 }
3216
3217 CopyCRDSet(1, 0); // get the best structure!!!
3218 PopCRDSets(1); // remove the extra crd-sets.
3219
3220 DiscardCurrEng(); // the coordinates are changed -> disable plotting.
3221 //SetupPlotting(); // is this an alternative??? need to check this...
3222
3223 UpdateAllGraphicsViews(true);
3224
3225 stringstream str1;
3226 str1 << _("lowest energy found = ") << ss.GetMinEnergy() << " kJ/mol" << endl << ends;
3227 PrintToLog(str1.str().c_str());
3228
3229 ostringstream oss2;
3230 oss2 << _("SYSTEMATIC SEARCH is ready");
3231 // if (cancel) oss2 << _(" (cancelled)");
3232 oss2 << "." << endl << ends;
3233
3234 PrintToLog(oss2.str().c_str());
3235 }
3236
DoMonteCarloSearch(i32s n_init_steps,i32s n_simul_steps,i32s optsteps,bool)3237 void model::DoMonteCarloSearch(i32s n_init_steps, i32s n_simul_steps, i32s optsteps, bool)
3238 {
3239 // eng is not really needed here, but just check that it's available...
3240 // see also project::DoEnergyPlot1D and project::DoEnergyPlot2D
3241 engine * eng = GetCurrentSetup()->GetCurrentEngine();
3242 if (eng == NULL) GetCurrentSetup()->CreateCurrentEngine();
3243 eng = GetCurrentSetup()->GetCurrentEngine();
3244 if (eng == NULL) return;
3245
3246
3247 if (cs_vector.size() < 2)
3248 {
3249 PushCRDSets(1);
3250 SetCRDSetVisible(1, false);
3251 }
3252
3253 monte_carlo_search mcs(this, 0, 0, 1, n_init_steps, n_simul_steps, optsteps);
3254
3255 while (true)
3256 {
3257 int zzz = mcs.TakeStep();
3258 // cout << zzz << endl;
3259
3260 UpdateAllGraphicsViews(true);
3261 if (zzz < 0) break;
3262 }
3263
3264 CopyCRDSet(1, 0); // get the best structure!!!
3265 PopCRDSets(1); // remove the extra crd-sets.
3266
3267 DiscardCurrEng(); // the coordinates are changed -> disable plotting.
3268 //SetupPlotting(); // is this an alternative??? need to check this...
3269
3270 UpdateAllGraphicsViews(true);
3271
3272 stringstream str1;
3273 str1 << _("lowest energy found = ") << mcs.GetMinEnergy() << " kJ/mol" << endl << ends;
3274 PrintToLog(str1.str().c_str());
3275
3276 ostringstream oss2;
3277 oss2 << _("MONTE CARLO SEARCH is ready");
3278 // if (cancel) oss2 << _(" (cancelled)");
3279 oss2 << "." << endl << ends;
3280
3281 PrintToLog(oss2.str().c_str());
3282 }
3283
3284 /*##############################################*/
3285 /*##############################################*/
3286
3287 // what comes to PDB chain id's here, we follow the "pdb-select principle" by converting
3288 // the empty id ' ' with '_' to improve readability...
3289
readpdb_ReadMData(const char * filename)3290 readpdb_mdata * model::readpdb_ReadMData(const char * filename)
3291 {
3292 cout << _("reading PDB metadata from file ") << filename << endl;
3293
3294 readpdb_mdata * mdata = new readpdb_mdata;
3295
3296 char buffer[1024];
3297
3298 ifstream file(filename, ios::in);
3299 if (file.fail())
3300 {
3301 cout << _("file \"") << filename << _("\" not found.") << endl;
3302
3303 file.close();
3304 return mdata;
3305 }
3306
3307 while (!file.eof())
3308 {
3309 char record_id[8];
3310 for (i32s n1 = 0;n1 < 6;n1++)
3311 {
3312 char tchar = (char) file.get();
3313 if (tchar != ' ') record_id[n1] = tchar;
3314 else record_id[n1] = 0;
3315 }
3316
3317 record_id[6] = 0; // terminate the record string...
3318
3319 if (!strcmp("SEQRES", record_id))
3320 {
3321 i32s blocknum; file >> blocknum;
3322 file.get(); char chn_id = (char) file.get();
3323 i32s chn_length; file >> chn_length;
3324
3325 if (chn_id == ' ') chn_id = '_'; // fix the empty chain id...
3326
3327 i32s chn_num;
3328 if (blocknum == 1) // this is a new record...
3329 {
3330 cout << _("found a new chain ") << mdata->chn_vector.size();
3331 cout << " '" << chn_id << _("' with ") << chn_length << _(" residues.") << endl;
3332
3333 readpdb_mdata_chain * new_data = new readpdb_mdata_chain;
3334 new_data->chn_id = chn_id; new_data->seqres = new char[chn_length + 1];
3335
3336 for (i32s n1 = 0;n1 < chn_length;n1++) new_data->seqres[n1] = '?';
3337 new_data->seqres[chn_length] = 0; // terminate the string!!!
3338
3339 chn_num = mdata->chn_vector.size();
3340 mdata->chn_vector.push_back(new_data);
3341 }
3342 else // this is an old record, find it...
3343 {
3344 chn_num = 0;
3345 while (chn_num < ((i32s) mdata->chn_vector.size()))
3346 {
3347 if (mdata->chn_vector[chn_num]->chn_id == chn_id) break;
3348 else chn_num++;
3349 }
3350
3351 if (chn_num == ((i32s) mdata->chn_vector.size()))
3352 {
3353 assertion_failed(__FILE__, __LINE__, "readpdb_ReadMData : unknown chain found!");
3354 }
3355 }
3356
3357 chn_length = strlen(mdata->chn_vector[chn_num]->seqres);
3358 i32s counter = (blocknum - 1) * 13;
3359
3360 for (i32s n1 = 0;n1 < 13;n1++)
3361 {
3362 char residue[16];
3363 file >> residue;
3364
3365 char symbol = '?';
3366 if (!strcmp("ALA", residue)) symbol = 'A';
3367 if (!strcmp("ARG", residue)) symbol = 'R';
3368 if (!strcmp("ASN", residue)) symbol = 'N';
3369 if (!strcmp("ASP", residue)) symbol = 'D';
3370 if (!strcmp("CYS", residue)) symbol = 'C';
3371 if (!strcmp("GLN", residue)) symbol = 'Q';
3372 if (!strcmp("GLU", residue)) symbol = 'E';
3373 if (!strcmp("GLY", residue)) symbol = 'G';
3374 if (!strcmp("HIS", residue)) symbol = 'H';
3375 if (!strcmp("ILE", residue)) symbol = 'I';
3376 if (!strcmp("LEU", residue)) symbol = 'L';
3377 if (!strcmp("LYS", residue)) symbol = 'K';
3378 if (!strcmp("MET", residue)) symbol = 'M';
3379 if (!strcmp("PHE", residue)) symbol = 'F';
3380 if (!strcmp("PRO", residue)) symbol = 'P';
3381 if (!strcmp("SER", residue)) symbol = 'S';
3382 if (!strcmp("THR", residue)) symbol = 'T';
3383 if (!strcmp("TRP", residue)) symbol = 'W';
3384 if (!strcmp("TYR", residue)) symbol = 'Y';
3385 if (!strcmp("VAL", residue)) symbol = 'V';
3386
3387 mdata->chn_vector[chn_num]->seqres[counter++] = symbol;
3388 if (counter == chn_length) break;
3389 }
3390 }
3391
3392 file.getline(buffer, sizeof(buffer)); // move to the next line...
3393 }
3394
3395 file.close();
3396
3397 // ready...
3398
3399 if (mdata->chn_vector.empty()) cout << _("WARNING : no chains found!!!") << endl;
3400 cout << _("done.") << endl;
3401
3402 return mdata;
3403 }
3404
3405 // chn_num must be either NOT_DEFINED (which is -1 meaning all chains) or the chain index in mdata...
3406 // chn_num must be either NOT_DEFINED (which is -1 meaning all chains) or the chain index in mdata...
3407 // chn_num must be either NOT_DEFINED (which is -1 meaning all chains) or the chain index in mdata...
3408
3409 //#define READPDB_ENABLE_MULTIPLE_CRDSETS // enable this to read all data from NMR entries...
3410
readpdb_ReadData(const char * filename,readpdb_mdata * mdata,i32s chn_num)3411 void model::readpdb_ReadData(const char * filename, readpdb_mdata * mdata, i32s chn_num)
3412 {
3413 cout << _("reading PDB data from file ") << filename << endl;
3414
3415 ifstream file(filename, ios::in);
3416 if (file.fail())
3417 {
3418 cout << _("file \"") << filename << _("\" not found.") << endl;
3419
3420 file.close();
3421 return;
3422 }
3423
3424 vector<atom *> all_added_atoms;
3425
3426 // read the whole file to temporary arrays...
3427
3428 vector<readpdb_data_atom> atom_data;
3429 vector<readpdb_data_ssbond> ssbond_data;
3430
3431 i32s model_counter = 0;
3432 i32s atom_counter = NOT_DEFINED;
3433
3434 while (!file.eof())
3435 {
3436 char record_id[8];
3437 for (i32s n1 = 0;n1 < 6;n1++)
3438 {
3439 char tchar = (char) file.get();
3440 if (tchar != ' ') record_id[n1] = tchar;
3441 else record_id[n1] = 0;
3442 }
3443
3444 record_id[6] = 0; // terminate the record string...
3445
3446 if (!strcmp("ATOM", record_id))
3447 {
3448 i32s serial_number; file >> serial_number;
3449
3450 char buffer[20];
3451 for (i32s n1 = 0;n1 < 18;n1++) buffer[n1] = (char) file.get();
3452 buffer[5] = ' '; buffer[15] = ' ';
3453
3454 char chn_id; i32s res_num; char res_name[5]; char atm_name[5]; fGL crd[3];
3455
3456 istringstream str(buffer);
3457 str >> atm_name >> res_name; str.get(); chn_id = (char) str.get(); str >> res_num;
3458 file >> crd[0]; crd[0] /= 10.0; file >> crd[1]; crd[1] /= 10.0; file >> crd[2]; crd[2] /= 10.0;
3459
3460 if (chn_id == ' ') chn_id = '_'; // fix the empty chain id...
3461
3462 bool test1 = (chn_num == NOT_DEFINED);
3463 bool test2 = test1 ? false : (chn_id == mdata->chn_vector[chn_num]->chn_id);
3464
3465 if (test1 || test2)
3466 {
3467 if (model_counter == 0)
3468 {
3469 readpdb_data_atom new_data;
3470 new_data.ref = NULL;
3471
3472 strcpy(new_data.atm_name, atm_name);
3473 strcpy(new_data.res_name, res_name);
3474
3475 new_data.chn_id = chn_id;
3476 new_data.res_num = res_num;
3477
3478 new_data.crd[model_counter][0] = crd[0];
3479 new_data.crd[model_counter][1] = crd[1];
3480 new_data.crd[model_counter][2] = crd[2];
3481
3482 atom_data.push_back(new_data);
3483 }
3484 else
3485 {
3486 bool test1 = !strcmp(atom_data[atom_counter].res_name, res_name);
3487 bool test2 = !strcmp(atom_data[atom_counter].atm_name, atm_name);
3488
3489 if (!test1) assertion_failed(__FILE__, __LINE__, "res_name");
3490 if (!test2) assertion_failed(__FILE__, __LINE__, "atm_name");
3491
3492 atom_data[atom_counter].crd[model_counter][0] = crd[0];
3493 atom_data[atom_counter].crd[model_counter][1] = crd[1];
3494 atom_data[atom_counter].crd[model_counter][2] = crd[2];
3495
3496 atom_counter++;
3497 }
3498 }
3499 }
3500
3501 if (!strcmp("SSBOND", record_id))
3502 {
3503 i32s serial_number; file >> serial_number;
3504 readpdb_data_ssbond new_data; new_data.ref = NULL;
3505
3506 for (i32s n1 = 0;n1 < 2;n1++)
3507 {
3508 char buffer[20];
3509 file >> buffer; file.get();
3510 char chn_id = (char) file.get();
3511
3512 if (chn_id == ' ') chn_id = '_'; // fix the empty chain id...
3513
3514 new_data.chn_id = chn_id;
3515
3516 for (i32s n2 = 0;n2 < 6;n2++) buffer[n2] = (char) file.get();
3517 istringstream str(buffer); str >> new_data.res_num;
3518
3519 ssbond_data.push_back(new_data);
3520 }
3521 }
3522
3523 // how to deal with multiple coordinate sets (often present in NMR entries)???
3524 // how to deal with multiple coordinate sets (often present in NMR entries)???
3525 // how to deal with multiple coordinate sets (often present in NMR entries)???
3526
3527 #ifdef READPDB_ENABLE_MULTIPLE_CRDSETS
3528
3529 // in this caseat the moment, we will read multiple records up to READPDB_MAX_CRDSETS.
3530 // if READPDB_MAX_CRDSETS is exceeded, the program stops; must recompile then...
3531
3532 if (!strcmp("MODEL", record_id))
3533 {
3534 i32s model_number; file >> model_number;
3535
3536 if (model_number > 1)
3537 {
3538 model_counter++;
3539
3540 if (model_counter == READPDB_MAX_CRDSETS) // this should never happen...
3541 {
3542 assertion_failed(__FILE__, __LINE__, "READPDB_MAX_CRDSETS exceeded!");
3543 }
3544
3545 atom_counter = 0;
3546 }
3547 }
3548
3549 #else // READPDB_ENABLE_MULTIPLE_CRDSETS
3550
3551 // in this case, we just look for "ENDMDL" record, which says this coordinate set is
3552 // now ready. so, we read only the first coordinate set. the "MODEL"/"ENDMDL" records
3553 // are only after "SSBOND", so it should not have big effects...
3554
3555 if (!strcmp("ENDMDL", record_id))
3556 {
3557 cout << _("ENDMDL record found, skipping the rest of this file...")<< endl;
3558 break;
3559 }
3560
3561 #endif // READPDB_ENABLE_MULTIPLE_CRDSETS
3562
3563 char buffer[1024];
3564 file.getline(buffer, sizeof(buffer)); // move to the next line...
3565 }
3566
3567 file.close();
3568
3569 if (atom_data.empty())
3570 {
3571 cout << _("no atoms found!!!") << endl;
3572 return;
3573 }
3574
3575 // we have now read all the crd-sets; check if we need more space for them...
3576
3577 i32s csets = cs_vector.size(); i32s new_csets = ((model_counter + 1) - csets);
3578 cout << _("there were ") << csets << _(" old crd-sets, creating ") << new_csets << _(" new...") << endl;
3579
3580 PushCRDSets(new_csets);
3581
3582 for (i32u n1 = 0;n1 < cs_vector.size();n1++) SetCRDSetVisible(n1, true);
3583
3584 // count the chains and relate them to mdata, create new records to mdata if necessary...
3585
3586 vector<i32s> chn_index;
3587 char tmp_chn_id; i32u tmp_index;
3588
3589 tmp_chn_id = atom_data.front().chn_id;
3590
3591 tmp_index = 0;
3592 while (tmp_index < mdata->chn_vector.size())
3593 {
3594 if (tmp_chn_id == mdata->chn_vector[tmp_index]->chn_id) break;
3595 else tmp_index++;
3596 }
3597
3598 if (tmp_index == mdata->chn_vector.size())
3599 {
3600 readpdb_mdata_chain * unknown = new readpdb_mdata_chain;
3601 unknown->chn_id = tmp_chn_id; unknown->seqres = NULL;
3602 mdata->chn_vector.push_back(unknown);
3603 }
3604
3605 chn_index.push_back(tmp_index);
3606
3607 for (i32u n1 = 1;n1 < atom_data.size();n1++)
3608 {
3609 tmp_chn_id = atom_data[n1].chn_id;
3610 if (tmp_chn_id == mdata->chn_vector[chn_index.back()]->chn_id) continue;
3611 else
3612 {
3613 tmp_index = 0;
3614 while (tmp_index < mdata->chn_vector.size())
3615 {
3616 if (tmp_chn_id == mdata->chn_vector[tmp_index]->chn_id) break;
3617 else tmp_index++;
3618 }
3619
3620 if (tmp_index == mdata->chn_vector.size())
3621 {
3622 readpdb_mdata_chain * unknown = new readpdb_mdata_chain;
3623 unknown->chn_id = tmp_chn_id; unknown->seqres = NULL;
3624 mdata->chn_vector.push_back(unknown);
3625 }
3626
3627 chn_index.push_back(tmp_index);
3628 }
3629 }
3630
3631 // now create the chains, checking validity of the residues...
3632
3633 for (i32u n1 = 0;n1 < chn_index.size();n1++)
3634 {
3635 char current_chn_id = mdata->chn_vector[chn_index[n1]]->chn_id;
3636
3637 i32s current_chn_length;
3638 if (mdata->chn_vector[chn_index[n1]]->seqres != NULL)
3639 {
3640 current_chn_length = strlen(mdata->chn_vector[chn_index[n1]]->seqres);
3641 }
3642 else
3643 {
3644 current_chn_length = NOT_DEFINED;
3645 }
3646
3647 i32s previous_residue = 0;
3648
3649 i32s range1[2];
3650
3651 range1[0] = 0;
3652 while (atom_data[range1[0]].chn_id != current_chn_id) range1[0]++;
3653
3654 range1[1] = range1[0];
3655 while (range1[1] < (i32s) atom_data.size() && atom_data[range1[1]].chn_id == current_chn_id) range1[1]++;
3656
3657 vector<i32s> res_data;
3658
3659 i32s range2[2];
3660 range2[0] = range1[0];
3661 while (range2[0] < range1[1])
3662 {
3663 i32s residue = atom_data[range2[0]].res_num;
3664 const char * res_name = atom_data[range2[0]].res_name;
3665
3666 range2[1] = range2[0];
3667 while (range2[1] < (i32s) atom_data.size() && atom_data[range2[1]].res_num == residue) range2[1]++;
3668
3669 bool standard = false;
3670 if (!strcmp("ALA", res_name)) standard = true;
3671 if (!strcmp("ARG", res_name)) standard = true;
3672 if (!strcmp("ASN", res_name)) standard = true;
3673 if (!strcmp("ASP", res_name)) standard = true;
3674 if (!strcmp("CYS", res_name)) standard = true;
3675 if (!strcmp("GLN", res_name)) standard = true;
3676 if (!strcmp("GLU", res_name)) standard = true;
3677 if (!strcmp("GLY", res_name)) standard = true;
3678 if (!strcmp("HIS", res_name)) standard = true;
3679 if (!strcmp("ILE", res_name)) standard = true;
3680 if (!strcmp("LEU", res_name)) standard = true;
3681 if (!strcmp("LYS", res_name)) standard = true;
3682 if (!strcmp("MET", res_name)) standard = true;
3683 if (!strcmp("PHE", res_name)) standard = true;
3684 if (!strcmp("PRO", res_name)) standard = true;
3685 if (!strcmp("SER", res_name)) standard = true;
3686 if (!strcmp("THR", res_name)) standard = true;
3687 if (!strcmp("TRP", res_name)) standard = true;
3688 if (!strcmp("TYR", res_name)) standard = true;
3689 if (!strcmp("VAL", res_name)) standard = true;
3690
3691 bool skip = false;
3692 if (standard)
3693 {
3694 if (readpdb_ReadData_sub1(atom_data, range2, "N", true) == NOT_DEFINED) skip = true;
3695 if (readpdb_ReadData_sub1(atom_data, range2, "CA", true) == NOT_DEFINED) skip = true;
3696 if (readpdb_ReadData_sub1(atom_data, range2, "C", true) == NOT_DEFINED) skip = true;
3697 if (readpdb_ReadData_sub1(atom_data, range2, "O", true) == NOT_DEFINED) skip = true;
3698 }
3699 else
3700 {
3701 cout << _("could not recognize this residue: ") << res_name << endl;
3702 skip = true;
3703 }
3704
3705 if (skip)
3706 {
3707 cout << _("skipping broken residue ") << residue << " " << res_name << endl;
3708 atom_data.erase(atom_data.begin() + range2[0], atom_data.begin() + range2[1]);
3709 range1[1] += (range2[0] - range2[1]);
3710 }
3711 else
3712 {
3713 // ok, there seems to be no big problems with this residue.
3714 // first we make some bookkeeping with the original sequence...
3715
3716 if ((previous_residue + 1) != residue)
3717 {
3718 for (i32s n2 = (previous_residue + 1);n2 < residue;n2++)
3719 {
3720 // the variables use "pascal/fortran-style" counting, beginning 1,
3721 // but we will store the values in "c-style" way, beginning 0!!!
3722
3723 mdata->chn_vector[chn_index[n1]]->missing_residues.push_back(n2 - 1);
3724 }
3725 }
3726
3727 previous_residue = residue;
3728 bool alpha_carbon_is_found = false;
3729
3730 // ...and then we just create the atoms.
3731
3732 res_data.push_back(range2[0]);
3733 for (i32s n2 = range2[0];n2 < range2[1];n2++)
3734 {
3735 char buffer[32];
3736 strcpy(buffer, atom_data[n2].atm_name);
3737 buffer[1] = 0; // suppose single-char elements!!!
3738
3739 element el = element(buffer);
3740
3741 if (el.GetAtomicNumber() == 1) continue;
3742 atom newatom(el, NULL, cs_vector.size());
3743
3744 // since 20060809 ; determine formal charges!!!
3745 // since 20060809 ; determine formal charges!!!
3746 //////////////////////////////////////////////////
3747
3748 // if the data here seems to be in conflict with builder/amino.txt file, then change this...
3749 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
3750
3751 // check if we are at N-terminus.
3752 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
3753 // alpha_carbons.size() and residue are equal???
3754 if (!mdata->chn_vector[chn_index[n1]]->alpha_carbons.size())
3755 {
3756 if (!strcmp(atom_data[n2].atm_name, "N"))
3757 {
3758 newatom.formal_charge = +1;
3759 }
3760 }
3761
3762 // there is no check for C-terminus since it is handled later...
3763 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
3764
3765 if (!strcmp(atom_data[n2].res_name, "ARG"))
3766 {
3767 if (!strcmp(atom_data[n2].atm_name, "NH1"))
3768 {
3769 newatom.formal_charge = +1;
3770 }
3771 }
3772
3773 if (!strcmp(atom_data[n2].res_name, "ASP"))
3774 {
3775 if (!strcmp(atom_data[n2].atm_name, "OD2"))
3776 {
3777 newatom.formal_charge = -1;
3778 }
3779 }
3780
3781 if (!strcmp(atom_data[n2].res_name, "GLU"))
3782 {
3783 if (!strcmp(atom_data[n2].atm_name, "OE2"))
3784 {
3785 newatom.formal_charge = -1;
3786 }
3787 }
3788
3789 if (!strcmp(atom_data[n2].res_name, "LYS"))
3790 {
3791 if (!strcmp(atom_data[n2].atm_name, "NZ"))
3792 {
3793 newatom.formal_charge = +1;
3794 }
3795 }
3796
3797 //////////////////////////////////////////////////
3798 // since 20060809 ; determine formal charges!!!
3799 // since 20060809 ; determine formal charges!!!
3800
3801 for (i32u n3 = 0;n3 < cs_vector.size();n3++)
3802 {
3803 fGL x = atom_data[n2].crd[n3][0];
3804 fGL y = atom_data[n2].crd[n3][1];
3805 fGL z = atom_data[n2].crd[n3][2];
3806 newatom.SetCRD(n3, x, y, z);
3807 }
3808
3809 AddAtom_lg(newatom);
3810 all_added_atoms.push_back(& atom_list.back());
3811
3812 atom_data[n2].ref = (& atom_list.back());
3813
3814 // if this atom was an alpha-carbon, then add it to the list... but alternative
3815 // locations might cause trouble here; make sure that only one c-alpha is added!!!
3816
3817 if (!strcmp(atom_data[n2].atm_name, "CA") && !alpha_carbon_is_found)
3818 {
3819 mdata->chn_vector[chn_index[n1]]->alpha_carbons.push_back(atom_data[n2].ref);
3820 alpha_carbon_is_found = true;
3821 }
3822 }
3823
3824 readpdb_ReadData_sub2(atom_data, range2, "N", "CA", 'S');
3825 readpdb_ReadData_sub2(atom_data, range2, "CA", "C", 'S');
3826 readpdb_ReadData_sub2(atom_data, range2, "C", "O", 'D');
3827
3828 // if the data here seems to be in conflict with builder/amino.txt file, then change this...
3829 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
3830
3831 if (!strcmp("ALA", res_name))
3832 {
3833 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3834 }
3835
3836 if (!strcmp("ARG", res_name))
3837 {
3838 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3839 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
3840 readpdb_ReadData_sub2(atom_data, range2, "CG", "CD", 'S');
3841 readpdb_ReadData_sub2(atom_data, range2, "CD", "NE", 'S');
3842 readpdb_ReadData_sub2(atom_data, range2, "NE", "CZ", 'S');
3843 readpdb_ReadData_sub2(atom_data, range2, "CZ", "NH1", 'D');
3844 readpdb_ReadData_sub2(atom_data, range2, "CZ", "NH2", 'S');
3845 }
3846
3847 if (!strcmp("ASN", res_name))
3848 {
3849 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3850 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
3851 readpdb_ReadData_sub2(atom_data, range2, "CG", "OD1", 'D');
3852 readpdb_ReadData_sub2(atom_data, range2, "CG", "ND2", 'S');
3853 }
3854
3855 if (!strcmp("ASP", res_name))
3856 {
3857 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3858 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
3859 readpdb_ReadData_sub2(atom_data, range2, "CG", "OD1", 'D');
3860 readpdb_ReadData_sub2(atom_data, range2, "CG", "OD2", 'S');
3861 }
3862
3863 if (!strcmp("CYS", res_name))
3864 {
3865 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3866 readpdb_ReadData_sub2(atom_data, range2, "CB", "SG", 'S');
3867 }
3868
3869 if (!strcmp("GLN", res_name))
3870 {
3871 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3872 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
3873 readpdb_ReadData_sub2(atom_data, range2, "CG", "CD", 'S');
3874 readpdb_ReadData_sub2(atom_data, range2, "CD", "OE1", 'D');
3875 readpdb_ReadData_sub2(atom_data, range2, "CD", "NE2", 'S');
3876 }
3877
3878 if (!strcmp("GLU", res_name))
3879 {
3880 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3881 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
3882 readpdb_ReadData_sub2(atom_data, range2, "CG", "CD", 'S');
3883 readpdb_ReadData_sub2(atom_data, range2, "CD", "OE1", 'D');
3884 readpdb_ReadData_sub2(atom_data, range2, "CD", "OE2", 'S');
3885 }
3886
3887 // GLY -> we don't have to do anything in this case...
3888 // GLY -> we don't have to do anything in this case...
3889 // GLY -> we don't have to do anything in this case...
3890
3891 if (!strcmp("HIS", res_name)) // HIE form...
3892 {
3893 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3894 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
3895 readpdb_ReadData_sub2(atom_data, range2, "CG", "ND1", 'S');
3896 readpdb_ReadData_sub2(atom_data, range2, "CG", "CD2", 'D');
3897 readpdb_ReadData_sub2(atom_data, range2, "ND1", "CE1", 'D');
3898 readpdb_ReadData_sub2(atom_data, range2, "CD2", "NE2", 'S');
3899 readpdb_ReadData_sub2(atom_data, range2, "CE1", "NE2", 'S');
3900 }
3901
3902 if (!strcmp("ILE", res_name))
3903 {
3904 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3905 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG1", 'S');
3906 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG2", 'S');
3907 readpdb_ReadData_sub2(atom_data, range2, "CG1", "CD1", 'S');
3908 }
3909
3910 if (!strcmp("LEU", res_name))
3911 {
3912 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3913 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
3914 readpdb_ReadData_sub2(atom_data, range2, "CG", "CD1", 'S');
3915 readpdb_ReadData_sub2(atom_data, range2, "CG", "CD2", 'S');
3916 }
3917
3918 if (!strcmp("LYS", res_name))
3919 {
3920 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3921 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
3922 readpdb_ReadData_sub2(atom_data, range2, "CG", "CD", 'S');
3923 readpdb_ReadData_sub2(atom_data, range2, "CD", "CE", 'S');
3924 readpdb_ReadData_sub2(atom_data, range2, "CE", "NZ", 'S');
3925 }
3926
3927 if (!strcmp("MET", res_name))
3928 {
3929 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3930 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
3931 readpdb_ReadData_sub2(atom_data, range2, "CG", "SD", 'S');
3932 readpdb_ReadData_sub2(atom_data, range2, "SD", "CE", 'S');
3933 }
3934
3935 if (!strcmp("PHE", res_name))
3936 {
3937 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3938 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
3939 readpdb_ReadData_sub2(atom_data, range2, "CG", "CD1", 'D');
3940 readpdb_ReadData_sub2(atom_data, range2, "CG", "CD2", 'S');
3941 readpdb_ReadData_sub2(atom_data, range2, "CD1", "CE1", 'S');
3942 readpdb_ReadData_sub2(atom_data, range2, "CD2", "CE2", 'D');
3943 readpdb_ReadData_sub2(atom_data, range2, "CE1", "CZ", 'D');
3944 readpdb_ReadData_sub2(atom_data, range2, "CE2", "CZ", 'S');
3945 }
3946
3947 if (!strcmp("PRO", res_name))
3948 {
3949 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3950 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
3951 readpdb_ReadData_sub2(atom_data, range2, "CG", "CD", 'S');
3952 readpdb_ReadData_sub2(atom_data, range2, "CD", "N", 'S');
3953 }
3954
3955 if (!strcmp("SER", res_name))
3956 {
3957 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3958 readpdb_ReadData_sub2(atom_data, range2, "CB", "OG", 'S');
3959 }
3960
3961 if (!strcmp("THR", res_name))
3962 {
3963 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3964 readpdb_ReadData_sub2(atom_data, range2, "CB", "OG1", 'S');
3965 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG2", 'S');
3966 }
3967
3968 if (!strcmp("TRP", res_name))
3969 {
3970 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3971 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
3972 readpdb_ReadData_sub2(atom_data, range2, "CG", "CD1", 'D');
3973 readpdb_ReadData_sub2(atom_data, range2, "CG", "CD2", 'S');
3974 readpdb_ReadData_sub2(atom_data, range2, "CD1", "NE1", 'S');
3975 readpdb_ReadData_sub2(atom_data, range2, "CD2", "CE2", 'D');
3976 readpdb_ReadData_sub2(atom_data, range2, "CD2", "CE3", 'S');
3977 readpdb_ReadData_sub2(atom_data, range2, "NE1", "CE2", 'S');
3978 readpdb_ReadData_sub2(atom_data, range2, "CE2", "CZ2", 'S');
3979 readpdb_ReadData_sub2(atom_data, range2, "CE3", "CZ3", 'D');
3980 readpdb_ReadData_sub2(atom_data, range2, "CZ2", "CH2", 'D');
3981 readpdb_ReadData_sub2(atom_data, range2, "CZ3", "CH2", 'S');
3982 }
3983
3984 if (!strcmp("TYR", res_name))
3985 {
3986 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
3987 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG", 'S');
3988 readpdb_ReadData_sub2(atom_data, range2, "CG", "CD1", 'D');
3989 readpdb_ReadData_sub2(atom_data, range2, "CG", "CD2", 'S');
3990 readpdb_ReadData_sub2(atom_data, range2, "CD1", "CE1", 'S');
3991 readpdb_ReadData_sub2(atom_data, range2, "CD2", "CE2", 'D');
3992 readpdb_ReadData_sub2(atom_data, range2, "CE1", "CZ", 'D');
3993 readpdb_ReadData_sub2(atom_data, range2, "CE2", "CZ", 'S');
3994 readpdb_ReadData_sub2(atom_data, range2, "CZ", "OH", 'S');
3995 }
3996
3997 if (!strcmp("VAL", res_name))
3998 {
3999 readpdb_ReadData_sub2(atom_data, range2, "CA", "CB", 'S');
4000 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG1", 'S');
4001 readpdb_ReadData_sub2(atom_data, range2, "CB", "CG2", 'S');
4002 }
4003
4004 range2[0] = range2[1];
4005 }
4006 }
4007
4008 // the end of sequence bookkeeping : check if there were missing n-terminal residues.
4009
4010 if (previous_residue < current_chn_length)
4011 {
4012 for (i32s n2 = previous_residue;n2 < current_chn_length;n2++)
4013 {
4014 mdata->chn_vector[chn_index[n1]]->missing_residues.push_back(n2);
4015 }
4016 }
4017
4018 if (!mdata->chn_vector[chn_index[n1]]->missing_residues.empty())
4019 {
4020 char out1 = mdata->chn_vector[chn_index[n1]]->chn_id;
4021 i32u out2 = mdata->chn_vector[chn_index[n1]]->missing_residues.size();
4022
4023 cout << _("at chain '") << out1 << _("' there were ") << out2 << _(" missing residues:") << endl;
4024
4025 for (i32u n2 = 0;n2 < mdata->chn_vector[chn_index[n1]]->missing_residues.size();n2++)
4026 {
4027 cout << mdata->chn_vector[chn_index[n1]]->missing_residues[n2] << " ";
4028 }
4029
4030 cout << endl;
4031 }
4032
4033 // add connecting bonds between residues...
4034
4035 res_data.push_back(range1[1]);
4036 for (i32s n2 = 0;n2 < ((i32s) res_data.size()) - 2;n2++)
4037 {
4038 i32s r1[2] = { res_data[n2], res_data[n2 + 1] };
4039 i32s r2[2] = { res_data[n2 + 1], res_data[n2 + 2] };
4040
4041 i32s ind1 = readpdb_ReadData_sub1(atom_data, r1, "C", false);
4042 i32s ind2 = readpdb_ReadData_sub1(atom_data, r2, "N", false);
4043
4044 if (ind1 == NOT_DEFINED || ind2 == NOT_DEFINED)
4045 {
4046 assertion_failed(__FILE__, __LINE__, "res connection failed.");
4047 }
4048 else
4049 {
4050 bondtype bt = bondtype('S');
4051 bond newbond(atom_data[ind1].ref, atom_data[ind2].ref, bt);
4052 AddBond(newbond);
4053 }
4054 }
4055
4056 // check the c-terminal residue...
4057
4058 if (res_data.size() < 2) continue;
4059 i32s last[2] = { res_data[res_data.size() - 2], res_data[res_data.size() - 1] };
4060
4061 i32s ind1 = readpdb_ReadData_sub1(atom_data, last, "C", false);
4062 i32s ind2 = readpdb_ReadData_sub1(atom_data, last, "O", false);
4063
4064 if (ind1 == NOT_DEFINED || ind2 == NOT_DEFINED)
4065 {
4066 assertion_failed(__FILE__, __LINE__, "could not find c-term group.");
4067 }
4068 else
4069 {
4070 bondtype btD = bondtype('D');
4071 bond tmpbond(atom_data[ind1].ref, atom_data[ind2].ref, btD);
4072 iter_bl it1 = find(bond_list.begin(), bond_list.end(), tmpbond);
4073 if (it1 != bond_list.end()) (* it1).bt = btD; // should be OK already!
4074
4075 bondtype btS = bondtype('S');
4076 ind2 = readpdb_ReadData_sub1(atom_data, last, "OXT", false);
4077 if (ind2 != NOT_DEFINED)
4078 {
4079 bond newbond(atom_data[ind1].ref, atom_data[ind2].ref, btS);
4080 AddBond(newbond);
4081
4082 atom_data[ind2].ref->formal_charge = -1;
4083 }
4084 else
4085 {
4086 cout << _("missing terminal oxygen...") << endl;
4087
4088 i32s ind[3];
4089 ind[0] = readpdb_ReadData_sub1(atom_data, last, "CA", false);
4090 ind[1] = readpdb_ReadData_sub1(atom_data, last, "C", false);
4091 ind[2] = readpdb_ReadData_sub1(atom_data, last, "O", false);
4092
4093 const fGL * ref1; const fGL * ref2;
4094
4095 ref1 = atom_data[ind[1]].ref->GetCRD(0);
4096 ref2 = atom_data[ind[0]].ref->GetCRD(0);
4097 v3d<fGL> v1 = v3d<fGL>(ref1, ref2);
4098
4099 ref1 = atom_data[ind[1]].ref->GetCRD(0);
4100 ref2 = atom_data[ind[2]].ref->GetCRD(0);
4101 v3d<fGL> v2 = v3d<fGL>(ref1, ref2);
4102
4103 fGL tmp2 = v1.len(); fGL tmp3 = v1.spr(v2) / (tmp2 * tmp2);
4104 v1 = v1 * tmp3; v2 = v2 - v1;
4105
4106 fGL tmp4[3];
4107 const fGL * cdata = atom_data[ind[1]].ref->GetCRD(0);
4108 for (i32s n1 = 0;n1 < 3;n1++)
4109 {
4110 tmp4[n1] = cdata[n1];
4111 tmp4[n1] += v1.data[n1] - v2.data[n1];
4112 }
4113
4114 element el = element(8);
4115 atom newatom(el, tmp4, cs_vector.size());
4116
4117 newatom.formal_charge = -1;
4118
4119 AddAtom_lg(newatom);
4120 all_added_atoms.push_back(& atom_list.back());
4121
4122 bond newbond(atom_data[ind1].ref, (& atom_list.back()), btS);
4123 AddBond(newbond);
4124 }
4125 }
4126 }
4127
4128 // add disulphide bonds...
4129
4130 for (i32u n1 = 0;n1 < ssbond_data.size();n1++)
4131 {
4132 i32u counter = 0;
4133 while (counter < atom_data.size())
4134 {
4135 bool test1 = (ssbond_data[n1].chn_id == atom_data[counter].chn_id);
4136 bool test2 = (ssbond_data[n1].res_num == atom_data[counter].res_num);
4137 bool test3 = !strcmp(atom_data[counter].atm_name, "SG");
4138 if (test1 && test2 && test3) break; else counter++;
4139 }
4140
4141 if (counter == atom_data.size()) continue;
4142 ssbond_data[n1].ref = atom_data[counter].ref;
4143 }
4144
4145 for (i32u n1 = 0;n1 < ssbond_data.size();n1 += 2)
4146 {
4147 if (ssbond_data[n1 + 0].ref == NULL || ssbond_data[n1 + 1].ref == NULL)
4148 {
4149 cout << _("could not create bridge ");
4150 cout << ssbond_data[n1 + 0].chn_id << " " << ssbond_data[n1 + 0].res_num << " -> ";
4151 cout << ssbond_data[n1 + 1].chn_id << " " << ssbond_data[n1 + 1].res_num << endl;
4152 }
4153 else
4154 {
4155 bondtype bt = bondtype('S');
4156 bond newbond(ssbond_data[n1 + 0].ref, ssbond_data[n1 + 1].ref, bt);
4157 AddBond(newbond);
4158 }
4159 }
4160
4161 // add hydrogens for all added atoms...
4162
4163 for (i32u n1 = 0;n1 < all_added_atoms.size();n1++)
4164 {
4165 AddHydrogens(all_added_atoms[n1]);
4166 }
4167
4168 // ready...
4169
4170 cout << _("done.") << endl;
4171
4172 CenterCRDSet(0, true);
4173 }
4174
readpdb_ReadData_sub1(vector<readpdb_data_atom> & adata,i32s * range,const char * atom_name,bool flag)4175 i32s model::readpdb_ReadData_sub1(vector<readpdb_data_atom> & adata, i32s * range, const char * atom_name, bool flag)
4176 {
4177 for (i32s n1 = range[0];n1 < range[1];n1++)
4178 {
4179 if (!strcmp(adata[n1].atm_name, atom_name)) return n1;
4180 }
4181
4182 cout << _("atom ") << atom_name << _(" is missing...") << endl;
4183 return NOT_DEFINED;
4184 }
4185
readpdb_ReadData_sub2(vector<readpdb_data_atom> & adata,i32s * range,const char * at1,const char * at2,char btype)4186 void model::readpdb_ReadData_sub2(vector<readpdb_data_atom> & adata, i32s * range, const char * at1, const char * at2, char btype)
4187 {
4188 i32s ind1 = readpdb_ReadData_sub1(adata, range, at1, false);
4189 i32s ind2 = readpdb_ReadData_sub1(adata, range, at2, false);
4190
4191 if (ind1 == NOT_DEFINED || ind2 == NOT_DEFINED) return;
4192
4193 bondtype bt = bondtype(btype);
4194 bond newbond(adata[ind1].ref, adata[ind2].ref, bt);
4195 AddBond(newbond);
4196 }
4197
4198 // trajectory files-related stuff is here...
4199 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
4200 // todo : move all these trajectory-related things into a separate class? 20040610 TH
4201
OpenTrajectory(const char * fn)4202 void model::OpenTrajectory(const char * fn)
4203 {
4204 if (!trajfile)
4205 {
4206 traj_num_atoms = GetAtomCount();
4207 /* for (iter_al it1 = GetAtomsBegin();it1 != GetAtomsEnd();it1++)
4208 {
4209 if (!((* it1).flags & ATOMFLAG_IS_HIDDEN)) traj_num_atoms++;
4210 } */
4211
4212 trajfile = new ifstream(fn, ios::in | ios::binary);
4213 trajfile->seekg(8, ios::beg); // skip the file id...
4214
4215 int natoms;
4216 trajfile->read((char *) & natoms, sizeof(natoms));
4217
4218 if (natoms != traj_num_atoms)
4219 {
4220 ErrorMessage(_("The trajectory is incompatible with the current structure/setup!!!"));
4221 PrintToLog(_("incompatible file : different number of atoms!\n"));
4222 CloseTrajectory(); return;
4223 }
4224
4225 trajfile->read((char *) & total_traj_frames, sizeof(total_traj_frames));
4226
4227 stringstream str;
4228 str << _("the trajectory file contains ") << total_traj_frames << _(" frames.") << endl << ends;
4229
4230 PrintToLog(str.str().c_str());
4231
4232 current_traj_frame = 0;
4233 }
4234 else PrintToLog(_("trajectory file is already open!\n"));
4235 }
4236
ReadTrajectoryFrame(void)4237 void model::ReadTrajectoryFrame(void)
4238 {
4239 i32s place = 8 + 2 * sizeof(int); // skip the header...
4240 place += ((2 + 3) + 3 * traj_num_atoms) * sizeof(float) * current_traj_frame; // get the correct frame...
4241 place += 2 * sizeof(float); // skip epot and ekin...
4242
4243 trajfile->seekg(place, ios::beg);
4244
4245 float tmp;
4246
4247 float boundary[3];
4248 trajfile->read((char *) & tmp, sizeof(tmp)); boundary[0] = tmp;
4249 trajfile->read((char *) & tmp, sizeof(tmp)); boundary[1] = tmp;
4250 trajfile->read((char *) & tmp, sizeof(tmp)); boundary[2] = tmp;
4251
4252 if (boundary[0] >= 0.0)
4253 {
4254 saved_periodic_box_HALFdim[0] = boundary[0];
4255 saved_periodic_box_HALFdim[1] = boundary[1];
4256 saved_periodic_box_HALFdim[2] = boundary[2];
4257 }
4258 else if (boundary[1] >= 0.0)
4259 {
4260 saved_boundary_potential_rad_solute = boundary[1];
4261 saved_boundary_potential_rad_solvent = boundary[2];
4262 }
4263
4264 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
4265 {
4266 // if ((* it1).flags & ATOMFLAG_IS_HIDDEN) continue; // currently all coordinates are written...
4267
4268 fGL cdata[3];
4269 for (i32s t4 = 0;t4 < 3;t4++)
4270 {
4271 trajfile->read((char *) & tmp, sizeof(tmp));
4272 cdata[t4] = tmp;
4273 }
4274
4275 (* it1).SetCRD(0, cdata[0], cdata[1], cdata[2]);
4276 }
4277 }
4278
CloseTrajectory(void)4279 void model::CloseTrajectory(void)
4280 {
4281 if (trajfile != NULL)
4282 {
4283 trajfile->close();
4284 delete trajfile;
4285
4286 trajfile = NULL;
4287 }
4288 }
4289
GetTotalFrames(void)4290 i32s model::GetTotalFrames(void)
4291 {
4292 return total_traj_frames;
4293 }
4294
GetTrajectoryFile(void)4295 ifstream * model::GetTrajectoryFile(void)
4296 {
4297 return trajfile;
4298 }
4299
GetCurrentFrame(void)4300 i32s model::GetCurrentFrame(void)
4301 {
4302 return current_traj_frame;
4303 }
4304
SetCurrentFrame(i32s p1)4305 void model::SetCurrentFrame(i32s p1)
4306 {
4307 current_traj_frame = p1;
4308 }
4309
WriteTrajectoryHeader(ofstream & ofile,int total_frames)4310 void model::WriteTrajectoryHeader(ofstream & ofile, int total_frames)
4311 {
4312 const char file_id[10] = "traj_v11";
4313 const int number_of_atoms = GetAtomCount();
4314
4315 ofile.write((char *) file_id, 8); // file id, 8 chars.
4316 ofile.write((char *) & number_of_atoms, sizeof(number_of_atoms)); // number of atoms, int.
4317 ofile.write((char *) & total_frames, sizeof(total_frames)); // total number of frames, int.
4318 }
4319
WriteTrajectoryFrame(ofstream & ofile,moldyn * dyn)4320 void model::WriteTrajectoryFrame(ofstream & ofile, moldyn * dyn)
4321 {
4322 const float ekin = dyn->GetEKin();
4323 const float epot = dyn->GetEPot();
4324
4325 float boundary[3] = { -1.0, -1.0, -1.0 };
4326
4327 engine_bp * eng_bp = dynamic_cast<engine_bp *>(dyn->eng);
4328 if (eng_bp != NULL)
4329 {
4330 // the 1st value is not modified...
4331 boundary[1] = eng_bp->bp_rad_solute;
4332 boundary[2] = eng_bp->bp_rad_solvent;
4333 }
4334
4335 engine_pbc * eng_pbc = dynamic_cast<engine_pbc *>(dyn->eng);
4336 if (eng_pbc != NULL)
4337 {
4338 boundary[0] = eng_pbc->box_HALFdim[0];
4339 boundary[1] = eng_pbc->box_HALFdim[1];
4340 boundary[2] = eng_pbc->box_HALFdim[2];
4341 }
4342
4343 ofile.write((char *) & ekin, sizeof(ekin)); // kinetic energy, float.
4344 ofile.write((char *) & epot, sizeof(epot)); // potential energy, float.
4345
4346 ofile.write((char *) boundary, sizeof(float) * 3);
4347
4348 for (iter_al itx = GetAtomsBegin();itx != GetAtomsEnd();itx++)
4349 {
4350 const fGL * cdata = (* itx).GetCRD(0);
4351 for (i32s t4 = 0;t4 < 3;t4++) // all coordinates, float.
4352 {
4353 float t1a = cdata[t4];
4354 ofile.write((char *) & t1a, sizeof(t1a));
4355 }
4356 }
4357 }
4358
EvaluateBFact(void)4359 void model::EvaluateBFact(void)
4360 {
4361 if (!trajfile)
4362 {
4363 PrintToLog(_("EvaluateBFact() : trajectory file not open!\n"));
4364 return;
4365 }
4366
4367 vector<atom *> av; // pick all selected atoms here...
4368 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
4369 {
4370 if ((* it1).flags & ATOMFLAG_USER_SELECTED) av.push_back(& (* it1));
4371 }
4372
4373 if (!av.size())
4374 {
4375 PrintToLog(_("EvaluateBFact() : no selected atoms!\n"));
4376 return;
4377 }
4378
4379 fGL * avrg_str = new fGL[av.size() * 3];
4380 fGL * b_tab = new fGL[av.size()];
4381
4382 for (i32u ac = 0;ac < av.size();ac++)
4383 {
4384 avrg_str[ac * 3 + 0] = 0.0;
4385 avrg_str[ac * 3 + 1] = 0.0;
4386 avrg_str[ac * 3 + 2] = 0.0;
4387
4388 b_tab[ac] = 0.0;
4389 }
4390
4391 for (i32s n1 = 0;n1 < GetTotalFrames();n1++)
4392 {
4393 SetCurrentFrame(n1);
4394 ReadTrajectoryFrame();
4395
4396 for (i32u ac = 0;ac < av.size();ac++)
4397 {
4398 const fGL * tmpc = av[ac]->GetCRD(0);
4399
4400 avrg_str[ac * 3 + 0] += tmpc[0];
4401 avrg_str[ac * 3 + 1] += tmpc[1];
4402 avrg_str[ac * 3 + 2] += tmpc[2];
4403 }
4404 }
4405
4406 for (i32u ac = 0;ac < av.size();ac++)
4407 {
4408 avrg_str[ac * 3 + 0] /= (fGL) GetTotalFrames();
4409 avrg_str[ac * 3 + 1] /= (fGL) GetTotalFrames();
4410 avrg_str[ac * 3 + 2] /= (fGL) GetTotalFrames();
4411 }
4412
4413 for (i32s n1 = 0;n1 < GetTotalFrames();n1++)
4414 {
4415 SetCurrentFrame(n1);
4416 ReadTrajectoryFrame();
4417
4418 for (i32u ac = 0;ac < av.size();ac++)
4419 {
4420 const fGL * tmpc = av[ac]->GetCRD(0);
4421
4422 fGL dx = avrg_str[ac * 3 + 0] - tmpc[0];
4423 fGL dy = avrg_str[ac * 3 + 1] - tmpc[1];
4424 fGL dz = avrg_str[ac * 3 + 2] - tmpc[2];
4425
4426 b_tab[ac] += dx * dx;
4427 b_tab[ac] += dy * dy;
4428 b_tab[ac] += dz * dz;
4429 }
4430 }
4431
4432 for (i32u ac = 0;ac < av.size();ac++)
4433 {
4434 b_tab[ac] /= (fGL) GetTotalFrames();
4435 }
4436
4437 for (i32u ac = 0;ac < av.size();ac++)
4438 {
4439 ostringstream txts;
4440 txts << "atom " << av[ac]->index << " ";
4441 txts << "displacement " << b_tab[ac] << " nm^2 = " << (b_tab[ac] * 100.0) << " �^2 ";
4442 txts << " -> Bfact " << 78.957 * (b_tab[ac] * 100.0) << endl << ends;
4443
4444 PrintToLog(txts.str().c_str());// 20060223 assertion g_utf8_validate() failed???
4445 cout << txts.str().c_str();
4446 }
4447
4448 delete[] avrg_str;
4449 delete[] b_tab;
4450 }
4451
EvaluateDiffConst(double dt)4452 void model::EvaluateDiffConst(double dt)
4453 {
4454 if (!trajfile)
4455 {
4456 PrintToLog(_("EvaluateDiffConst() : trajectory file not open!\n"));
4457 return;
4458 }
4459
4460 vector<atom *> av; // pick all selected atoms here...
4461 for (iter_al it1 = atom_list.begin();it1 != atom_list.end();it1++)
4462 {
4463 if ((* it1).flags & ATOMFLAG_USER_SELECTED) av.push_back(& (* it1));
4464 }
4465
4466 if (!av.size())
4467 {
4468 PrintToLog(_("EvaluateDiffConst() : no selected atoms!\n"));
4469 return;
4470 }
4471
4472 fGL * init_loc = new fGL[av.size() * 3];
4473 double * dc_tab = new double[av.size()];
4474
4475 SetCurrentFrame(0);
4476 ReadTrajectoryFrame();
4477
4478 for (i32u ac = 0;ac < av.size();ac++)
4479 {
4480 const fGL * tmpc = av[ac]->GetCRD(0);
4481
4482 init_loc[ac * 3 + 0] = tmpc[0];
4483 init_loc[ac * 3 + 1] = tmpc[1];
4484 init_loc[ac * 3 + 2] = tmpc[2];
4485
4486 dc_tab[ac] = 0.0;
4487 }
4488
4489 double time = 0.0;
4490 for (i32s n1 = 1;n1 < GetTotalFrames();n1++)
4491 {
4492 time += dt;
4493
4494 SetCurrentFrame(n1);
4495 ReadTrajectoryFrame();
4496
4497 for (i32u ac = 0;ac < av.size();ac++)
4498 {
4499 const fGL * tmpc = av[ac]->GetCRD(0);
4500
4501 double dist = 0.0; double tmpd;
4502 tmpd = tmpc[0] - init_loc[ac * 3 + 0]; dist += tmpd * tmpd;
4503 tmpd = tmpc[1] - init_loc[ac * 3 + 1]; dist += tmpd * tmpd;
4504 tmpd = tmpc[2] - init_loc[ac * 3 + 2]; dist += tmpd * tmpd;
4505
4506 // convert nm^2 into cm^2 ; 10^-14
4507 // convert fs into s ; 10^-15
4508
4509 double dc = (dist * 1.0e-14) / (time * 1.0e-15);
4510 dc_tab[ac] += dc;
4511 }
4512 }
4513
4514 for (i32u ac = 0;ac < av.size();ac++)
4515 {
4516 dc_tab[ac] /= (double) (GetTotalFrames() - 1);
4517 }
4518
4519 for (i32u ac = 0;ac < av.size();ac++)
4520 {
4521 ostringstream txts;
4522 txts << "atom " << av[ac]->index << " ";
4523 txts << "diffconst " << (dc_tab[ac] * 1.0e+5) << " * 10^-5 cm^2/s" << endl << ends;
4524
4525 PrintToLog(txts.str().c_str());// 20060223 assertion g_utf8_validate() failed???
4526 cout << txts.str().c_str();
4527 }
4528
4529 delete[] init_loc;
4530 delete[] dc_tab;
4531 }
4532
4533 // the ecomp stuff is here...
4534 // ^^^^^^^^^^^^^^^^^^^^^^^^^^
4535
ecomp_AddGroup(const char * gn)4536 int model::ecomp_AddGroup(const char * gn)
4537 {
4538 if (gn == NULL) return NOT_DEFINED;
4539
4540 const int new_grp_i = (int) ecomp_grp_names.size();
4541
4542 char * copy_of_grp_name = new char[strlen(gn) + 1];
4543 strcpy(copy_of_grp_name, gn);
4544
4545 ecomp_grp_names.push_back(copy_of_grp_name);
4546
4547 return new_grp_i;
4548 }
4549
ecomp_DeleteGroups(void)4550 void model::ecomp_DeleteGroups(void)
4551 {
4552 // here we will delete all groups added by the user, but will not
4553 // delete the "default" group with an index zero.
4554
4555 // it would be just confusing to allow deletion of some groups,
4556 // since it would either change the indices or create gaps...
4557
4558 while (ecomp_grp_names.size() > 1)
4559 {
4560 delete[] ecomp_grp_names.back();
4561 ecomp_grp_names.pop_back();
4562 }
4563 }
4564
4565 /*################################################################################################*/
4566
crd_set(void)4567 crd_set::crd_set(void)
4568 {
4569 description = NULL;
4570
4571 accum_weight = 1.0; // accum_value is not updated here?!?!?!?!?!
4572 visible = false;
4573 }
4574
crd_set(const crd_set & p1)4575 crd_set::crd_set(const crd_set & p1)
4576 {
4577 description = NULL;
4578 SetDescription(p1.description);
4579
4580 accum_weight = p1.accum_weight;
4581 accum_value = p1.accum_value;
4582 visible = p1.visible;
4583 }
4584
~crd_set(void)4585 crd_set::~crd_set(void)
4586 {
4587 if (description != NULL) delete[] description;
4588 }
4589
SetDescription(const char * p1)4590 void crd_set::SetDescription(const char * p1)
4591 {
4592 if (description != NULL) delete[] description;
4593
4594 if (p1 != NULL)
4595 {
4596 description = new char[strlen(p1) + 1];
4597 strcpy(description, p1);
4598 }
4599 else description = NULL;
4600 }
4601
4602 /*################################################################################################*/
4603
4604 #include "libghemicalconfig.h"
4605
libghemical_init(const char * tmp_path)4606 void libghemical_init(const char * tmp_path)
4607 {
4608 static singleton_cleaner<sequencebuilder> amino_cleaner;
4609 static singleton_cleaner<sequencebuilder> nucleic_cleaner;
4610
4611 static singleton_cleaner<default_tables> deftab_cleaner;
4612 static singleton_cleaner<tripos52_tables> tripos52_cleaner;
4613
4614 static int counter = 0;
4615 if (counter++ != 0)
4616 {
4617 assertion_failed(__FILE__, __LINE__, "libghemical_init() was called more than once!");
4618 }
4619
4620 #ifdef ENABLE_NLS
4621 bindtextdomain(GETTEXT_PACKAGE, LOCALE_DIR);
4622 bind_textdomain_codeset(GETTEXT_PACKAGE, "UTF-8");
4623 #endif // ENABLE_NLS
4624
4625 // in windows we need to get the installation path at runtime ; the old (pre-2008-july) solution was
4626 // to use a registry key but with gtk+/glib it's better to ask it directly from the application...
4627
4628 strcpy(model::libdata_path, tmp_path);
4629
4630 // now when model::libdata_path is defined, we can start using it...
4631 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
4632
4633 model::amino_builder = new sequencebuilder(chn_info::amino_acid, AMINO_BUILDER_FILE);
4634 amino_cleaner.SetInstance(model::amino_builder);
4635
4636 model::nucleic_builder = new sequencebuilder(chn_info::nucleic_acid, NUCLEIC_BUILDER_FILE);
4637 nucleic_cleaner.SetInstance(model::nucleic_builder);
4638
4639 default_tables::instance = default_tables::GetInstance();
4640 deftab_cleaner.SetInstance(default_tables::instance);
4641
4642 tripos52_tables::instance = tripos52_tables::GetInstance();
4643 tripos52_cleaner.SetInstance(tripos52_tables::instance);
4644
4645 // ready...
4646 // ^^^^^^^^
4647 }
4648
4649 /*################################################################################################*/
4650
4651 // eof
4652