1 // TAB_MM_TRIPOS52.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 "tab_mm_tripos52.h"
23
24 #include "typerule.h"
25 #include "utility.h"
26
27 #include "local_i18n.h"
28
29 #include <cstring>
30 #include <fstream>
31 #include <sstream>
32 #include <iomanip>
33 using namespace std;
34
35 /*################################################################################################*/
36
37 tripos52_tables * tripos52_tables::instance = NULL;
38 //singleton_cleaner<tripos52_tables> tripos52_cleaner(tripos52_tables::GetInstance());
39
tripos52_tables(void)40 tripos52_tables::tripos52_tables(void)
41 {
42 ifstream file;
43 file.unsetf(ios::dec | ios::oct | ios::hex);
44
45 const char * fn;
46 char buffer[1024];
47
48 ostream * ostr = NULL; // do not print output.
49 // ostream * ostr = & cout; // print output to cout.
50
51 /*##############################################*/
52 /*##############################################*/
53
54 fn = "param_mm/tripos52/atomtypes.txt";
55 model::OpenLibDataFile(file, false, fn);
56
57 if (ostr != NULL) (* ostr) << _("reading file \"") << fn << "\": ";
58
59 while (file.peek() != '#') // #end
60 {
61 if (file.peek() == '0') // 0x????
62 {
63 tripos52_at newat;
64 file >> newat.atomtype;
65
66 while (file.peek() != '(') file.get();
67 newat.tr = new typerule(& file, ostr);
68
69 while (file.get() != '\"');
70 file.getline(buffer, sizeof(buffer), '\"');
71 newat.description = new char[strlen(buffer) + 1];
72 strcpy(newat.description, buffer);
73
74 at_vector.push_back(newat);
75 }
76
77 file.getline(buffer, sizeof(buffer));
78 }
79
80 file.close();
81
82 if (ostr != NULL) (* ostr) << _("found ") << at_vector.size() << _(" atomtypes.") << endl;
83
84 /*##############################################*/
85 /*##############################################*/
86
87 fn = "param_mm/tripos52/parameters1.txt";
88 model::OpenLibDataFile(file, false, fn);
89
90 if (ostr != NULL) (* ostr) << _("reading file \"") << fn << "\": ";
91
92 while (file.peek() != '#') // #end
93 {
94 if (file.peek() == '0') // 0x????
95 {
96 tripos52_bs tmp; char bt[16];
97 file >> tmp.atmtp[0] >> tmp.atmtp[1] >> bt;
98 file >> tmp.param[0] >> tmp.param[1];
99
100 tmp.bndtp = bondtype(bt[0]);
101 bs_vector.push_back(tmp);
102 }
103
104 file.getline(buffer, sizeof(buffer));
105 }
106
107 file.close();
108
109 if (ostr != NULL) (* ostr) << _("found ") << bs_vector.size() << _(" bs-terms.") << endl;
110
111 /*##############################################*/
112 /*##############################################*/
113
114 fn = "param_mm/tripos52/parameters2.txt";
115 model::OpenLibDataFile(file, false, fn);
116
117 if (ostr != NULL) (* ostr) << _("reading file \"") << fn << "\": ";
118
119 while (file.peek() != '#') // #end
120 {
121 if (file.peek() == '0') // 0x????
122 {
123 tripos52_ab tmp; char bt[16];
124 file >> tmp.atmtp[0] >> tmp.atmtp[1] >> tmp.atmtp[2] >> bt;
125 file >> tmp.param[0] >> tmp.param[1];
126
127 for (i32s n1 = 0;n1 < 2;n1++) tmp.bndtp[n1] = bondtype(bt[n1]);
128 ab_vector.push_back(tmp);
129 }
130
131 file.getline(buffer, sizeof(buffer));
132 }
133
134 file.close();
135
136 if (ostr != NULL) (* ostr) << _("found ") << ab_vector.size() << _(" ab-terms.") << endl;
137
138 /*##############################################*/
139 /*##############################################*/
140
141 fn = "param_mm/tripos52/parameters3.txt";
142 model::OpenLibDataFile(file, false, fn);
143
144 if (ostr != NULL) (* ostr) << _("reading file \"") << fn << "\": ";
145
146 while (file.peek() != '#') // #end
147 {
148 if (file.peek() == '0') // 0x????
149 {
150 tripos52_tr tmp; char bt[16];
151 file >> tmp.atmtp[0] >> tmp.atmtp[1] >> tmp.atmtp[2] >> tmp.atmtp[3] >> bt;
152 file >> tmp.k >> tmp.s;
153
154 for (i32s n1 = 0;n1 < 3;n1++) tmp.bndtp[n1] = bondtype(bt[n1]);
155 tr_vector.push_back(tmp);
156 }
157
158 file.getline(buffer, sizeof(buffer));
159 }
160
161 file.close();
162
163 if (ostr != NULL) (* ostr) << _("found ") << tr_vector.size() << _(" tr-terms.") << endl;
164
165 /*##############################################*/
166 /*##############################################*/
167
168 fn = "param_mm/tripos52/parameters4.txt";
169 model::OpenLibDataFile(file, false, fn);
170
171 if (ostr != NULL) (* ostr) << _("reading file \"") << fn << "\": ";
172
173 while (file.peek() != '#') // #end
174 {
175 if (file.peek() == '0') // 0x????
176 {
177 tripos52_lj tmp;
178 file >> tmp.type;
179 file >> tmp.r >> tmp.k;
180
181 lj_vector.push_back(tmp);
182 }
183
184 file.getline(buffer, sizeof(buffer));
185 }
186
187 file.close();
188
189 if (ostr != NULL) (* ostr) << _("found ") << lj_vector.size() << _(" lj-datasets.") << endl;
190
191 /*##############################################*/
192 /*##############################################*/
193
194 fn = "param_mm/tripos52/parameters5.txt";
195 model::OpenLibDataFile(file, false, fn);
196
197 if (ostr != NULL) (* ostr) << _("reading file \"") << fn << "\": ";
198
199 while (file.peek() != '#') // #end
200 {
201 if (file.peek() == '0') // 0x????
202 {
203 tripos52_ci tmp; char bt[16];
204 file >> tmp.atmtp[0] >> tmp.atmtp[1] >> bt;
205 file >> tmp.delta;
206
207 tmp.bndtp = bondtype(bt[0]);
208 ci_vector.push_back(tmp);
209 }
210
211 file.getline(buffer, sizeof(buffer));
212 }
213
214 file.close();
215
216 if (ostr != NULL) (* ostr) << _("found ") << ci_vector.size() << _(" ci-datasets.") << endl;
217 }
218
~tripos52_tables(void)219 tripos52_tables::~tripos52_tables(void)
220 {
221 for (i32u n1 = 0;n1 < at_vector.size();n1++)
222 {
223 delete at_vector[n1].tr;
224 delete[] at_vector[n1].description;
225 }
226 }
227
GetInstance(void)228 tripos52_tables * tripos52_tables::GetInstance(void)
229 {
230 if (instance != NULL) return instance;
231 else return (instance = new tripos52_tables());
232 }
233
PrintAllTypeRules(ostream & p1)234 void tripos52_tables::PrintAllTypeRules(ostream & p1)
235 {
236 for (i32u n1 = 0;n1 < at_vector.size();n1++)
237 {
238 p1 << (n1 + 1) << ": 0x" << hex << setw(4) << setfill('0') << at_vector[n1].atomtype << dec;
239 p1 << " (" << (* at_vector[n1].tr) << ") \"" << at_vector[n1].description << "\"" << endl;
240 }
241
242 p1 << at_vector.size() << _(" entries.") << endl;
243 }
244
UpdateTypes(setup * su)245 i32u tripos52_tables::UpdateTypes(setup * su)
246 {
247 model * mdl = su->GetModel();
248
249 if (mdl->verbosity >= 3)
250 {
251 ostringstream str;
252 str << _("Setting up atom types and formal charges...") << endl << ends;
253
254 mdl->PrintToLog(str.str().c_str());
255 }
256
257 i32u errors = 0;
258
259 for (iter_al it1 = mdl->GetAtomsBegin();it1 != mdl->GetAtomsEnd();it1++)
260 {
261 i32u at_range[2];
262
263 at_range[0] = 0;
264 while (true)
265 {
266 if (at_range[0] == at_vector.size()) break;
267 if ((at_vector[at_range[0]].atomtype >> 8) == (* it1).el.GetAtomicNumber()) break;
268
269 at_range[0]++;
270 }
271
272 at_range[1] = at_range[0];
273 while (true)
274 {
275 if (at_range[1] == at_vector.size()) break;
276 if ((at_vector[at_range[1]].atomtype >> 8) != (* it1).el.GetAtomicNumber()) break;
277
278 at_range[1]++;
279 }
280
281 i32s index = NOT_DEFINED;
282 for (i32u n1 = at_range[0];n1 < at_range[1];n1++)
283 {
284 bool flag = at_vector[n1].tr->Check(mdl, & (* it1), 0);
285 if (flag) index = n1;
286 }
287
288 if (index != NOT_DEFINED)
289 {
290 (* it1).atmtp = at_vector[index].atomtype;
291
292 // no formal charges defined in tripos5.2???
293 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
294
295 (* it1).charge = 0.0;
296 }
297 else
298 {
299 if (mdl->verbosity >= 2)
300 {
301 ostringstream str;
302 str << _("WARNING : could not determine atomtype (atom index = ") << (* it1).index << ")." << endl << ends;
303
304 mdl->PrintToLog(str.str().c_str());
305 }
306
307 (* it1).atmtp = NOT_DEFINED;
308 (* it1).charge = 0.0;
309
310 (* it1).flags |= ATOMFLAG_USER_SELECTED;
311 errors++;
312 }
313 }
314
315 return errors;
316 }
317
UpdateCharges(setup * su)318 i32u tripos52_tables::UpdateCharges(setup * su)
319 {
320 model * mdl = su->GetModel();
321
322 if (mdl->verbosity >= 3)
323 {
324 ostringstream str;
325 str << _("Setting up partial charges...") << endl << ends;
326 mdl->PrintToLog(str.str().c_str());
327 }
328
329 i32u errors = 0;
330
331 for (iter_bl it1 = mdl->GetBondsBegin();it1 != mdl->GetBondsEnd();it1++)
332 {
333 f64 delta = tripos52_tables::GetInstance()->GetChargeInc(& (* it1), mdl);
334
335 (* it1).atmr[0]->charge -= delta;
336 (* it1).atmr[1]->charge += delta;
337 }
338
339 return errors;
340 }
341
GetChargeInc(bond * ref,model * mdl)342 f64 tripos52_tables::GetChargeInc(bond * ref, model * mdl)
343 {
344 i32s atmtp[2];
345 for (i32s n1 = 0;n1 < 2;n1++)
346 {
347 atmtp[n1] = ref->atmr[n1]->atmtp;
348 }
349
350 for (i32u n1 = 0;n1 < ci_vector.size();n1++)
351 {
352 if (ci_vector[n1].bndtp.GetValue() != ref->bt.GetValue()) continue;
353
354 for (i32s dir = 0;dir < 2;dir++)
355 {
356 bool test1 = (ci_vector[n1].atmtp[0] == atmtp[dir]);
357 bool test2 = (ci_vector[n1].atmtp[1] == atmtp[!dir]);
358
359 if (test1 && test2)
360 {
361 if (!dir) return +ci_vector[n1].delta;
362 else return -ci_vector[n1].delta;
363 }
364 }
365 }
366
367 if (mdl != NULL && mdl->verbosity >= 2)
368 {
369 ostringstream str;
370 str << _("WARNING : there was no record for the following ci: ");
371 str << "0x" << hex << setw(4) << setfill('0') << atmtp[0] << dec << " ";
372 str << "0x" << hex << setw(4) << setfill('0') << atmtp[1] << dec << " ";
373 str << ref->bt.GetValue() << endl << ends;
374
375 mdl->PrintToLog(str.str().c_str());
376 }
377
378 return 0.0;
379 }
380
Init(eng1_mm * eng,mm_tripos52_bt1 * ref,i32s bt)381 bool tripos52_tables::Init(eng1_mm * eng, mm_tripos52_bt1 * ref, i32s bt)
382 {
383 atom ** atmtab = eng->GetSetup()->GetMMAtoms();
384
385 i32s atmtp[2];
386 for (i32s n1 = 0;n1 < 2;n1++)
387 {
388 atmtp[n1] = atmtab[ref->atmi[n1]]->atmtp;
389 }
390
391 for (i32u n1 = 0;n1 < bs_vector.size();n1++)
392 {
393 if (bs_vector[n1].bndtp.GetValue() != bt) continue;
394
395 bool flag = false;
396 for (i32s dir = 0;dir < 2;dir++)
397 {
398 bool test1 = (bs_vector[n1].atmtp[0] == atmtp[dir]);
399 bool test2 = (bs_vector[n1].atmtp[1] == atmtp[!dir]);
400
401 if (test1 && test2) flag = true;
402
403 bool wc1 = (bs_vector[n1].atmtp[0] == 0xFFFF);
404 bool wc2 = (bs_vector[n1].atmtp[1] == 0xFFFF);
405
406 if (wc1 && test2) flag = true;
407 if (test1 && wc2) flag = true;
408 if (wc1 && wc2) flag = true;
409
410 if (flag) break;
411 }
412
413 // we will convert units to [nm] and [kJ/mol]...
414
415 if (flag)
416 {
417 ref->opt = 0.1 * bs_vector[n1].param[0];
418 ref->fc = 418.68 * bs_vector[n1].param[1];
419
420 return true;
421 }
422 }
423
424 model * mdl = eng->GetSetup()->GetModel();
425 if (mdl->verbosity >= 2)
426 {
427 ostringstream str;
428 str << _("WARNING : unknown bst: ");
429 str << "0x" << hex << setw(4) << setfill('0') << atmtp[0] << dec << " ";
430 str << "0x" << hex << setw(4) << setfill('0') << atmtp[1] << dec << " ";
431 str << bt << endl << ends;
432
433 mdl->PrintToLog(str.str().c_str());
434 }
435
436 ref->opt = 0.1 * 1.100;
437 ref->fc = 418.68 * 500.0;
438
439 return false;
440 }
441
Init(eng1_mm * eng,mm_tripos52_bt2 * ref,i32s * bt)442 bool tripos52_tables::Init(eng1_mm * eng, mm_tripos52_bt2 * ref, i32s * bt)
443 {
444 atom ** atmtab = eng->GetSetup()->GetMMAtoms();
445
446 i32s atmtp[3];
447 for (i32s n1 = 0;n1 < 3;n1++)
448 {
449 atmtp[n1] = atmtab[ref->atmi[n1]]->atmtp;
450 }
451
452 for (i32u n1 = 0;n1 < ab_vector.size();n1++)
453 {
454 if (ab_vector[n1].atmtp[1] != atmtp[1]) continue;
455
456 // bondtype checking not yet implemented....
457 // it will be basically similar to the above, but must be moved into the dir-loop!!!
458 // it will be basically similar to the above, but must be moved into the dir-loop!!!
459 // it will be basically similar to the above, but must be moved into the dir-loop!!!
460
461 bool flag = false;
462 for (i32s dir = 0;dir < 2;dir++)
463 {
464 i32s index[2];
465 index[0] = (!dir ? 0 : 2);
466 index[1] = (!dir ? 2 : 0);
467
468 bool test1 = (ab_vector[n1].atmtp[0] == atmtp[index[0]]);
469 bool test2 = (ab_vector[n1].atmtp[2] == atmtp[index[1]]);
470
471 if (test1 && test2) flag = true;
472
473 bool wc1 = (ab_vector[n1].atmtp[0] == 0xFFFF);
474 bool wc2 = (ab_vector[n1].atmtp[2] == 0xFFFF);
475
476 if (wc1 && test2) flag = true;
477 if (test1 && wc2) flag = true;
478 if (wc1 && wc2) flag = true;
479
480 if (flag) break;
481 }
482
483 // we will convert units to [rad] and [kJ/mol]...
484
485 if (flag)
486 {
487 ref->opt = M_PI * ab_vector[n1].param[0] / 180.0;
488 ref->fc = 13744.5 * ab_vector[n1].param[1];
489
490 return true;
491 }
492 }
493
494 model * mdl = eng->GetSetup()->GetModel();
495 if (mdl->verbosity >= 2)
496 {
497 ostringstream str;
498 str << _("WARNING : unknown abn: ") << hex;
499 str << "0x" << hex << setw(4) << setfill('0') << atmtp[0] << dec << " ";
500 str << "0x" << hex << setw(4) << setfill('0') << atmtp[1] << dec << " ";
501 str << "0x" << hex << setw(4) << setfill('0') << atmtp[2] << dec << " ";
502 str << bt[0] << " " << bt[1] << endl << ends;
503
504 mdl->PrintToLog(str.str().c_str());
505 }
506
507 ref->opt = M_PI * 120.0 / 180.0;
508 ref->fc = 13744.5 * 0.020;
509
510 return false;
511 }
512
Init(eng1_mm * eng,mm_tripos52_bt3 * ref,i32s * bt)513 bool tripos52_tables::Init(eng1_mm * eng, mm_tripos52_bt3 * ref, i32s * bt)
514 {
515 atom ** atmtab = eng->GetSetup()->GetMMAtoms();
516
517 i32s atmtp[4];
518 for (i32s n1 = 0;n1 < 4;n1++)
519 {
520 atmtp[n1] = atmtab[ref->atmi[n1]]->atmtp;
521 }
522
523 for (i32u n1 = 0;n1 < tr_vector.size();n1++)
524 {
525 if (tr_vector[n1].bndtp[1].GetValue() != bt[1]) continue;
526
527 // proper bondtype checking not yet implemented....
528 // it will be basically similar to the above, but must be moved into the dir-loop!!!
529 // it will be basically similar to the above, but must be moved into the dir-loop!!!
530 // it will be basically similar to the above, but must be moved into the dir-loop!!!
531
532 bool flag = false;
533 for (i32s dir = 0;dir < 2;dir++)
534 {
535 i32s index[4];
536 index[0] = (!dir ? 0 : 3);
537 index[1] = (!dir ? 1 : 2);
538 index[2] = (!dir ? 2 : 1);
539 index[3] = (!dir ? 3 : 0);
540
541 bool test1 = (tr_vector[n1].atmtp[0] == atmtp[index[0]]);
542 bool test2 = (tr_vector[n1].atmtp[1] == atmtp[index[1]]);
543 bool test3 = (tr_vector[n1].atmtp[2] == atmtp[index[2]]);
544 bool test4 = (tr_vector[n1].atmtp[3] == atmtp[index[3]]);
545
546 if (test1 && test2 && test3 && test4) flag = true;
547
548 bool wc1 = (tr_vector[n1].atmtp[0] == 0xFFFF);
549 bool wc2 = (tr_vector[n1].atmtp[3] == 0xFFFF);
550
551 if (wc1 && test2 && test3 && test4) flag = true;
552 if (test1 && test2 && test3 && wc2) flag = true;
553 if (wc1 && test2 && test3 && wc2) flag = true;
554
555 if (flag) break;
556 }
557
558 // we will convert units to [kJ/mol]...
559
560 if (flag)
561 {
562 ref->k = 4.1868 * tr_vector[n1].k;
563 ref->s = tr_vector[n1].s;
564
565 return true;
566 }
567 }
568
569 model * mdl = eng->GetSetup()->GetModel();
570 if (mdl->verbosity >= 2)
571 {
572 ostringstream str;
573 str << _("WARNING : unknown tor: ") << hex;
574 str << "0x" << hex << setw(4) << setfill('0') << atmtp[0] << dec << " ";
575 str << "0x" << hex << setw(4) << setfill('0') << atmtp[1] << dec << " ";
576 str << "0x" << hex << setw(4) << setfill('0') << atmtp[2] << dec << " ";
577 str << "0x" << hex << setw(4) << setfill('0') << atmtp[3] << dec << " ";
578 str << bt[0] << " " << bt[1] << " " << bt[2] << endl << ends;
579
580 mdl->PrintToLog(str.str().c_str());
581 }
582
583 ref->k = 0.0;
584 ref->s = +1.0;
585
586 return false;
587 }
588
Init(eng1_mm * eng,mm_tripos52_nbt1 * ref,bool is14)589 bool tripos52_tables::Init(eng1_mm * eng, mm_tripos52_nbt1 * ref, bool is14)
590 {
591 atom ** atmtab = eng->GetSetup()->GetMMAtoms();
592
593 i32s atmtp[2];
594 for (i32s n1 = 0;n1 < 2;n1++)
595 {
596 atmtp[n1] = atmtab[ref->atmi[n1]]->atmtp;
597 }
598
599 i32u index[2];
600 for (i32s n1 = 0;n1 < 2;n1++)
601 {
602 index[n1] = 0;
603
604 while (index[n1] < lj_vector.size())
605 {
606 if (lj_vector[index[n1]].type == atmtp[n1]) break;
607 else index[n1]++;
608 }
609
610 if (index[n1] == lj_vector.size())
611 {
612 model * mdl = eng->GetSetup()->GetModel();
613 if (mdl->verbosity >= 2)
614 {
615 ostringstream str;
616 str << _("WARNING : bad atomtype ; using hydrogen instead...") << endl << ends;
617
618 mdl->PrintToLog(str.str().c_str());
619 }
620
621 index[n1] = 0;
622 }
623 }
624
625 // we will convert units to [nm] and [kJ/mol]...
626
627 f64 energy = 4.1868 * sqrt(lj_vector[index[0]].k * lj_vector[index[1]].k);
628 f64 optdist = 0.1 * (lj_vector[index[0]].r + lj_vector[index[1]].r);
629
630 f64 charge1 = atmtab[ref->atmi[0]]->charge;
631 f64 charge2 = atmtab[ref->atmi[1]]->charge;
632 ref->qq = 138.9354518 * charge1 * charge2;
633
634 if (is14)
635 {
636 energy *= 0.5;
637 ref->qq *= 0.5;
638 }
639
640 ref->k1 = optdist * pow(1.0 * energy, 1.0 / 12.0);
641 ref->k2 = optdist * pow(2.0 * energy, 1.0 / 6.0);
642
643 return true;
644 }
645
646 /*################################################################################################*/
647
648 // eof
649