1 // ENGINE.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 "engine.h"
23
24 #include "v3d.h"
25 #include "model.h"
26
27 #include "local_i18n.h"
28 #include "notice.h"
29
30 #include <vector>
31 #include <sstream>
32 using namespace std;
33
34 /*################################################################################################*/
35
setup(model * p1)36 setup::setup(model * p1)
37 {
38 mdl = p1;
39
40 // initialize the tables...
41
42 has_setup_tables = false;
43
44 atmtab = NULL; natm = NOT_DEFINED;
45
46 qm_atmtab = NULL; qm_natm = NOT_DEFINED;
47 qm_bndtab = NULL; qm_nbnd = NOT_DEFINED;
48
49 mm_atmtab = NULL; mm_natm = NOT_DEFINED;
50 mm_bndtab = NULL; mm_nbnd = NOT_DEFINED;
51
52 boundary_bndtab = NULL; boundary_nbnd = NOT_DEFINED;
53
54 sf_atmtab = NULL; sf_natm = NOT_DEFINED;
55
56 current_eng = NULL;
57 current_eng_index = 0;
58 }
59
~setup(void)60 setup::~setup(void)
61 {
62 DiscardCurrentEngine();
63 DiscardSetupInfo();
64 }
65
DiscardSetupInfo(void)66 void setup::DiscardSetupInfo(void)
67 {
68 if (atmtab != NULL) { delete[] atmtab; atmtab = NULL;}
69
70 if (qm_atmtab != NULL) { delete[] qm_atmtab; qm_atmtab = NULL; }
71 if (qm_bndtab != NULL) { delete[] qm_bndtab; qm_bndtab = NULL; }
72
73 if (mm_atmtab != NULL) { delete[] mm_atmtab; mm_atmtab = NULL; }
74 if (mm_bndtab != NULL) { delete[] mm_bndtab; mm_bndtab = NULL; }
75
76 if (boundary_bndtab != NULL) { delete[] boundary_bndtab; boundary_bndtab = NULL; }
77
78 if (sf_atmtab != NULL) { delete[] sf_atmtab; sf_atmtab = NULL; }
79
80 has_setup_tables = false;
81 }
82
UpdateSetupInfo(void)83 void setup::UpdateSetupInfo(void)
84 {
85 // discard the previous tables, if needed...
86
87 DiscardSetupInfo();
88
89 // clear all atom flags, and then bring them up-to-date.
90
91 for (iter_al it1 = GetModel()->GetAtomsBegin();it1 != GetModel()->GetAtomsEnd();it1++)
92 {
93 (* it1).flags &= (~ ATOMFLAG_IS_QM_ATOM);
94 (* it1).flags &= (~ ATOMFLAG_IS_MM_ATOM);
95 (* it1).flags &= (~ ATOMFLAG_IS_SF_ATOM);
96 (* it1).flags &= (~ ATOMFLAG_IS_HIDDEN);
97 }
98
99 UpdateAtomFlags();
100
101 // start initializing the tables...
102
103 natm = 0; i32s nbnd = 0;
104
105 qm_natm = 0; qm_nbnd = 0;
106 mm_natm = 0; mm_nbnd = 0;
107 boundary_nbnd = 0;
108
109 sf_natm = 0;
110
111 // count the atoms and bonds that are taken into computations.
112 // count the atoms and bonds that are taken into computations.
113 // count the atoms and bonds that are taken into computations.
114
115 iter_al ita = mdl->GetAtomsBegin();
116 while (ita != mdl->GetAtomsEnd())
117 {
118 bool test1 = ((* ita).flags & (ATOMFLAG_IS_QM_ATOM | ATOMFLAG_IS_MM_ATOM | ATOMFLAG_IS_SF_ATOM));
119 bool test2 = ((* ita).flags & ATOMFLAG_IS_HIDDEN);
120 bool skip = (!test1 || test2);
121
122 if (!skip)
123 {
124 natm++;
125
126 if ((* ita).flags & ATOMFLAG_IS_QM_ATOM)
127 {
128 qm_natm++;
129 }
130
131 if ((* ita).flags & ATOMFLAG_IS_MM_ATOM)
132 {
133 mm_natm++;
134 }
135
136 if ((* ita).flags & ATOMFLAG_IS_SF_ATOM)
137 {
138 sf_natm++;
139 }
140 }
141
142 ita++;
143 }
144
145 iter_bl itb = mdl->GetBondsBegin();
146 while (itb != mdl->GetBondsEnd())
147 {
148 bool testA1 = ((* itb).atmr[0]->flags & (ATOMFLAG_IS_QM_ATOM | ATOMFLAG_IS_MM_ATOM));
149 bool testA2 = ((* itb).atmr[0]->flags & ATOMFLAG_IS_HIDDEN);
150 bool skipA = (!testA1 || testA2);
151
152 bool testB1 = ((* itb).atmr[1]->flags & (ATOMFLAG_IS_QM_ATOM | ATOMFLAG_IS_MM_ATOM));
153 bool testB2 = ((* itb).atmr[1]->flags & ATOMFLAG_IS_HIDDEN);
154 bool skipB = (!testB1 || testB2);
155
156 if (!skipA && !skipB)
157 {
158 nbnd++;
159
160 if (((* itb).atmr[0]->flags & ATOMFLAG_IS_QM_ATOM) && ((* itb).atmr[1]->flags & ATOMFLAG_IS_QM_ATOM))
161 {
162 qm_nbnd++;
163 }
164
165 if (((* itb).atmr[0]->flags & ATOMFLAG_IS_MM_ATOM) && ((* itb).atmr[1]->flags & ATOMFLAG_IS_MM_ATOM))
166 {
167 mm_nbnd++;
168 }
169
170 bool tqmmm1 = (((* itb).atmr[0]->flags & ATOMFLAG_IS_QM_ATOM) && ((* itb).atmr[1]->flags & ATOMFLAG_IS_MM_ATOM));
171 bool tqmmm2 = (((* itb).atmr[0]->flags & ATOMFLAG_IS_MM_ATOM) && ((* itb).atmr[1]->flags & ATOMFLAG_IS_QM_ATOM));
172 if (tqmmm1 || tqmmm2)
173 {
174 boundary_nbnd++;
175 }
176 }
177
178 itb++;
179 }
180
181 // bonds (in contrast to atoms) should only be counted in a single category; check!!!
182
183 if (qm_nbnd + mm_nbnd + boundary_nbnd != nbnd)
184 {
185 ostringstream msg;
186 msg << "bond counting mismatch : " << qm_nbnd << " + " << mm_nbnd << " + " << boundary_nbnd << " != " << nbnd << ends;
187 assertion_failed(__FILE__, __LINE__, msg.str().c_str());
188 }
189
190 // allocate the tables.
191
192 atmtab = new ref_atom[natm];
193
194 qm_atmtab = new ref_atom[qm_natm];
195 qm_bndtab = new ref_bond[qm_nbnd];
196
197 mm_atmtab = new ref_atom[mm_natm];
198 mm_bndtab = new ref_bond[mm_nbnd];
199
200 boundary_bndtab = new ref_bond[boundary_nbnd];
201
202 sf_atmtab = new ref_atom[sf_natm];
203
204 // fill the tables; BE SURE TO USE THE SAME RULES AS IN THE COUNTING STAGE!!!
205 // fill the tables; BE SURE TO USE THE SAME RULES AS IN THE COUNTING STAGE!!!
206 // fill the tables; BE SURE TO USE THE SAME RULES AS IN THE COUNTING STAGE!!!
207
208 i32u ac = 0;
209 i32u qm_ac = 0;
210 i32u mm_ac = 0;
211 i32u sf_ac = 0;
212
213 ita = mdl->GetAtomsBegin();
214 while (ita != mdl->GetAtomsEnd())
215 {
216 bool test1 = ((* ita).flags & (ATOMFLAG_IS_QM_ATOM | ATOMFLAG_IS_MM_ATOM | ATOMFLAG_IS_SF_ATOM));
217 bool test2 = ((* ita).flags & ATOMFLAG_IS_HIDDEN);
218 bool skip = (!test1 || test2);
219
220 if (!skip)
221 {
222 (* ita).varind = ac; // SET THE VARIABLE INDEX!!!
223 atmtab[ac++] = & (* ita);
224
225 if ((* ita).flags & ATOMFLAG_IS_QM_ATOM)
226 {
227 qm_atmtab[qm_ac++] = & (* ita);
228 }
229
230 if ((* ita).flags & ATOMFLAG_IS_MM_ATOM)
231 {
232 mm_atmtab[mm_ac++] = & (* ita);
233 }
234
235 if ((* ita).flags & ATOMFLAG_IS_SF_ATOM)
236 {
237 sf_atmtab[sf_ac++] = & (* ita);
238 }
239 }
240 else (* ita).varind = NOT_DEFINED;
241
242 ita++;
243 }
244
245 i32u qm_bc = 0;
246 i32u mm_bc = 0;
247 i32u boundary_bc = 0;
248
249 itb = mdl->GetBondsBegin();
250 while (itb != mdl->GetBondsEnd())
251 {
252 bool testA1 = ((* itb).atmr[0]->flags & (ATOMFLAG_IS_QM_ATOM | ATOMFLAG_IS_MM_ATOM));
253 bool testA2 = ((* itb).atmr[0]->flags & ATOMFLAG_IS_HIDDEN);
254 bool skipA = (!testA1 || testA2);
255
256 bool testB1 = ((* itb).atmr[1]->flags & (ATOMFLAG_IS_QM_ATOM | ATOMFLAG_IS_MM_ATOM));
257 bool testB2 = ((* itb).atmr[1]->flags & ATOMFLAG_IS_HIDDEN);
258 bool skipB = (!testB1 || testB2);
259
260 if (!skipA && !skipB)
261 {
262 if (((* itb).atmr[0]->flags & ATOMFLAG_IS_QM_ATOM) && ((* itb).atmr[1]->flags & ATOMFLAG_IS_QM_ATOM))
263 {
264 qm_bndtab[qm_bc++] = & (* itb);
265 }
266
267 if (((* itb).atmr[0]->flags & ATOMFLAG_IS_MM_ATOM) && ((* itb).atmr[1]->flags & ATOMFLAG_IS_MM_ATOM))
268 {
269 mm_bndtab[mm_bc++] = & (* itb);
270 }
271
272 bool tqmmm1 = (((* itb).atmr[0]->flags & ATOMFLAG_IS_QM_ATOM) && ((* itb).atmr[1]->flags & ATOMFLAG_IS_MM_ATOM));
273 bool tqmmm2 = (((* itb).atmr[0]->flags & ATOMFLAG_IS_MM_ATOM) && ((* itb).atmr[1]->flags & ATOMFLAG_IS_QM_ATOM));
274 if (tqmmm1 || tqmmm2)
275 {
276 boundary_bndtab[boundary_bc++] = & (* itb);
277 }
278 }
279
280 itb++;
281 }
282
283 // finally check that the counts match!!!
284 // finally check that the counts match!!!
285 // finally check that the counts match!!!
286
287 if ((i32s) ac != natm)
288 {
289 ostringstream msg;
290 msg << "atom mismatch : " << ac << " != " << natm << ends;
291 assertion_failed(__FILE__, __LINE__, msg.str().c_str());
292 }
293
294 if ((i32s) qm_ac != qm_natm)
295 {
296 ostringstream msg;
297 msg << "qm-atom mismatch : " << qm_ac << " != " << qm_natm << ends;
298 assertion_failed(__FILE__, __LINE__, msg.str().c_str());
299 }
300
301 if ((i32s) qm_bc != qm_nbnd)
302 {
303 ostringstream msg;
304 msg << "qm-bond mismatch : " << qm_bc << " != " << qm_nbnd << ends;
305 assertion_failed(__FILE__, __LINE__, msg.str().c_str());
306 }
307
308 if ((i32s) mm_ac != mm_natm)
309 {
310 ostringstream msg;
311 msg << "mm-atom mismatch : " << mm_ac << " != " << mm_natm << ends;
312 assertion_failed(__FILE__, __LINE__, msg.str().c_str());
313 }
314
315 if ((i32s) mm_bc != mm_nbnd)
316 {
317 ostringstream msg;
318 msg << "mm-bond mismatch : " << mm_bc << " != " << mm_nbnd << ends;
319 assertion_failed(__FILE__, __LINE__, msg.str().c_str());
320 }
321
322 if ((i32s) boundary_bc != boundary_nbnd)
323 {
324 ostringstream msg;
325 msg << "boundary-bond mismatch : " << boundary_bc << " != " << boundary_nbnd << ends;
326 assertion_failed(__FILE__, __LINE__, msg.str().c_str());
327 }
328
329 if ((i32s) sf_ac != sf_natm)
330 {
331 ostringstream msg;
332 msg << "sf-atom mismatch : " << sf_ac << " != " << sf_natm << ends;
333 assertion_failed(__FILE__, __LINE__, msg.str().c_str());
334 }
335
336 //////////////////////////////////////////////////
337 //////////////////////////////////////////////////
338 //cout << "atoms : " << qm_ac << " " << mm_ac << " " << sf_ac << " bonds : " << qm_bc << " " << mm_bc << endl; //CGI!
339 //////////////////////////////////////////////////
340 //////////////////////////////////////////////////
341
342 has_setup_tables = true;
343 }
344
CreateCurrentEngine(void)345 void setup::CreateCurrentEngine(void)
346 {
347 DiscardCurrentEngine();
348 current_eng = CreateEngineByIndex(current_eng_index);
349 }
350
DiscardCurrentEngine(void)351 void setup::DiscardCurrentEngine(void)
352 {
353 if (current_eng != NULL)
354 {
355 delete current_eng;
356 current_eng = NULL;
357 }
358 }
359
CreateEngineByIDNumber(i32u id)360 engine * setup::CreateEngineByIDNumber(i32u id)
361 {
362 i32u index = 0;
363 while (index < GetEngineCount())
364 {
365 if (GetEngineIDNumber(index) == id) break;
366 else index++;
367 }
368
369 if (index < GetEngineCount())
370 {
371 return CreateEngineByIndex(index);
372 }
373 else
374 {
375 cout << "WARNING : setup::CreateEngineByIDNumber() failed!" << endl;
376 return NULL;
377 }
378 }
379
380 /*################################################################################################*/
381
engine(setup * p1,i32u p2)382 engine::engine(setup * p1, i32u p2)
383 {
384 stp = p1;
385 if (!stp->HasSetupTables()) assertion_failed(__FILE__, __LINE__, "no setup tables");
386
387 // here it is very important that we take the atom count from setup::GetAtomCount(), not model::GetAtomCount()!!!
388 // it is possible that these are not at all the same!!! some atoms might be hidden from the calculations.
389
390 natm = stp->GetAtomCount();
391
392 crd = new f64[natm * 3];
393
394 if (p2 > 0) d1 = new f64[natm * 3];
395 else d1 = NULL;
396
397 if (p2 > 1) d2 = new f64[natm * natm * 9];
398 else d2 = NULL;
399
400 virial[0] = 0.0;
401 virial[1] = 0.0;
402 virial[2] = 0.0;
403
404 update_neighbor_list = false;
405
406 ECOMPcycles = 0;
407 ECOMPstore_size = 0;
408 ECOMPstore = NULL;
409
410 if (stp->GetModel()->ecomp_enabled)
411 {
412 const int n = stp->GetModel()->ecomp_grp_names.size();
413 ECOMPstore_size = n * (n + 1) / 2;
414
415 ECOMPstore = new ecomp_data[ECOMPstore_size];
416
417 ecomp_Reset();
418 }
419 }
420
~engine(void)421 engine::~engine(void)
422 {
423 delete[] crd;
424 crd = NULL;
425
426 if (d1 != NULL)
427 {
428 delete[] d1;
429 d1 = NULL;
430 }
431
432 if (d2 != NULL)
433 {
434 delete[] d2;
435 d2 = NULL;
436 }
437
438 if (ECOMPstore != NULL)
439 {
440 delete[] ECOMPstore;
441 ECOMPstore = NULL;
442 }
443 }
444
Check(i32u)445 void engine::Check(i32u)
446 {
447 const f64 delta = 0.000001; // the finite difference step...
448
449 Compute(1);
450 f64 tmp1 = energy;
451
452 f64 tmp2; f64 old;
453 for (i32s n1 = 0;n1 < natm;n1++)
454 {
455 for (i32u n2 = 0;n2 < 3;n2++)
456 {
457 old = crd[n1 * 3 + n2];
458 crd[n1 * 3 + n2] = old + delta;
459 Compute(0); tmp2 = (energy - tmp1) / delta;
460 crd[n1 * 3 + n2] = old;
461
462 cout << n1 << ((char) ('x' + n2)) << " ";
463 cout << "a = " << d1[n1 * 3 + n2] << " ";
464 cout << "n = " << tmp2 << endl;
465
466 if ((n1 % 5) == 4) cin >> old;
467 }
468 }
469 }
470
GetGradientVectorLength(void)471 f64 engine::GetGradientVectorLength(void)
472 {
473 f64 sum = 0.0;
474
475 for (i32s n1 = 0;n1 < natm;n1++)
476 {
477 for (i32u n2 = 0;n2 < 3;n2++)
478 {
479 f64 tmp1 = d1[n1 * 3 + n2];
480 sum += tmp1 * tmp1;
481 }
482 }
483
484 return sqrt(sum);
485 }
486
ScaleCRD(f64 kx,f64 ky,f64 kz)487 void engine::ScaleCRD(f64 kx, f64 ky, f64 kz)
488 {
489 atom ** glob_atmtab = GetSetup()->GetAtoms();
490 for (i32s n1 = 0;n1 < natm;n1++)
491 {
492 if (glob_atmtab[n1]->flags & ATOMFLAG_USER_LOCKED) continue;
493
494 crd[n1 * 3 + 0] *= kx;
495 crd[n1 * 3 + 1] *= ky;
496 crd[n1 * 3 + 2] *= kz;
497 }
498 }
499
DoSHAKE(bool)500 void engine::DoSHAKE(bool)
501 {
502 }
503
GetEnergy(void)504 f64 engine::GetEnergy(void)
505 {
506 return energy;
507 }
508
ReadCRD(int i)509 f64 engine::ReadCRD(int i)
510 {
511 if (i < 0) return 0.0;
512 if (i >= natm * 3) return 0.0;
513
514 return crd[i];
515 }
516
ecomp_AddCycle(void)517 void engine::ecomp_AddCycle(void)
518 {
519 ECOMPcycles++;
520 }
521
ecomp_AddStore2(int gA,int gB,int sc,double value)522 void engine::ecomp_AddStore2(int gA, int gB, int sc, double value)
523 {
524 const int gLO = (gA < gB ? gA : gB);
525 const int gHI = (gA > gB ? gA : gB);
526
527 const int index = (gHI * (gHI + 1) / 2) + gLO;
528
529 if (index >= ECOMPstore_size) assertion_failed(__FILE__, __LINE__, "index overflow");
530
531 ECOMPstore[index][sc] += value;
532 }
533
ecomp_AddStoreX(vector<int> & gv,int sc,double value)534 void engine::ecomp_AddStoreX(vector<int> & gv, int sc, double value)
535 {
536 if (gv.size() < 1) return;
537
538 if (gv.size() == 1)
539 {
540 ecomp_AddStore2(gv[0], gv[0], sc, value);
541 return;
542 }
543
544 if (gv.size() == 2)
545 {
546 ecomp_AddStore2(gv[0], gv[1], sc, value);
547 return;
548 }
549
550 vector<int> uniq;
551 uniq.push_back(gv[0]);
552
553 for (int n1 = 1;n1 < (int) gv.size();n1++)
554 {
555 bool is_uniq = true;
556 for (int n2 = 0;n2 < (int) uniq.size();n2++)
557 {
558 if (uniq[n2] == gv[n1])
559 {
560 is_uniq = false;
561 break;
562 }
563 }
564
565 if (is_uniq) uniq.push_back(gv[n1]);
566 }
567
568 if (uniq.size() == 1)
569 {
570 ecomp_AddStore2(uniq[0], uniq[0], sc, value);
571 }
572 else
573 {
574 const int s = uniq.size();
575 const int n = s * (s - 1) / 2;
576 const double valX = value / (double) n;
577
578 for (int x = 0;x < s - 1;x++)
579 {
580 for (int y = x + 1;y < s;y++)
581 {
582 ecomp_AddStore2(uniq[x], uniq[y], sc, valX);
583 }
584 }
585 }
586 }
587
ecomp_Reset(void)588 void engine::ecomp_Reset(void)
589 {
590 ECOMPcycles = 0;
591
592 for (int i = 0;i < ECOMPstore_size;i++)
593 {
594 ECOMPstore[i][ECOMP_DATA_IND_B_bs] = 0.0;
595 ECOMPstore[i][ECOMP_DATA_IND_B_ab] = 0.0;
596 ECOMPstore[i][ECOMP_DATA_IND_B_ti] = 0.0;
597 ECOMPstore[i][ECOMP_DATA_IND_NB_lj] = 0.0;
598 ECOMPstore[i][ECOMP_DATA_IND_NB_es] = 0.0;
599 }
600 }
601
ecomp_ReadStore(int gA,int gB,int sc)602 double engine::ecomp_ReadStore(int gA, int gB, int sc)
603 {
604 const int gLO = (gA < gB ? gA : gB);
605 const int gHI = (gA > gB ? gA : gB);
606
607 const int index = (gHI * (gHI + 1) / 2) + gLO;
608
609 if (index >= ECOMPstore_size) assertion_failed(__FILE__, __LINE__, "index overflow");
610
611 return (ECOMPstore[index][sc] / (double) ECOMPcycles);
612 }
613
614 /*################################################################################################*/
615
value_VDWSurf(engine * eng,fGL * crd,fGL * grad)616 fGL value_VDWSurf(engine * eng, fGL * crd, fGL * grad)
617 {
618 return eng->GetVDWSurf(crd, grad);
619 }
620
value_ESP(engine * eng,fGL * crd,fGL * grad)621 fGL value_ESP(engine * eng, fGL * crd, fGL * grad)
622 {
623 return eng->GetESP(crd, grad);
624 }
625
value_ElDens(engine * eng,fGL * crd,fGL * grad)626 fGL value_ElDens(engine * eng, fGL * crd, fGL * grad)
627 {
628 return eng->GetElDens(crd, grad);
629 }
630
value_Orbital(engine * eng,fGL * crd,fGL * grad)631 fGL value_Orbital(engine * eng, fGL * crd, fGL * grad)
632 {
633 return eng->GetOrbital(crd, grad);
634 }
635
value_OrbDens(engine * eng,fGL * crd,fGL * grad)636 fGL value_OrbDens(engine * eng, fGL * crd, fGL * grad)
637 {
638 return eng->GetOrbDens(crd, grad);
639 }
640
641 // this #include macro is located here, because the eng1_sf class derived from the
642 // engine class and therefore it only can be introduced after the engine class itself...
643
644 #include "eng1_sf.h"
645
CopyCRD(model * p1,engine * p2,i32u p3)646 void CopyCRD(model * p1, engine * p2, i32u p3)
647 {
648 if (p3 >= p1->cs_vector.size()) assertion_failed(__FILE__, __LINE__, "cs overflow");
649
650 atom ** glob_atmtab = p2->GetSetup()->GetAtoms();
651 for (i32s n1 = 0;n1 < p2->GetSetup()->GetAtomCount();n1++)
652 {
653 const fGL * cdata = glob_atmtab[n1]->GetCRD(p3);
654
655 p2->crd[n1 * 3 + 0] = cdata[0];
656 p2->crd[n1 * 3 + 1] = cdata[1];
657 p2->crd[n1 * 3 + 2] = cdata[2];
658 }
659
660 // the rest is SF-related...
661 // the rest is SF-related...
662 // the rest is SF-related...
663
664 eng1_sf * esf = dynamic_cast<eng1_sf *>(p2);
665 setup1_sf * ssf = dynamic_cast<setup1_sf *>(p2->GetSetup());
666 if (esf != NULL && ssf != NULL)
667 {
668 i32s bt3_counter = 0;
669
670 for (i32u cc = 0;cc < ssf->chn_vector.size();cc++)
671 {
672 for (i32s rc = 1;rc < ((i32s) ssf->chn_vector[cc].res_vector.size()) - 2;rc++)
673 {
674 const fGL * prev = ssf->chn_vector[cc].res_vector[rc - 1].atmr[0]->GetCRD(p3);
675 const fGL * curr = ssf->chn_vector[cc].res_vector[rc + 0].atmr[0]->GetCRD(p3);
676 const fGL * next = ssf->chn_vector[cc].res_vector[rc + 1].atmr[0]->GetCRD(p3);
677
678 atom * ref_to_C = ssf->chn_vector[cc].res_vector[rc + 0].peptide[2];
679 atom * ref_to_O = ssf->chn_vector[cc].res_vector[rc + 0].peptide[3];
680
681 v3d<fGL> v1(curr, prev); v3d<fGL> v2(curr, next);
682 v3d<fGL> v3(ref_to_O->GetCRD(p3), ref_to_C->GetCRD(p3));
683 fGL pbdd = v1.tor(v2, v3);
684
685 if (bt3_counter >= (i32s) esf->bt3_vector.size()) assertion_failed(__FILE__, __LINE__, "bt3_counter overflow");
686
687 esf->bt3_vector[bt3_counter++].pbdd = pbdd;
688 }
689 }
690 }
691
692 // ready!!!
693 // ready!!!
694 // ready!!!
695 }
696
CopyCRD(engine * p1,model * p2,i32u p3)697 void CopyCRD(engine * p1, model * p2, i32u p3)
698 {
699 if (p3 >= p2->cs_vector.size()) assertion_failed(__FILE__, __LINE__, "cs overflow");
700
701 atom ** glob_atmtab = p1->GetSetup()->GetAtoms();
702 for (i32s n1 = 0;n1 < p1->GetSetup()->GetAtomCount();n1++)
703 {
704 fGL x = p1->crd[n1 * 3 + 0];
705 fGL y = p1->crd[n1 * 3 + 1];
706 fGL z = p1->crd[n1 * 3 + 2];
707
708 glob_atmtab[n1]->SetCRD(p3, x, y, z);
709 }
710
711 // the rest is SF-related...
712 // the rest is SF-related...
713 // the rest is SF-related...
714
715 eng1_sf * esf = dynamic_cast<eng1_sf *>(p1);
716 setup1_sf * ssf = dynamic_cast<setup1_sf *>(p1->GetSetup());
717 if (esf != NULL && ssf != NULL)
718 {
719 i32s bt3_counter = 0;
720
721 for (i32u cc = 0;cc < ssf->chn_vector.size();cc++)
722 {
723 for (i32s rc = 1;rc < ((i32s) ssf->chn_vector[cc].res_vector.size()) - 2;rc++)
724 {
725 const fGL * prev = ssf->chn_vector[cc].res_vector[rc - 1].atmr[0]->GetCRD(p3);
726 const fGL * curr = ssf->chn_vector[cc].res_vector[rc + 0].atmr[0]->GetCRD(p3);
727 const fGL * next = ssf->chn_vector[cc].res_vector[rc + 1].atmr[0]->GetCRD(p3);
728
729 v3d<fGL> v1(curr, prev); v3d<fGL> v2(curr, next);
730 v3d<fGL> v3 = v1.vpr(v2); v3 = v3 / v3.len(); // this is vector c in the JCC 2001 paper.
731 v3d<fGL> v4 = v2.vpr(v3); v4 = v4 / v4.len(); // this is vector n in the JCC 2001 paper.
732
733 if (bt3_counter >= (i32s) esf->bt3_vector.size()) assertion_failed(__FILE__, __LINE__, "bt3_counter overflow");
734
735 fGL pbdd = esf->bt3_vector[bt3_counter++].pbdd;
736 v3d<fGL> v5 = (v4 * cos(pbdd)) - (v3 * sin(pbdd));
737
738 v2 = v2 / v2.len(); const fGL scale = 0.3785;
739 v3d<fGL> vC = (v2 * (+0.384 * scale)) + (v5 * (-0.116 * scale));
740 v3d<fGL> vO = (v2 * (+0.442 * scale)) + (v5 * (-0.449 * scale));
741 v3d<fGL> vN = (v2 * (+0.638 * scale)) + (v5 * (+0.109 * scale));
742
743 atom * ref_to_C = ssf->chn_vector[cc].res_vector[rc + 0].peptide[2];
744 atom * ref_to_O = ssf->chn_vector[cc].res_vector[rc + 0].peptide[3];
745 atom * ref_to_N = ssf->chn_vector[cc].res_vector[rc + 1].peptide[0];
746
747 fGL x; fGL y; fGL z;
748
749 x = curr[0] + vC[0]; y = curr[1] + vC[1]; z = curr[2] + vC[2]; ref_to_C->SetCRD(p3, x, y, z);
750 x = curr[0] + vO[0]; y = curr[1] + vO[1]; z = curr[2] + vO[2]; ref_to_O->SetCRD(p3, x, y, z);
751 x = curr[0] + vN[0]; y = curr[1] + vN[1]; z = curr[2] + vN[2]; ref_to_N->SetCRD(p3, x, y, z);
752 }
753 }
754 }
755
756 // ready!!!
757 // ready!!!
758 // ready!!!
759 }
760
761 /*################################################################################################*/
762
engine_bp(setup * p1,i32u p2)763 engine_bp::engine_bp(setup * p1, i32u p2) : engine(p1, p2)
764 {
765 use_bp = GetSetup()->GetModel()->use_boundary_potential;
766 bp_rad_solute = GetSetup()->GetModel()->saved_boundary_potential_rad_solute;
767 bp_rad_solvent = GetSetup()->GetModel()->saved_boundary_potential_rad_solvent;
768
769 nd_eval = NULL;
770 rdf_eval = NULL;
771 }
772
~engine_bp(void)773 engine_bp::~engine_bp(void)
774 {
775 if (nd_eval != NULL) delete nd_eval;
776 if (rdf_eval != NULL) delete rdf_eval;
777 }
778
779 /*################################################################################################*/
780
engine_pbc(setup * p1,i32u p2)781 engine_pbc::engine_pbc(setup * p1, i32u p2) : engine(p1, p2)
782 {
783 f64 * tmp = GetSetup()->GetModel()->saved_periodic_box_HALFdim;
784
785 box_HALFdim[0] = tmp[0];
786 box_HALFdim[1] = tmp[1];
787 box_HALFdim[2] = tmp[2];
788
789 tmp = NULL;
790
791 num_mol = 0;
792
793 // count the molecules present in the full atom set ; since the "mol"
794 // level is the highest criteria in sorting, atoms in a molecule should
795 // be adjacent -> a continuous range of pointers.
796
797 if (!GetSetup()->GetModel()->IsGroupsSorted()) assertion_failed(__FILE__, __LINE__, "not_sorted");
798
799 // here we calculate the molecule locations precisely, but a simple trigger atom could be used as well...
800
801 i32s previous = -123; // what is the safest setting here??? NOT_DEFINED might be used there???
802
803 atom ** atmtab = GetSetup()->GetAtoms();
804 for (i32s index = 0;index < GetSetup()->GetAtomCount();index++)
805 {
806 if (atmtab[index]->id[0] != previous)
807 {
808 num_mol++;
809 previous = atmtab[index]->id[0];
810 }
811 }
812
813 mrange = new i32s[num_mol + 1];
814
815 mrange[0] = 0; i32s a_index = 0; // a_index counts LOCAL atom indices.
816 for (i32s n1 = 0;n1 < num_mol;n1++)
817 {
818 i32s m_index = atmtab[a_index]->id[0];
819
820 // m_index counts atom::id[0] molecule numbers.
821 // m_index MAY APPEAR DISCONTINUOUS IF eng1_sf!
822 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
823
824 while (a_index < GetSetup()->GetAtomCount() && atmtab[a_index]->id[0] == m_index) a_index++;
825 mrange[n1 + 1] = a_index;
826 }
827 }
828
~engine_pbc(void)829 engine_pbc::~engine_pbc(void)
830 {
831 delete[] mrange;
832 }
833
CheckLocations(void)834 void engine_pbc::CheckLocations(void)
835 {
836 atom ** atmtab = GetSetup()->GetAtoms();
837 for (i32s n1 = 0;n1 < num_mol;n1++)
838 {
839 f64 sum[3] = { 0.0, 0.0, 0.0 };
840 f64 ac = (f64) (mrange[n1 + 1] - mrange[n1]);
841 for (i32s n2 = mrange[n1];n2 < mrange[n1 + 1];n2++)
842 {
843 i32u index = atmtab[n2]->varind;
844 for (i32s n3 = 0;n3 < 3;n3++)
845 {
846 sum[n3] += crd[index * 3 + n3];
847 }
848 }
849
850 for (i32s n2 = 0;n2 < 3;n2++)
851 {
852 f64 test = sum[n2] / ac;
853
854 if (test < -box_HALFdim[n2])
855 {
856 for (i32s n3 = mrange[n1];n3 < mrange[n1 + 1];n3++)
857 {
858 i32u index = atmtab[n3]->varind;
859 crd[index * 3 + n2] += 2.0 * box_HALFdim[n2];
860 }
861 }
862 else if (test > +box_HALFdim[n2])
863 {
864 for (i32s n3 = mrange[n1];n3 < mrange[n1 + 1];n3++)
865 {
866 i32u index = atmtab[n3]->varind;
867 crd[index * 3 + n2] -= 2.0 * box_HALFdim[n2];
868 }
869 }
870 }
871 }
872 }
873
874 // TODO :
875 // GetVDWSurf() for engine_pbc???
876 // GetESP() for engine_pbc???
877
878 /*################################################################################################*/
879
number_density_evaluator(engine_bp * p1,bool p2,i32s p3)880 number_density_evaluator::number_density_evaluator(engine_bp * p1, bool p2, i32s p3)
881 {
882 eng = p1;
883 linear = p2;
884 classes = p3;
885
886 if (!eng->use_bp) assertion_failed(__FILE__, __LINE__, "use_bp is false");
887
888 upper_limits = new f64[classes];
889 class_volumes = new f64[classes];
890 UpdateClassLimits();
891
892 counts = new i32u[classes + 1];
893 ResetCounters();
894 }
895
~number_density_evaluator(void)896 number_density_evaluator::~number_density_evaluator(void)
897 {
898 delete[] upper_limits;
899 delete[] class_volumes;
900
901 delete[] counts;
902 }
903
UpdateClassLimits(void)904 void number_density_evaluator::UpdateClassLimits(void)
905 {
906 if (linear)
907 {
908 f64 prev_radius = 0.0;
909 for (i32s n1 = 0;n1 < classes;n1++)
910 {
911 f64 next_radius = eng->bp_rad_solvent * ((f64) (n1 + 1)) / (f64) classes;
912 upper_limits[n1] = next_radius;
913
914 f64 volume1 = 4.0 * M_PI * prev_radius * prev_radius * prev_radius / 3.0;
915 f64 volume2 = 4.0 * M_PI * next_radius * next_radius * next_radius / 3.0;
916
917 class_volumes[n1] = volume2 - volume1;
918
919 prev_radius = next_radius; // this is the last operation...
920 }
921 }
922 else
923 {
924 f64 volume1 = 4.0 * M_PI * eng->bp_rad_solvent * eng->bp_rad_solvent * eng->bp_rad_solvent / 3.0;
925 f64 volume2 = volume1 / (f64) classes;
926
927 f64 prev_radius = 0.0;
928 for (i32s n1 = 0;n1 < classes;n1++)
929 {
930 f64 tmp1 = volume2 + 4.0 * M_PI * prev_radius * prev_radius * prev_radius / 3.0;
931 f64 tmp2 = tmp1 / (4.0 * M_PI / 3.0);
932
933 f64 next_radius = pow(tmp2, 1.0 / 3.0);
934
935 upper_limits[n1] = next_radius;
936 class_volumes[n1] = volume2;
937
938 prev_radius = next_radius; // this is the last operation...
939 }
940 }
941 }
942
ResetCounters(void)943 void number_density_evaluator::ResetCounters(void)
944 {
945 cycles = 0;
946 for (i32s n1 = 0;n1 < classes + 1;n1++)
947 {
948 counts[n1] = 0;
949 }
950 }
951
PrintResults(ostream & str)952 void number_density_evaluator::PrintResults(ostream & str)
953 {
954 str << "ND : ";
955 for (i32s n1 = 0;n1 < classes;n1++)
956 {
957 f64 tmp1 = ((f64) counts[n1]) / ((f64) cycles);
958 f64 tmp2 = tmp1 / class_volumes[n1];
959
960 str << tmp2 << " ";
961 }
962
963 f64 tmp1 = ((f64) counts[classes]) / ((f64) cycles);
964 str << _("(outside bp_radius = ") << tmp1 << ")." << endl;
965
966 ResetCounters();
967 }
968
969 /*################################################################################################*/
970
radial_density_function_evaluator(engine_bp * p1,i32s p2,f64 gb,f64 ge,f64 cb,f64 ce)971 radial_density_function_evaluator::radial_density_function_evaluator(engine_bp * p1, i32s p2, f64 gb, f64 ge, f64 cb, f64 ce)
972 {
973 eng = p1;
974 classes = p2;
975
976 graph_begin = gb;
977 graph_end = ge;
978
979 count_begin = gb;
980 count_end = ge;
981
982 if (count_begin < 0.0)
983 {
984 if (!eng->use_bp) assertion_failed(__FILE__, __LINE__, "use_bp is false");
985
986 if (!eng->nd_eval) assertion_failed(__FILE__, __LINE__, "nd_eval is NULL");
987
988 f64 graph_width = graph_end - graph_begin;
989 f64 count_width = count_end - count_begin;
990 if (count_width < graph_width) assertion_failed(__FILE__, __LINE__, "bad width");
991 }
992
993 upper_limits = new f64[classes];
994 class_volumes = new f64[classes];
995
996 f64 prev_radius = graph_begin;
997 for (i32s n1 = 0;n1 < classes;n1++)
998 {
999 f64 next_radius = graph_begin + (graph_end - graph_begin) * ((f64) (n1 + 1)) / (f64) classes;
1000 upper_limits[n1] = next_radius;
1001
1002 f64 volume1 = 4.0 * M_PI * prev_radius * prev_radius * prev_radius / 3.0;
1003 f64 volume2 = 4.0 * M_PI * next_radius * next_radius * next_radius / 3.0;
1004
1005 class_volumes[n1] = volume2 - volume1;
1006
1007 prev_radius = next_radius; // this is the last operation...
1008 }
1009
1010 counts = new i32u[classes];
1011 ResetCounters();
1012 }
1013
~radial_density_function_evaluator(void)1014 radial_density_function_evaluator::~radial_density_function_evaluator(void)
1015 {
1016 delete[] upper_limits;
1017 delete[] class_volumes;
1018
1019 delete[] counts;
1020 }
1021
ResetCounters(void)1022 void radial_density_function_evaluator::ResetCounters(void)
1023 {
1024 cycles = 0;
1025 for (i32s n1 = 0;n1 < classes;n1++)
1026 {
1027 counts[n1] = 0;
1028 }
1029 }
1030
PrintResults(ostream & str)1031 void radial_density_function_evaluator::PrintResults(ostream & str)
1032 {
1033 str << "RDF : ";
1034 for (i32s n1 = 0;n1 < classes;n1++)
1035 {
1036 f64 tmp1 = ((f64) counts[n1]) / ((f64) cycles);
1037 f64 tmp2 = tmp1 / class_volumes[n1];
1038
1039 str << tmp2 << " ";
1040 } str << endl;
1041
1042 ResetCounters();
1043 }
1044
1045 /*################################################################################################*/
1046
1047 // eof
1048