1 // ENG1_MM_DEFAULT.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 "eng1_mm_default.h"
23
24 #include "v3d.h"
25
26 #include "eng1_mm.h"
27 #include "tab_mm_default.h"
28
29 #include "local_i18n.h"
30 #include "notice.h"
31
32 #include <algorithm>
33 #include <sstream>
34 using namespace std;
35
36 /*################################################################################################*/
37
eng1_mm_default_bt(setup * p1,i32u p2)38 eng1_mm_default_bt::eng1_mm_default_bt(setup * p1, i32u p2) :
39 engine(p1, p2),
40 eng1_mm(p1, p2)
41 {
42 if (GetSetup()->GetModel()->GetConstD_count() > 0)
43 {
44 GetSetup()->GetModel()->WarningMessage(CONSTRAINTS_NOT_SUPPORTED);
45 }
46
47 default_tables::GetInstance()->UpdateTypes(GetSetup());
48 default_tables::GetInstance()->UpdateCharges(GetSetup()); // for AMBER stuff, this affects types as well...
49
50 // set default sizes to containers. does it have any practical effects???
51
52 bt1_vector.reserve(250);
53 bt2_vector.reserve(500);
54 bt3_vector.reserve(500);
55
56 atom ** atmtab = GetSetup()->GetMMAtoms();
57 bond ** bndtab = GetSetup()->GetMMBonds();
58
59 atom ** glob_atmtab = GetSetup()->GetAtoms();
60
61 ostream * ostr = NULL; // do not print output.
62 // ostream * ostr = & cout; // print output to cout.
63
64 /*##############################################*/
65 /*##############################################*/
66
67 // create bt1-terms...
68
69 if (ostr != NULL) (* ostr) << _("creating bt1-terms: ");
70 i32s bt1_err = 0;
71
72 // 2008-12-17 ; an old comment suggested this test is obsolete? ; PROBABLY NOT!!! TH
73 if (!GetSetup()->GetModel()->IsIndexClean()) assertion_failed(__FILE__, __LINE__, "is_index_clean");
74
75 for (i32s n1 = 0;n1 < GetSetup()->GetMMBondCount();n1++)
76 {
77 i32s local_ind1 = 0;
78 while (local_ind1 < GetSetup()->GetMMAtomCount())
79 {
80 if (atmtab[local_ind1] == bndtab[n1]->atmr[0]) break;
81 else local_ind1++;
82 }
83 if (local_ind1 >= GetSetup()->GetMMAtomCount()) assertion_failed(__FILE__, __LINE__, "local_ind1");
84
85 i32s local_ind2 = 0;
86 while (local_ind2 < GetSetup()->GetMMAtomCount())
87 {
88 if (atmtab[local_ind2] == bndtab[n1]->atmr[1]) break;
89 else local_ind2++;
90 }
91 if (local_ind2 >= GetSetup()->GetMMAtomCount()) assertion_failed(__FILE__, __LINE__, "local_ind2");
92
93 mm_default_bt1 newbt1;
94 newbt1.atmi[0] = local_ind1; // this is a local index...
95 newbt1.atmi[1] = local_ind2; // this is a local index...
96
97 bool success = false;
98 if (dynamic_cast<setup1_mm *>(GetSetup())->GetExceptions())
99 {
100 i32s bt = bndtab[n1]->bt.GetValue();
101
102 success = default_tables::GetInstance()->e_Init(this, & newbt1, bt);
103 }
104
105 if (success != true)
106 {
107 default_bs_query query; query.strict = false;
108 query.atmtp[0] = atmtab[newbt1.atmi[0]]->atmtp;
109 query.atmtp[1] = atmtab[newbt1.atmi[1]]->atmtp;
110 query.bndtp = bndtab[n1]->bt.GetValue();
111
112 default_tables::GetInstance()->DoParamSearch(& query, GetSetup()->GetModel());
113 if (query.index != NOT_DEFINED) success = true;
114
115 newbt1.opt = query.opt;
116 newbt1.fc = query.fc;
117 }
118
119 // write temporary index information into the bond objects; this is needed below and is just for convenience.
120 // write temporary index information into the bond objects; this is needed below and is just for convenience.
121 // write temporary index information into the bond objects; this is needed below and is just for convenience.
122
123 bndtab[n1]->tmp_bt1_index = bt1_vector.size(); // the bond objects are modified here!!!
124
125 bt1_err += !success;
126 bt1_vector.push_back(newbt1);
127 }
128
129 if (ostr != NULL) (* ostr) << bt1_vector.size() << _(" terms, ") << bt1_err << _(" errors.") << endl;
130
131 bt1data = new mm_bt1_data[bt1_vector.size()];
132
133 /*##############################################*/
134 /*##############################################*/
135
136 // create bt2-terms...
137
138 if (ostr != NULL) (* ostr) << _("creating bt2-terms: ");
139 i32s bt2_err = 0;
140
141 for (i32s n1 = 0;n1 < GetSetup()->GetMMAtomCount();n1++)
142 {
143 if (atmtab[n1]->GetConnRecCount() < 2) continue;
144
145 // central atom is known, now find all possible combinations of bonds...
146
147 bool dir[2];
148 iter_cl ita; iter_cl rnga[2];
149 iter_cl itb; iter_cl rngb[2];
150
151 rngb[1] = atmtab[n1]->GetConnRecsEnd();
152
153 rnga[0] = atmtab[n1]->GetConnRecsBegin();
154 rnga[1] = rngb[1]; rnga[1]--;
155
156 for (ita = rnga[0];ita != rnga[1];ita++)
157 {
158 rngb[0] = ita; rngb[0]++;
159 dir[0] = (atmtab[n1] == (* ita).bndr->atmr[0]);
160
161 for (itb = rngb[0];itb != rngb[1];itb++)
162 {
163 dir[1] = (atmtab[n1] == (* itb).bndr->atmr[0]);
164
165 mm_default_bt2 newbt2;
166 newbt2.index1[0] = (* ita).bndr->tmp_bt1_index; newbt2.dir1[0] = dir[0];
167 newbt2.index1[1] = (* itb).bndr->tmp_bt1_index; newbt2.dir1[1] = dir[1];
168
169 newbt2.atmi[0] = bt1_vector[(* ita).bndr->tmp_bt1_index].get_atmi(1, dir[0]);
170 newbt2.atmi[2] = bt1_vector[(* itb).bndr->tmp_bt1_index].get_atmi(1, dir[1]);
171
172 newbt2.atmi[1] = bt1_vector[(* ita).bndr->tmp_bt1_index].get_atmi(0, dir[0]);
173
174 bool success = false;
175 if (dynamic_cast<setup1_mm *>(GetSetup())->GetExceptions())
176 {
177 i32s bt[2];
178 bt[0] = (* ita).bndr->bt.GetValue();
179 bt[1] = (* itb).bndr->bt.GetValue();
180
181 success = default_tables::GetInstance()->e_Init(this, & newbt2, bt);
182 }
183
184 if (success != true)
185 {
186 default_ab_query query; query.strict = false;
187 query.atmtp[0] = atmtab[newbt2.atmi[0]]->atmtp;
188 query.atmtp[1] = atmtab[newbt2.atmi[1]]->atmtp;
189 query.atmtp[2] = atmtab[newbt2.atmi[2]]->atmtp;
190 query.bndtp[0] = (* ita).bndr->bt.GetValue();
191 query.bndtp[1] = (* itb).bndr->bt.GetValue();
192
193 default_tables::GetInstance()->DoParamSearch(& query, GetSetup()->GetModel());
194 if (query.index != NOT_DEFINED) success = true;
195
196 newbt2.opt = query.opt;
197 newbt2.fc = query.fc;
198 }
199
200 bt2_err += !success;
201 bt2_vector.push_back(newbt2);
202 }
203 }
204 }
205
206 if (ostr != NULL) (* ostr) << bt2_vector.size() << _(" terms, ") << bt2_err << _(" errors.") << endl;
207
208 bt2data = new mm_bt2_data[bt2_vector.size()];
209
210 /*##############################################*/
211 /*##############################################*/
212
213 // create bt3-terms...
214
215 if (ostr != NULL) (* ostr) << _("creating bt3-terms: ");
216 i32s bt3_err = 0;
217
218 for (i32s n1 = 0;n1 < GetSetup()->GetMMAtomCount();n1++)
219 {
220 if (atmtab[n1]->GetConnRecCount() < 2) continue;
221
222 // skip if default geometry is marked linear for this atomtype...
223 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
224 const default_at * tp1 = default_tables::GetInstance()->GetAtomType(atmtab[n1]->atmtp);
225 // if (tp1 == NULL) { cout << "GetAtomType() failed! " << atmtab[n1]->atmtp << endl; exit(EXIT_FAILURE); }
226 if (tp1 != NULL && (tp1->flags & 3) == 1) continue; // linear geometry -> undefined torsion.
227
228 for (iter_cl it2 = atmtab[n1]->GetConnRecsBegin();it2 != atmtab[n1]->GetConnRecsEnd();it2++)
229 {
230 bool another = (atmtab[n1] == (* it2).bndr->atmr[0]);
231 atom * atmr = (* it2).bndr->atmr[another];
232
233 if (atmr->GetConnRecCount() < 2) continue; // no torsion...
234 if (atmr->varind > atmtab[n1]->varind) continue; // skip duplicates...
235
236 // skip if default geometry is marked linear for this atomtype...
237 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
238 const default_at * tp2 = default_tables::GetInstance()->GetAtomType(atmr->atmtp);
239 // if (tp2 == NULL) { cout << "GetAtomType() failed! " << atmr->atmtp << endl; exit(EXIT_FAILURE); }
240 if (tp2 != NULL && (tp2->flags & 3) == 1) continue; // linear geometry -> undefined torsion.
241
242 // central atoms are known, now find all possible combinations of bonds...
243
244 vector<i32s> ind1a; vector<bool> dir1a; // search for the 1st group...
245 for (iter_cl it3 = atmtab[n1]->GetConnRecsBegin();it3 != atmtab[n1]->GetConnRecsEnd();it3++)
246 {
247 if ((* it3) == (* it2)) continue;
248
249 ind1a.push_back((* it3).bndr->tmp_bt1_index);
250 dir1a.push_back(atmtab[n1] == (* it3).bndr->atmr[0]);
251 }
252
253 vector<i32s> ind2a; vector<bool> dir2a;
254 for (i32u n2 = 0;n2 < ind1a.size();n2++)
255 {
256 i32s tmp1 = 0; i32s tmp2 = NOT_DEFINED;
257 while (true)
258 {
259 i32s bt1[2] = { bt2_vector[tmp1].index1[0], bt2_vector[tmp1].index1[1] };
260 i32s bt2[2] = { bt2_vector[tmp1].dir1[0], bt2_vector[tmp1].dir1[1] };
261
262 if (bt1[0] == ind1a[n2] && bt2[0] == dir1a[n2] && bt1[1] == (* it2).bndr->tmp_bt1_index && bt2[1] == another) tmp2 = true;
263 if (bt1[1] == ind1a[n2] && bt2[1] == dir1a[n2] && bt1[0] == (* it2).bndr->tmp_bt1_index && bt2[0] == another) tmp2 = false;
264
265 if (tmp2 < 0) tmp1++; else break;
266 }
267
268 ind2a.push_back(tmp1);
269 dir2a.push_back(tmp2);
270 }
271
272 vector<i32s> ind1b; vector<bool> dir1b; // search for the 2nd group...
273 for (iter_cl it3 = atmr->GetConnRecsBegin();it3 != atmr->GetConnRecsEnd();it3++)
274 {
275 if ((* it3) == (* it2)) continue;
276
277 ind1b.push_back((* it3).bndr->tmp_bt1_index);
278 dir1b.push_back(atmr == (* it3).bndr->atmr[0]);
279 }
280
281 vector<i32s> ind2b; vector<bool> dir2b;
282 for (i32u n2 = 0;n2 < ind1b.size();n2++)
283 {
284 i32s tmp1 = 0; i32s tmp2 = NOT_DEFINED;
285 while (true)
286 {
287 i32s bt1[2] = { bt2_vector[tmp1].index1[0], bt2_vector[tmp1].index1[1] };
288 i32s bt2[2] = { bt2_vector[tmp1].dir1[0], bt2_vector[tmp1].dir1[1] };
289
290 if (bt1[0] == ind1b[n2] && bt2[0] == dir1b[n2] && bt1[1] == (* it2).bndr->tmp_bt1_index && bt2[1] != another) tmp2 = false;
291 if (bt1[1] == ind1b[n2] && bt2[1] == dir1b[n2] && bt1[0] == (* it2).bndr->tmp_bt1_index && bt2[0] != another) tmp2 = true;
292
293 if (tmp2 < 0) tmp1++; else break;
294 }
295
296 ind2b.push_back(tmp1);
297 dir2b.push_back(tmp2);
298 }
299
300 // now finally create the terms!!!!!
301
302 for (i32u n2 = 0;n2 < ind2a.size();n2++)
303 {
304 for (i32u n3 = 0;n3 < ind2b.size();n3++)
305 {
306 mm_default_bt3 newbt3;
307 newbt3.index2[0] = ind2a[n2];
308 newbt3.index2[1] = ind2b[n3];
309
310 newbt3.index1[0] = bt2_vector[newbt3.index2[0]].get_index(0, dir2a[n2]);
311 newbt3.dir1[0] = bt2_vector[newbt3.index2[0]].get_dir(0, dir2a[n2]);
312
313 newbt3.index1[1] = bt2_vector[newbt3.index2[0]].get_index(1, dir2a[n2]);
314 newbt3.dir1[1] = bt2_vector[newbt3.index2[0]].get_dir(1, dir2a[n2]);
315
316 newbt3.index1[2] = bt2_vector[newbt3.index2[1]].get_index(0, dir2b[n3]);
317 newbt3.dir1[2] = bt2_vector[newbt3.index2[1]].get_dir(0, dir2b[n3]);
318
319 newbt3.index1[3] = bt2_vector[newbt3.index2[1]].get_index(1, dir2b[n3]);
320 newbt3.dir1[3] = bt2_vector[newbt3.index2[1]].get_dir(1, dir2b[n3]);
321
322 newbt3.atmi[0] = bt1_vector[newbt3.index1[0]].get_atmi(1, newbt3.dir1[0]);
323 newbt3.atmi[1] = bt1_vector[newbt3.index1[0]].get_atmi(0, newbt3.dir1[0]);
324 newbt3.atmi[2] = bt1_vector[newbt3.index1[3]].get_atmi(0, newbt3.dir1[3]);
325 newbt3.atmi[3] = bt1_vector[newbt3.index1[3]].get_atmi(1, newbt3.dir1[3]);
326
327 // easiest way to get bondtypes is to find the bonds...
328 // easiest way to get bondtypes is to find the bonds...
329 // easiest way to get bondtypes is to find the bonds...
330
331 bond tmpb1 = bond(atmtab[newbt3.atmi[0]], atmtab[newbt3.atmi[1]], bondtype());
332 iter_bl itb1 = find(GetSetup()->GetModel()->GetBondsBegin(), GetSetup()->GetModel()->GetBondsEnd(), tmpb1);
333
334 bond tmpb2 = bond(atmtab[newbt3.atmi[1]], atmtab[newbt3.atmi[2]], bondtype());
335 iter_bl itb2 = find(GetSetup()->GetModel()->GetBondsBegin(), GetSetup()->GetModel()->GetBondsEnd(), tmpb2);
336
337 bond tmpb3 = bond(atmtab[newbt3.atmi[2]], atmtab[newbt3.atmi[3]], bondtype());
338 iter_bl itb3 = find(GetSetup()->GetModel()->GetBondsBegin(), GetSetup()->GetModel()->GetBondsEnd(), tmpb3);
339
340 bool success = false;
341 if (dynamic_cast<setup1_mm *>(GetSetup())->GetExceptions())
342 {
343 i32s bt[3];
344 bt[0] = (* itb1).bt.GetValue();
345 bt[1] = (* itb2).bt.GetValue();
346 bt[2] = (* itb3).bt.GetValue();
347
348 success = default_tables::GetInstance()->e_Init(this, & newbt3, bt);
349 }
350
351 if (success != true)
352 {
353 default_tr_query query; query.strict = false;
354 query.atmtp[0] = atmtab[newbt3.atmi[0]]->atmtp;
355 query.atmtp[1] = atmtab[newbt3.atmi[1]]->atmtp;
356 query.atmtp[2] = atmtab[newbt3.atmi[2]]->atmtp;
357 query.atmtp[3] = atmtab[newbt3.atmi[3]]->atmtp;
358
359 query.bndtp[0] = (* itb1).bt.GetValue();
360 query.bndtp[1] = (* itb2).bt.GetValue();
361 query.bndtp[2] = (* itb3).bt.GetValue();
362
363 default_tables::GetInstance()->DoParamSearch(& query, GetSetup()->GetModel());
364 if (query.index != NOT_DEFINED) success = true;
365
366 newbt3.fc1 = query.fc1;
367 newbt3.fc2 = query.fc2;
368 newbt3.fc3 = query.fc3;
369 newbt3.fc4 = 0.0;
370 }
371
372 newbt3.constraint = false;
373
374 bt3_err += !success;
375 bt3_vector.push_back(newbt3);
376 }
377 }
378 }
379 }
380
381 if (ostr != NULL) (* ostr) << bt3_vector.size() << _(" terms, ")<< bt3_err << _(" errors.") << endl;
382
383 /*##############################################*/
384 /*##############################################*/
385
386 // create bt4-terms...
387
388 if (ostr != NULL) (* ostr) << _("creating bt4-terms: ");
389 i32s bt4_err = 0;
390
391 for (i32s n1 = 0;n1 < GetSetup()->GetMMAtomCount();n1++)
392 {
393 const default_at * tp = default_tables::GetInstance()->GetAtomType(atmtab[n1]->atmtp);
394 if (!tp || (tp != NULL && !(tp->flags & 16))) continue;
395
396 if (atmtab[n1]->GetConnRecCount() != 3) continue;
397
398 crec * crtab[3];
399 iter_cl iter = atmtab[n1]->GetConnRecsBegin();
400 crtab[0] = & (* iter); iter++;
401 crtab[1] = & (* iter); iter++;
402 crtab[2] = & (* iter);
403
404 for (i32s n2 = 0;n2 < 3;n2++)
405 {
406 mm_default_bt4 newbt4;
407
408 newbt4.index1[0] = crtab[(n2 + 0) % 3]->bndr->tmp_bt1_index;
409 newbt4.dir1[0] = (bt1_vector[newbt4.index1[0]].get_atmi(0, true) == n1);
410
411 newbt4.index1[1] = crtab[(n2 + 1) % 3]->bndr->tmp_bt1_index;
412 newbt4.dir1[1] = (bt1_vector[newbt4.index1[1]].get_atmi(0, true) == n1);
413
414 newbt4.index1[2] = crtab[(n2 + 2) % 3]->bndr->tmp_bt1_index;
415 newbt4.dir1[2] = (bt1_vector[newbt4.index1[2]].get_atmi(0, true) == n1);
416
417 newbt4.index2 = 0;
418 while (newbt4.index2 < (i32s) bt2_vector.size())
419 {
420 bool test[4];
421
422 newbt4.dir2 = true;
423 test[0] = (bt2_vector[newbt4.index2].get_index(0, newbt4.dir2) == newbt4.index1[0]);
424 test[1] = (bt2_vector[newbt4.index2].get_dir(0, newbt4.dir2) == newbt4.dir1[0]);
425 test[2] = (bt2_vector[newbt4.index2].get_index(1, newbt4.dir2) == newbt4.index1[1]);
426 test[3] = (bt2_vector[newbt4.index2].get_dir(1, newbt4.dir2) == newbt4.dir1[1]);
427 if (test[0] && test[1] && test[2] && test[3]) break;
428
429 newbt4.dir2 = false;
430 test[0] = (bt2_vector[newbt4.index2].get_index(0, newbt4.dir2) == newbt4.index1[0]);
431 test[1] = (bt2_vector[newbt4.index2].get_dir(0, newbt4.dir2) == newbt4.dir1[0]);
432 test[2] = (bt2_vector[newbt4.index2].get_index(1, newbt4.dir2) == newbt4.index1[1]);
433 test[3] = (bt2_vector[newbt4.index2].get_dir(1, newbt4.dir2) == newbt4.dir1[1]);
434 if (test[0] && test[1] && test[2] && test[3]) break;
435
436 newbt4.index2++;
437 }
438
439 if (newbt4.index2 == (i32s) bt2_vector.size()) assertion_failed(__FILE__, __LINE__, "index2 overflow");
440
441 i32s iglob[4] =
442 {
443 crtab[(n2 + 0) % 3]->atmr->varind,
444 atmtab[n1]->varind,
445 crtab[(n2 + 1) % 3]->atmr->varind,
446 crtab[(n2 + 2) % 3]->atmr->varind
447 };
448
449 for (i32s n9 = 0;n9 < 4;n9++)
450 {
451 i32s index = 0;
452 while (index < GetSetup()->GetMMAtomCount())
453 {
454 if (glob_atmtab[iglob[n9]] == atmtab[index]) break;
455 else index++;
456 }
457
458 if (index >= GetSetup()->GetMMAtomCount()) assertion_failed(__FILE__, __LINE__, "iloc search failed");
459
460 newbt4.atmi[n9] = index;
461 }
462
463 if (newbt4.atmi[1] != n1) assertion_failed(__FILE__, __LINE__, "atmi[1] != n1");
464
465 bool success = false;
466 if (dynamic_cast<setup1_mm *>(GetSetup())->GetExceptions())
467 {
468 i32s bt[3];
469 bt[0] = crtab[(n2 + 0) % 3]->bndr->bt.GetValue();
470 bt[1] = crtab[(n2 + 1) % 3]->bndr->bt.GetValue();
471 bt[2] = crtab[(n2 + 2) % 3]->bndr->bt.GetValue();
472
473 success = default_tables::GetInstance()->e_Init(this, & newbt4, bt);
474 }
475
476 if (success != true)
477 {
478 default_op_query query; query.strict = false;
479 query.atmtp[0] = atmtab[newbt4.atmi[0]]->atmtp;
480 query.atmtp[1] = atmtab[newbt4.atmi[1]]->atmtp;
481 query.atmtp[2] = atmtab[newbt4.atmi[2]]->atmtp;
482 query.atmtp[3] = atmtab[newbt4.atmi[3]]->atmtp;
483
484 query.bndtp[0] = crtab[(n2 + 0) % 3]->bndr->bt.GetValue();
485 query.bndtp[1] = crtab[(n2 + 1) % 3]->bndr->bt.GetValue();
486 query.bndtp[2] = crtab[(n2 + 2) % 3]->bndr->bt.GetValue();
487
488 default_tables::GetInstance()->DoParamSearch(& query, GetSetup()->GetModel());
489 if (query.index != NOT_DEFINED) success = true;
490
491 newbt4.opt = query.opt;
492 newbt4.fc = query.fc;
493 }
494
495 bt4_err += !success;
496 bt4_vector.push_back(newbt4);
497 }
498 }
499
500 if (ostr != NULL) (* ostr) << bt4_vector.size() << _(" terms, ") << bt4_err << _(" errors.") << endl;
501
502 /*##############################################*/
503 /*##############################################*/
504
505 // report possible errors...
506
507 i32s total_err = bt1_err + bt2_err + bt3_err + bt4_err;
508 if (total_err && GetSetup()->GetModel()->verbosity >= 2)
509 {
510 ostringstream str;
511 str << _("WARNING : there were ") << total_err << _(" missing parameters in the bonded terms.") << endl << ends;
512 GetSetup()->GetModel()->PrintToLog(str.str().c_str());
513 }
514 }
515
~eng1_mm_default_bt(void)516 eng1_mm_default_bt::~eng1_mm_default_bt(void)
517 {
518 delete[] bt1data;
519 delete[] bt2data;
520 }
521
FindTorsion(atom * a1,atom * a2,atom * a3,atom * a4)522 i32s eng1_mm_default_bt::FindTorsion(atom * a1, atom * a2, atom * a3, atom * a4)
523 {
524 i32s iglob[4] = { a1->varind, a2->varind, a3->varind, a4->varind };
525 i32s iloc[4];
526
527 atom ** atmtab = GetSetup()->GetMMAtoms();
528 atom ** glob_atmtab = GetSetup()->GetAtoms();
529
530 for (i32s n1 = 0;n1 < 4;n1++)
531 {
532 i32s index = 0;
533 while (index < GetSetup()->GetMMAtomCount())
534 {
535 if (glob_atmtab[iglob[n1]] == atmtab[index]) break;
536 else index++;
537 }
538
539 if (index >= GetSetup()->GetMMAtomCount()) assertion_failed(__FILE__, __LINE__, "iloc search failed");
540
541 iloc[n1] = index;
542 }
543
544 for (i32s n1 = 0;n1 < (i32s) bt3_vector.size();n1++)
545 {
546 bool test; // since torsion is the same in both directions, we can ignore direction...
547
548 test = true;
549 if (bt3_vector[n1].atmi[0] != iloc[0]) test = false;
550 if (bt3_vector[n1].atmi[1] != iloc[1]) test = false;
551 if (bt3_vector[n1].atmi[2] != iloc[2]) test = false;
552 if (bt3_vector[n1].atmi[3] != iloc[3]) test = false;
553 if (test) return n1;
554
555 test = true;
556 if (bt3_vector[n1].atmi[0] != iloc[3]) test = false;
557 if (bt3_vector[n1].atmi[1] != iloc[2]) test = false;
558 if (bt3_vector[n1].atmi[2] != iloc[1]) test = false;
559 if (bt3_vector[n1].atmi[3] != iloc[0]) test = false;
560 if (test) return n1;
561 }
562
563 return NOT_DEFINED;
564 }
565
SetTorsionConstraint(atom * a1,atom * a2,atom * a3,atom * a4,f64 opt,f64 fc,bool lock_local_structure)566 bool eng1_mm_default_bt::SetTorsionConstraint(atom * a1, atom * a2, atom * a3, atom * a4, f64 opt, f64 fc, bool lock_local_structure)
567 {
568 const i32s ind_bt3 = FindTorsion(a1, a2, a3, a4);
569
570 if (ind_bt3 < 0) return false;
571 if (ind_bt3 >= (i32s) bt3_vector.size()) return false;
572
573 // check that opt is in valid range [-pi,+pi]!!!
574
575 while (opt > +M_PI) opt -= 2.0 * M_PI;
576 while (opt < -M_PI) opt += 2.0 * M_PI;
577
578 // measure the current torsion and set constraints for the other torsions, if requested...
579
580 if (lock_local_structure)
581 {
582 v3d<f64> v1a(& crd[l2g_mm[bt3_vector[ind_bt3].atmi[1]] * 3], & crd[l2g_mm[bt3_vector[ind_bt3].atmi[0]] * 3]);
583 v3d<f64> v1b(& crd[l2g_mm[bt3_vector[ind_bt3].atmi[1]] * 3], & crd[l2g_mm[bt3_vector[ind_bt3].atmi[2]] * 3]);
584 v3d<f64> v1c(& crd[l2g_mm[bt3_vector[ind_bt3].atmi[2]] * 3], & crd[l2g_mm[bt3_vector[ind_bt3].atmi[3]] * 3]);
585 f64 delta = opt - v1a.tor(v1b, v1c);
586
587 while (delta > +M_PI) delta -= 2.0 * M_PI;
588 while (delta < -M_PI) delta += 2.0 * M_PI;
589
590 i32s tmp1 = bt3_vector[ind_bt3].atmi[1];
591 i32s tmp2 = bt3_vector[ind_bt3].atmi[2];
592
593 for (i32s n1 = 0;n1 < (i32s) bt3_vector.size();n1++)
594 {
595 bool test1 = true;
596 if (bt3_vector[n1].atmi[1] != tmp1) test1 = false;
597 if (bt3_vector[n1].atmi[2] != tmp2) test1 = false;
598
599 bool test2 = true;
600 if (bt3_vector[n1].atmi[1] != tmp2) test2 = false;
601 if (bt3_vector[n1].atmi[2] != tmp1) test2 = false;
602
603 if (test1 || test2)
604 {
605 v3d<f64> v2a(& crd[l2g_mm[bt3_vector[n1].atmi[1]] * 3], & crd[l2g_mm[bt3_vector[n1].atmi[0]] * 3]);
606 v3d<f64> v2b(& crd[l2g_mm[bt3_vector[n1].atmi[1]] * 3], & crd[l2g_mm[bt3_vector[n1].atmi[2]] * 3]);
607 v3d<f64> v2c(& crd[l2g_mm[bt3_vector[n1].atmi[2]] * 3], & crd[l2g_mm[bt3_vector[n1].atmi[3]] * 3]);
608 f64 local = v2a.tor(v2b, v2c) + delta;
609
610 while (local > +M_PI) local -= 2.0 * M_PI;
611 while (local < -M_PI) local += 2.0 * M_PI;
612
613 bt3_vector[n1].constraint = true;
614 bt3_vector[n1].fc1 = local; // fc1 = opt
615 bt3_vector[n1].fc2 = fc; // fc2 = fc
616 }
617 }
618 }
619
620 // then set the requested constraint... if lock_local_structure was true, then
621 // this operation is in fact redundant, but perhaps a bit more accurate than above.
622
623 bt3_vector[ind_bt3].constraint = true;
624 bt3_vector[ind_bt3].fc1 = opt; // fc1 = opt
625 bt3_vector[ind_bt3].fc2 = fc; // fc2 = fc
626
627 return true;
628 }
629
RemoveTorsionConstraint(atom * a1,atom * a2,atom * a3,atom * a4)630 bool eng1_mm_default_bt::RemoveTorsionConstraint(atom * a1, atom * a2, atom * a3, atom * a4)
631 {
632 // this is not yet properly implemented...
633 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
634 // what to do if lock_local_structure was set??? then must unlock all...
635
636 return false;
637 }
638
ComputeBT1(i32u p1)639 void eng1_mm_default_bt::ComputeBT1(i32u p1)
640 {
641 energy_bt1 = 0.0;
642
643 atom ** atmtab = GetSetup()->GetMMAtoms();
644
645 // len -> length of the bond vector, in nanometers [nm].
646
647 // the bond is a vector, since it has unique begin and end points.
648 // if a bond vector needs to be reversed, it's begin and end points are swapped.
649 // all data for both forward and reverse vectors are calculated and stored...
650
651 // dlen[0] -> grad[0-2]: for atom 1 when direction = 0, for atom 0 when direction = 1
652 // dlen[1] -> grad[0-2]: for atom 0 when direction = 0, for atom 1 when direction = 1
653
654 // the end point gradient of a vector is always the same as components of the unit vector!!!
655
656 for (i32s n1 = 0;n1 < (i32s) bt1_vector.size();n1++)
657 {
658 i32s * atmi = bt1_vector[n1].atmi;
659
660 f64 t1a[3]; f64 t1b = 0.0;
661 for (i32s n2 = 0;n2 < 3;n2++)
662 {
663 f64 t9a = crd[l2g_mm[atmi[0]] * 3 + n2];
664 f64 t9b = crd[l2g_mm[atmi[1]] * 3 + n2];
665
666 t1a[n2] = t9a - t9b;
667 t1b += t1a[n2] * t1a[n2];
668 }
669
670 f64 t1c = sqrt(t1b);
671 bt1data[n1].len = t1c;
672 for (i32s n2 = 0;n2 < 3;n2++)
673 {
674 f64 t9a = t1a[n2] / t1c;
675
676 bt1data[n1].dlen[0][n2] = +t9a;
677 bt1data[n1].dlen[1][n2] = -t9a;
678 }
679
680 // f = a(x-b)^2
681 // df/dx = 2a(x-b)
682
683 f64 t2a = t1c - bt1_vector[n1].opt;
684 f64 t2b = bt1_vector[n1].fc * t2a * t2a;
685
686 energy_bt1 += t2b;
687
688 if (ECOMPstore != NULL)
689 {
690 const int iA = atmtab[atmi[0]]->ecomp_grp_i;
691 const int iB = atmtab[atmi[1]]->ecomp_grp_i;
692
693 ecomp_AddStore2(iA, iB, ECOMP_DATA_IND_B_bs, t2b);
694 }
695
696 if (p1 > 0)
697 {
698 f64 t2c = 2.0 * bt1_vector[n1].fc * t2a;
699
700 for (i32s n2 = 0;n2 < 3;n2++)
701 {
702 f64 t2d = bt1data[n1].dlen[0][n2] * t2c;
703
704 d1[l2g_mm[atmi[0]] * 3 + n2] += t2d;
705 d1[l2g_mm[atmi[1]] * 3 + n2] -= t2d;
706
707 if (do_virial)
708 {
709 //f64 crdA = crd[l2g_mm[atmi[0]] * 3 + n2];
710 //f64 crdB = crd[l2g_mm[atmi[1]] * 3 + n2];
711 //f64 rAB = crdA - crdB;
712 ///////////////////////////////////////////
713 f64 rAB = t1a[n2];
714
715 virial[n2] -= rAB * t2d; // F = -dE/dr
716 }
717 }
718 }
719 }
720 }
721
ComputeBT2(i32u p1)722 void eng1_mm_default_bt::ComputeBT2(i32u p1)
723 {
724 energy_bt2 = 0.0;
725
726 atom ** atmtab = GetSetup()->GetMMAtoms();
727
728 // ang -> cosine of the bond angle, in the usual range [-1.0, +1.0]
729
730 // we need directions also here... the angle consists of three points, say A-B-C.
731 // when we reverse the angle, we will swap the end points: now they will be C-B-A.
732
733 // dang[0] -> grad[0-2]: for atom 2 when direction = 0, for atom 0 when direction = 1
734 // dang[1] -> grad[0-2]: for atom 1 when direction = 0, for atom 1 when direction = 1
735 // dang[2] -> grad[0-2]: for atom 0 when direction = 0, for atom 2 when direction = 1
736
737 for (i32s n1 = 0;n1 < (i32s) bt2_vector.size();n1++)
738 {
739 i32s * atmi = bt2_vector[n1].atmi;
740
741 i32s * index1 = bt2_vector[n1].index1;
742 bool * dir1 = bt2_vector[n1].dir1;
743
744 f64 * t1a = bt1data[index1[0]].dlen[dir1[0]];
745 f64 * t1b = bt1data[index1[1]].dlen[dir1[1]];
746
747 f64 t1c = t1a[0] * t1b[0] + t1a[1] * t1b[1] + t1a[2] * t1b[2];
748
749 if (t1c < -1.0) t1c = -1.0; // domain check...
750 if (t1c > +1.0) t1c = +1.0; // domain check...
751
752 bt2data[n1].csa = t1c;
753
754 for (i32s n2 = 0;n2 < 3;n2++)
755 {
756 f64 t9a = (t1b[n2] - t1c * t1a[n2]) / bt1data[index1[0]].len;
757 f64 t9b = (t1a[n2] - t1c * t1b[n2]) / bt1data[index1[1]].len;
758
759 bt2data[n1].dcsa[0][n2] = t9a;
760 bt2data[n1].dcsa[1][n2] = -(t9a + t9b);
761 bt2data[n1].dcsa[2][n2] = t9b;
762 }
763
764 f64 t2a; f64 t2b; // t2a = f ; t2b = df/dx
765 if (bt2_vector[n1].opt > NEAR_LINEAR_LIMIT)
766 {
767 // f = a(1 + x)
768 // df/dx = a
769
770 f64 t3b = bt2_vector[n1].fc;
771 t2a = t3b * (1.0 + t1c);
772
773 t2b = t3b;
774 }
775 else
776 {
777 // f = a(acos(x)-b)^2
778 // df/dx = -2a(x-b)/sqrt(1-x*x)
779
780 f64 t3b = acos(t1c) - bt2_vector[n1].opt;
781 t2a = bt2_vector[n1].fc * t3b * t3b;
782
783 t2b = -2.0 * bt2_vector[n1].fc * t3b / sqrt(1.0 - t1c * t1c);
784 }
785
786 energy_bt2 += t2a;
787
788 if (ECOMPstore != NULL)
789 {
790 vector<int> gv;
791 gv.push_back(atmtab[atmi[0]]->ecomp_grp_i);
792 gv.push_back(atmtab[atmi[1]]->ecomp_grp_i);
793 gv.push_back(atmtab[atmi[2]]->ecomp_grp_i);
794
795 ecomp_AddStoreX(gv, ECOMP_DATA_IND_B_ab, t2a);
796 }
797
798 if (p1 > 0)
799 {
800 for (i32s n2 = 0;n2 < 3;n2++)
801 {
802 d1[l2g_mm[atmi[0]] * 3 + n2] += bt2data[n1].dcsa[0][n2] * t2b;
803 d1[l2g_mm[atmi[1]] * 3 + n2] += bt2data[n1].dcsa[1][n2] * t2b;
804 d1[l2g_mm[atmi[2]] * 3 + n2] += bt2data[n1].dcsa[2][n2] * t2b;
805 }
806 }
807 }
808 }
809
ComputeBT3(i32u p1)810 void eng1_mm_default_bt::ComputeBT3(i32u p1)
811 {
812 energy_bt3 = 0.0;
813
814 atom ** atmtab = GetSetup()->GetMMAtoms();
815
816 for (i32s n1 = 0;n1 < (i32s) bt3_vector.size();n1++)
817 {
818 i32s * atmi = bt3_vector[n1].atmi;
819
820 i32s * index2 = bt3_vector[n1].index2;
821 i32s * index1 = bt3_vector[n1].index1;
822 bool * dir1 = bt3_vector[n1].dir1;
823
824 f64 t1a[2] = { bt2data[index2[0]].csa, bt2data[index2[1]].csa };
825 f64 t1b[2] = { 1.0 - t1a[0] * t1a[0], 1.0 - t1a[1] * t1a[1] };
826
827 f64 t1c[2][3];
828 t1c[0][0] = bt1data[index1[0]].dlen[dir1[0]][0] - t1a[0] * bt1data[index1[1]].dlen[dir1[1]][0];
829 t1c[0][1] = bt1data[index1[0]].dlen[dir1[0]][1] - t1a[0] * bt1data[index1[1]].dlen[dir1[1]][1];
830 t1c[0][2] = bt1data[index1[0]].dlen[dir1[0]][2] - t1a[0] * bt1data[index1[1]].dlen[dir1[1]][2];
831 t1c[1][0] = bt1data[index1[3]].dlen[dir1[3]][0] - t1a[1] * bt1data[index1[2]].dlen[dir1[2]][0];
832 t1c[1][1] = bt1data[index1[3]].dlen[dir1[3]][1] - t1a[1] * bt1data[index1[2]].dlen[dir1[2]][1];
833 t1c[1][2] = bt1data[index1[3]].dlen[dir1[3]][2] - t1a[1] * bt1data[index1[2]].dlen[dir1[2]][2];
834
835 f64 t1d = t1c[0][0] * t1c[1][0] + t1c[0][1] * t1c[1][1] + t1c[0][2] * t1c[1][2];
836 f64 t1e = t1d / sqrt(t1b[0] * t1b[1]);
837
838 if (t1e < -1.0) t1e = -1.0; // domain check...
839 if (t1e > +1.0) t1e = +1.0; // domain check...
840
841 f64 t1f[4];
842 t1f[0] = acos(t1e);
843
844 // now we still have to determine the sign of the result...
845 // now we still have to determine the sign of the result...
846 // now we still have to determine the sign of the result...
847
848 f64 t1g[3];
849 t1g[0] = bt1data[index1[2]].dlen[dir1[2]][1] * bt1data[index1[3]].dlen[dir1[3]][2] -
850 bt1data[index1[2]].dlen[dir1[2]][2] * bt1data[index1[3]].dlen[dir1[3]][1];
851 t1g[1] = bt1data[index1[2]].dlen[dir1[2]][2] * bt1data[index1[3]].dlen[dir1[3]][0] -
852 bt1data[index1[2]].dlen[dir1[2]][0] * bt1data[index1[3]].dlen[dir1[3]][2];
853 t1g[2] = bt1data[index1[2]].dlen[dir1[2]][0] * bt1data[index1[3]].dlen[dir1[3]][1] -
854 bt1data[index1[2]].dlen[dir1[2]][1] * bt1data[index1[3]].dlen[dir1[3]][0];
855
856 f64 t1h = t1c[0][0] * t1g[0] + t1c[0][1] * t1g[1] + t1c[0][2] * t1g[2];
857 if (t1h < 0.0) t1f[0] = -t1f[0];
858
859 t1f[1] = t1f[0] + t1f[0]; // 2x
860 t1f[2] = t1f[1] + t1f[0]; // 3x
861 t1f[3] = t1f[2] + t1f[0]; // 4x
862
863 f64 t9a; f64 t9b;
864 if (bt3_vector[n1].constraint)
865 {
866 // Dx = x-b | for the following formulas...
867
868 // f = a(2PI-Dx)^4 | if Dx > +PI
869 // df/fx = -4a(2PI-Dx)^3
870
871 // f = a(2PI+Dx)^4 | if Dx < -PI
872 // df/fx = +4a(2PI+Dx)^3
873
874 // f = a(Dx)^4 | otherwise
875 // df/fx = 4a(Dx)^3
876
877 f64 t8a = t1f[0] - bt3_vector[n1].fc1;
878 if (t8a > +M_PI)
879 {
880 t8a = 2.0 * M_PI - t8a;
881 f64 t8b = t8a * t8a;
882
883 t9a = bt3_vector[n1].fc2 * t8b * t8b;
884 t9b = -4.0 * bt3_vector[n1].fc2 * t8b * t8a;
885 }
886 else if (t8a < -M_PI)
887 {
888 t8a = 2.0 * M_PI + t8a;
889 f64 t8b = t8a * t8a;
890
891 t9a = bt3_vector[n1].fc2 * t8b * t8b;
892 t9b = +4.0 * bt3_vector[n1].fc2 * t8b * t8a;
893 }
894 else
895 {
896 f64 t8b = t8a * t8a;
897
898 t9a = bt3_vector[n1].fc2 * t8b * t8b;
899 t9b = 4.0 * bt3_vector[n1].fc2 * t8b * t8a;
900 }
901 }
902 else
903 {
904 // f = a*cos(x)+b*cos(2x)+c*cos(3x) +d*cos(4x)
905 // df/dx = -a*sin(x)-2b*sin(2x)-3c*sin(3x) -4d*sin(4x)
906
907 f64 t8a = bt3_vector[n1].fc1 * cos(t1f[0]);
908 f64 t8b = bt3_vector[n1].fc2 * cos(t1f[1]);
909 f64 t8c = bt3_vector[n1].fc3 * cos(t1f[2]);
910 f64 t8d = bt3_vector[n1].fc4 * cos(t1f[3]);
911 t9a = t8a + t8b + t8c + t8d;
912
913 f64 t8r = bt3_vector[n1].fc1 * sin(t1f[0]);
914 f64 t8s = bt3_vector[n1].fc2 * sin(t1f[1]) * 2.0;
915 f64 t8t = bt3_vector[n1].fc3 * sin(t1f[2]) * 3.0;
916 f64 t8u = bt3_vector[n1].fc4 * sin(t1f[3]) * 4.0;
917 t9b = -(t8r + t8s + t8t + t8u);
918 }
919
920 energy_bt3 += t9a;
921
922 if (ECOMPstore != NULL)
923 {
924 vector<int> gv;
925 gv.push_back(atmtab[atmi[0]]->ecomp_grp_i);
926 gv.push_back(atmtab[atmi[1]]->ecomp_grp_i);
927 gv.push_back(atmtab[atmi[2]]->ecomp_grp_i);
928 gv.push_back(atmtab[atmi[3]]->ecomp_grp_i);
929
930 ecomp_AddStoreX(gv, ECOMP_DATA_IND_B_ti, t9a);
931 }
932
933 if (p1 > 0)
934 {
935 f64 t2a = bt1data[index1[0]].len * t1b[0];
936 f64 t2b = bt1data[index1[0]].len * t1a[0] / bt1data[index1[1]].len;
937
938 f64 t3a = bt1data[index1[3]].len * t1b[1];
939 f64 t3b = bt1data[index1[3]].len * t1a[1] / bt1data[index1[2]].len;
940
941 f64 t4c[3]; f64 t5c[3]; f64 t6a[3]; f64 t7a[3];
942 const i32s cp[3][3] = { { 0, 1, 2 }, { 1, 2, 0 }, { 2, 0, 1 } };
943
944 for (i32s n2 = 0;n2 < 3;n2++)
945 {
946 f64 t4a = bt1data[index1[0]].dlen[dir1[0]][cp[n2][1]] * bt1data[index1[1]].dlen[dir1[1]][cp[n2][2]];
947 f64 t4b = bt1data[index1[0]].dlen[dir1[0]][cp[n2][2]] * bt1data[index1[1]].dlen[dir1[1]][cp[n2][1]];
948 t4c[n2] = (t4a - t4b) / t2a;
949
950 f64 t5a = bt1data[index1[2]].dlen[dir1[2]][cp[n2][2]] * bt1data[index1[3]].dlen[dir1[3]][cp[n2][1]];
951 f64 t5b = bt1data[index1[2]].dlen[dir1[2]][cp[n2][1]] * bt1data[index1[3]].dlen[dir1[3]][cp[n2][2]];
952 t5c[n2] = (t5a - t5b) / t3a;
953
954 d1[l2g_mm[atmi[0]] * 3 + n2] += t4c[n2] * t9b;
955 d1[l2g_mm[atmi[3]] * 3 + n2] += t5c[n2] * t9b;
956
957 t6a[n2] = (t2b - 1.0) * t4c[n2] - t3b * t5c[n2];
958 t7a[n2] = (t3b - 1.0) * t5c[n2] - t2b * t4c[n2];
959
960 d1[l2g_mm[atmi[1]] * 3 + n2] += t6a[n2] * t9b;
961 d1[l2g_mm[atmi[2]] * 3 + n2] += t7a[n2] * t9b;
962 }
963 }
964 }
965 }
966
ComputeBT4(i32u p1)967 void eng1_mm_default_bt::ComputeBT4(i32u p1)
968 {
969 energy_bt4 = 0.0;
970
971 atom ** atmtab = GetSetup()->GetMMAtoms();
972
973 for (i32s n1 = 0;n1 < (i32s) bt4_vector.size();n1++)
974 {
975 i32s * atmi = bt4_vector[n1].atmi;
976
977 i32s index2 = bt4_vector[n1].index2;
978 bool dir2 = bt4_vector[n1].dir2;
979 i32s * index1 = bt4_vector[n1].index1;
980 bool * dir1 = bt4_vector[n1].dir1;
981
982 f64 t1a[3];
983 t1a[0] = bt1data[index1[0]].dlen[dir1[0]][1] * bt1data[index1[1]].dlen[dir1[1]][2] -
984 bt1data[index1[0]].dlen[dir1[0]][2] * bt1data[index1[1]].dlen[dir1[1]][1];
985 t1a[1] = bt1data[index1[0]].dlen[dir1[0]][2] * bt1data[index1[1]].dlen[dir1[1]][0] -
986 bt1data[index1[0]].dlen[dir1[0]][0] * bt1data[index1[1]].dlen[dir1[1]][2];
987 t1a[2] = bt1data[index1[0]].dlen[dir1[0]][0] * bt1data[index1[1]].dlen[dir1[1]][1] -
988 bt1data[index1[0]].dlen[dir1[0]][1] * bt1data[index1[1]].dlen[dir1[1]][0];
989
990 f64 t1b = 0.0;
991 for (i32s n2 = 0;n2 < 3;n2++)
992 {
993 t1b += bt1data[index1[2]].dlen[dir1[2]][n2] * t1a[n2];
994 }
995
996 f64 t1c = 1.0 - bt2data[index2].csa * bt2data[index2].csa;
997 f64 t1d = sqrt(t1c);
998
999 f64 t1e = t1b / t1d;
1000
1001 //cout << "t1e = " << t1e << " " << t1e * t1e << endl;
1002 //v3d<f64> v1(crd[l2g_mm[atmi[1]]], crd[l2g_mm[atmi[0]]]);
1003 //v3d<f64> v2(crd[l2g_mm[atmi[1]]], crd[l2g_mm[atmi[2]]]);
1004 //v3d<f64> v3(crd[l2g_mm[atmi[1]]], crd[l2g_mm[atmi[3]]]);
1005 //v3d<f64> v4 = v1.vpr(v2);
1006 //f64 ang = v3.ang(v4);
1007 //cout << "check = " << cos(ang) << " " << cos(ang) * cos(ang) << endl;
1008 //int zzz; cin >> zzz;
1009
1010 if (t1e < -1.0) t1e = -1.0; // domain check...
1011 if (t1e > +1.0) t1e = +1.0; // domain check...
1012
1013 // f = a(asin(x)-b)^2
1014 // df/dx = 2a(asin(x)-b)/sqrt(1-x^2)
1015
1016 f64 t1f = asin(t1e) - bt4_vector[n1].opt;
1017 f64 t2a = bt4_vector[n1].fc * t1f * t1f;
1018
1019 energy_bt4 += t2a;
1020
1021 if (ECOMPstore != NULL)
1022 {
1023 vector<int> gv;
1024 gv.push_back(atmtab[atmi[0]]->ecomp_grp_i);
1025 gv.push_back(atmtab[atmi[1]]->ecomp_grp_i);
1026 gv.push_back(atmtab[atmi[2]]->ecomp_grp_i);
1027 gv.push_back(atmtab[atmi[3]]->ecomp_grp_i);
1028
1029 ecomp_AddStoreX(gv, ECOMP_DATA_IND_B_ti, t2a);
1030 }
1031
1032 if (p1 > 0)
1033 {
1034 f64 t1g = 2.0 * bt4_vector[n1].fc * t1f / sqrt(1.0 - t1e * t1e);
1035 // f64 t1h = bt2data[index2].csa / t1d;
1036
1037 f64 t9b = 1.0 - bt2data[index2].csa * bt2data[index2].csa;
1038 f64 t9c = sqrt(t9b);
1039
1040 for (i32s n2 = 0;n2 < 3;n2++)
1041 {
1042 // this code is borrowed from eng1_sf.cpp file...
1043 // might not be optimal, just a quick way to go ahead.
1044
1045 f64 t6a[2]; // epsilon i
1046 t6a[0] = bt2data[index2].dcsa[dir2 ? 0 : 2][n2] * bt2data[index2].csa / t9b;
1047 t6a[1] = bt2data[index2].dcsa[dir2 ? 2 : 0][n2] * bt2data[index2].csa / t9b;
1048
1049 f64 t6b[2]; // sigma ii / r2X
1050 t6b[0] = (1.0 - bt1data[index1[0]].dlen[dir1[0]][n2] * bt1data[index1[0]].dlen[dir1[0]][n2]) / bt1data[index1[0]].len;
1051 t6b[1] = (1.0 - bt1data[index1[1]].dlen[dir1[1]][n2] * bt1data[index1[1]].dlen[dir1[1]][n2]) / bt1data[index1[1]].len;
1052
1053 i32s n3[2];
1054 n3[0] = (n2 + 1) % 3; // index j
1055 n3[1] = (n2 + 2) % 3; // index k
1056
1057 f64 t6c[2]; // sigma ij / r2X
1058 t6c[0] = -bt1data[index1[0]].dlen[dir1[0]][n2] * bt1data[index1[0]].dlen[dir1[0]][n3[0]] / bt1data[index1[0]].len;
1059 t6c[1] = -bt1data[index1[1]].dlen[dir1[1]][n2] * bt1data[index1[1]].dlen[dir1[1]][n3[0]] / bt1data[index1[1]].len;
1060
1061 f64 t6d[2]; // sigma ik / r2X
1062 t6d[0] = -bt1data[index1[0]].dlen[dir1[0]][n2] * bt1data[index1[0]].dlen[dir1[0]][n3[1]] / bt1data[index1[0]].len;
1063 t6d[1] = -bt1data[index1[1]].dlen[dir1[1]][n2] * bt1data[index1[1]].dlen[dir1[1]][n3[1]] / bt1data[index1[1]].len;
1064
1065 f64 t5a[2][3]; // components of dc/di
1066 t5a[0][n2] = (t6c[0] * bt1data[index1[1]].dlen[dir1[1]][n3[1]] -
1067 t6d[0] * bt1data[index1[1]].dlen[dir1[1]][n3[0]] + t1a[n2] * t6a[0]) / t9c;
1068 t5a[0][n3[0]] = (t6d[0] * bt1data[index1[1]].dlen[dir1[1]][n2] -
1069 t6b[0] * bt1data[index1[1]].dlen[dir1[1]][n3[1]] + t1a[n3[0]] * t6a[0]) / t9c;
1070 t5a[0][n3[1]] = (t6b[0] * bt1data[index1[1]].dlen[dir1[1]][n3[0]] -
1071 t6c[0] * bt1data[index1[1]].dlen[dir1[1]][n2] + t1a[n3[1]] * t6a[0]) / t9c;
1072 t5a[1][n2] = (t6d[1] * bt1data[index1[0]].dlen[dir1[0]][n3[0]] -
1073 t6c[1] * bt1data[index1[0]].dlen[dir1[0]][n3[1]] + t1a[n2] * t6a[1]) / t9c;
1074 t5a[1][n3[0]] = (t6b[1] * bt1data[index1[0]].dlen[dir1[0]][n3[1]] -
1075 t6d[1] * bt1data[index1[0]].dlen[dir1[0]][n2] + t1a[n3[0]] * t6a[1]) / t9c;
1076 t5a[1][n3[1]] = (t6c[1] * bt1data[index1[0]].dlen[dir1[0]][n2] -
1077 t6b[1] * bt1data[index1[0]].dlen[dir1[0]][n3[0]] + t1a[n3[1]] * t6a[1]) / t9c;
1078
1079 f64 tmp1a = t5a[0][0] * bt1data[index1[2]].dlen[dir1[2]][0] +
1080 t5a[0][1] * bt1data[index1[2]].dlen[dir1[2]][1] +
1081 t5a[0][2] * bt1data[index1[2]].dlen[dir1[2]][2];
1082
1083 f64 tmp3a = t5a[1][0] * bt1data[index1[2]].dlen[dir1[2]][0] +
1084 t5a[1][1] * bt1data[index1[2]].dlen[dir1[2]][1] +
1085 t5a[1][2] * bt1data[index1[2]].dlen[dir1[2]][2];
1086
1087 f64 tmp4a = 0.0;
1088 for (i32s n4 = 0;n4 < 3;n4++)
1089 {
1090 f64 tmp4b;
1091 if (n2 != n4) tmp4b = -bt1data[index1[2]].dlen[!dir1[2]][n2] * bt1data[index1[2]].dlen[!dir1[2]][n4];
1092 else tmp4b = 1.0 - bt1data[index1[2]].dlen[!dir1[2]][n4] * bt1data[index1[2]].dlen[!dir1[2]][n4];
1093 tmp4a += (tmp4b / bt1data[index1[2]].len) * (t1a[n4] / t1d);
1094 }
1095
1096 d1[l2g_mm[atmi[0]] * 3 + n2] += tmp1a * t1g;
1097 d1[l2g_mm[atmi[1]] * 3 + n2] -= (tmp1a + tmp3a + tmp4a) * t1g;
1098 d1[l2g_mm[atmi[2]] * 3 + n2] += tmp3a * t1g;
1099 d1[l2g_mm[atmi[3]] * 3 + n2] += tmp4a * t1g;
1100 }
1101 }
1102 }
1103 }
1104
1105 /*################################################################################################*/
1106
eng1_mm_default_nbt_bp(setup * p1,i32u p2)1107 eng1_mm_default_nbt_bp::eng1_mm_default_nbt_bp(setup * p1, i32u p2) :
1108 engine(p1, p2),
1109 eng1_mm(p1, p2),
1110 engine_bp(p1, p2)
1111 {
1112 atom ** atmtab = GetSetup()->GetMMAtoms();
1113 // bond ** bndtab = GetSetup()->GetMMBonds();
1114
1115 ostream * ostr = NULL; // do not print output.
1116 // ostream * ostr = & cout; // print output to cout.
1117
1118 // also check engine::bp_center!!!
1119 // also check engine::bp_center!!!
1120 // also check engine::bp_center!!!
1121 bp_fc_solute = 5000.0; // 50 kJ/(mol*�^2) = 5000 kJ/(mol*(nm)^2)
1122 bp_fc_solvent = 12500.0; // 125 kJ/(mol*�^2) = 12500 kJ/(mol*(nm)^2)
1123
1124 if (ostr != NULL && use_bp)
1125 {
1126 (* ostr) << _("use_bp ; ");
1127 (* ostr) << bp_rad_solute << " " << bp_fc_solute << " ; ";
1128 (* ostr) << bp_rad_solvent << " " << bp_fc_solvent << endl;
1129 }
1130
1131 if (ostr != NULL) (* ostr) << _("creating nbt1-terms: ");
1132 i32s nbt1_err = 0;
1133
1134 for (i32s ind1 = 0;ind1 < GetSetup()->GetMMAtomCount() - 1;ind1++)
1135 {
1136 for (i32s ind2 = ind1 + 1;ind2 < GetSetup()->GetMMAtomCount();ind2++)
1137 {
1138 i32s test = range_cr1[ind1];
1139 while (test < range_cr1[ind1 + 1])
1140 {
1141 if (cr1[test] == atmtab[ind2]) break;
1142 else test++;
1143 }
1144
1145 // if this is true, then the atoms are not 1-2 or 1-3 related.
1146 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1147 if (test == range_cr1[ind1 + 1])
1148 {
1149 test = range_cr2[ind1];
1150 while (test < range_cr2[ind1 + 1])
1151 {
1152 if (cr2[test] == atmtab[ind2]) break;
1153 else test++;
1154 }
1155
1156 // if this is true, then the atoms are 1-4 related.
1157 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1158 const bool is14 = (test != range_cr2[ind1 + 1]);
1159 // if (is14) cout << "DEBUG ; is 1-4 : " << ind1 << " " << ind2 << endl;
1160
1161 mm_default_nbt1 newnbt1;
1162 newnbt1.atmi[0] = ind1; // this is a local index...
1163 newnbt1.atmi[1] = ind2; // this is a local index...
1164
1165 bool success = false;
1166 if (dynamic_cast<setup1_mm *>(GetSetup())->GetExceptions())
1167 {
1168 success = default_tables::GetInstance()->e_Init(this, & newnbt1, is14);
1169 }
1170
1171 if (success != true)
1172 {
1173
1174 // see also eng1_mm_default_nbt_mim::UpdateTerms() ; should be the same!!!
1175 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1176
1177 bool errors = false;
1178 const default_at * at;
1179
1180 f64 r1 = 0.150; f64 e1 = 0.175; // default...
1181 at = default_tables::GetInstance()->GetAtomType(atmtab[ind1]->atmtp);
1182 if (at != NULL) { r1 = at->vdw_R; e1 = at->vdw_E; }
1183 else errors = true;
1184
1185 f64 r2 = 0.150; f64 e2 = 0.175; // default...
1186 at = default_tables::GetInstance()->GetAtomType(atmtab[ind2]->atmtp);
1187 if (at != NULL) { r2 = at->vdw_R; e2 = at->vdw_E; }
1188 else errors = true;
1189
1190 f64 optdist = r1 + r2;
1191 f64 energy = sqrt(e1 * e2);
1192
1193 f64 charge1 = atmtab[ind1]->charge;
1194 f64 charge2 = atmtab[ind2]->charge;
1195 newnbt1.qq = 138.9354518 * charge1 * charge2;
1196
1197 if (is14)
1198 {
1199 energy *= 0.5;
1200 newnbt1.qq *= 0.75;
1201 }
1202
1203 f64 tmp1 = optdist * pow(1.0 * energy, 1.0 / 12.0);
1204 f64 tmp2 = optdist * pow(2.0 * energy, 1.0 / 6.0);
1205
1206 newnbt1.kr = tmp1;
1207 newnbt1.kd = tmp2;
1208
1209 if (!errors) success = true;
1210 }
1211
1212 nbt1_err += !success;
1213 nbt1_vector.push_back(newnbt1);
1214 }
1215 }
1216 }
1217
1218 // report possible errors...
1219
1220 i32s total_err = nbt1_err;
1221 if (total_err && GetSetup()->GetModel()->verbosity >= 2)
1222 {
1223 ostringstream str;
1224 str << _("WARNING : there were ") << total_err << _(" missing parameters in the nonbonded terms.") << endl << ends;
1225 GetSetup()->GetModel()->PrintToLog(str.str().c_str());
1226 }
1227 }
1228
~eng1_mm_default_nbt_bp(void)1229 eng1_mm_default_nbt_bp::~eng1_mm_default_nbt_bp(void)
1230 {
1231 }
1232
ComputeNBT1(i32u p1)1233 void eng1_mm_default_nbt_bp::ComputeNBT1(i32u p1)
1234 {
1235 energy_nbt1a = 0.0;
1236 energy_nbt1b = 0.0;
1237 energy_nbt1c = 0.0;
1238 energy_nbt1d = 0.0;
1239
1240 atom ** atmtab = GetSetup()->GetMMAtoms();
1241
1242 // an additional pass for the boundary potential (optional).
1243 // an additional pass for the boundary potential (optional).
1244 // an additional pass for the boundary potential (optional).
1245
1246 if (use_bp)
1247 {
1248 if (nd_eval != NULL) nd_eval->AddCycle();
1249 for (i32s n1 = 0;n1 < GetSetup()->GetMMAtomCount();n1++)
1250 {
1251 f64 radius = bp_rad_solute;
1252 f64 fc = bp_fc_solute;
1253
1254 if (atmtab[n1]->flags & ATOMFLAG_IS_SOLVENT_ATOM)
1255 {
1256 radius = bp_rad_solvent;
1257 fc = bp_fc_solvent;
1258 }
1259
1260 f64 t1a[3]; f64 t1b = 0.0;
1261 for (i32s n2 = 0;n2 < 3;n2++)
1262 {
1263 f64 t9a = 0.0;//bp_center[n2];
1264 f64 t9b = crd[l2g_mm[n1] * 3 + n2];
1265
1266 t1a[n2] = t9a - t9b;
1267 t1b += t1a[n2] * t1a[n2];
1268 }
1269
1270 f64 t1c = sqrt(t1b);
1271
1272 if (nd_eval != NULL && (atmtab[n1]->flags & ATOMFLAG_MEASURE_ND_RDF)) nd_eval->AddValue(t1c);
1273
1274 if (rdf_eval != NULL && rdf_eval->count_begin > -0.5)
1275 {
1276 bool is_in_rdf_counting_window = true;
1277 if (t1c < rdf_eval->count_begin) is_in_rdf_counting_window = false;
1278 if (t1c >= rdf_eval->count_end) is_in_rdf_counting_window = false;
1279
1280 if (is_in_rdf_counting_window) atmtab[n1]->flags |= ATOMFLAG_COUNT_IN_RDF;
1281 else atmtab[n1]->flags &= (~ATOMFLAG_COUNT_IN_RDF);
1282 }
1283
1284 if (t1c < radius) continue;
1285
1286 // f = a(x-b)^2
1287 // df/dx = 2a(x-b)
1288
1289 f64 t2a = t1c - radius;
1290 f64 t2b = fc * t2a * t2a;
1291
1292 energy_bt1 += t2b;
1293
1294 // do not add to ECOMP!
1295
1296 if (p1 > 0)
1297 {
1298 f64 t2c = 2.0 * fc * t2a;
1299
1300 for (i32s n2 = 0;n2 < 3;n2++)
1301 {
1302 f64 t2d = (t1a[n2] / t1c) * t2c;
1303
1304 d1[l2g_mm[n1] * 3 + n2] -= t2d;
1305 }
1306 }
1307 }
1308 }
1309
1310 // the nonbonded terms begin...
1311 // the nonbonded terms begin...
1312 // the nonbonded terms begin...
1313
1314 if (rdf_eval != NULL) rdf_eval->AddCycle();
1315 for (i32s n1 = 0;n1 < (i32s) nbt1_vector.size();n1++)
1316 {
1317 i32s * atmi = nbt1_vector[n1].atmi;
1318
1319 f64 t1a[3]; f64 t1b = 0.0;
1320 for (i32s n2 = 0;n2 < 3;n2++)
1321 {
1322 f64 t2a = crd[l2g_mm[atmi[0]] * 3 + n2];
1323 f64 t2b = crd[l2g_mm[atmi[1]] * 3 + n2];
1324
1325 t1a[n2] = t2a - t2b;
1326 t1b += t1a[n2] * t1a[n2];
1327 }
1328
1329 f64 t1c = sqrt(t1b);
1330 if (rdf_eval != NULL)
1331 {
1332 bool flag = true;
1333 if (!(atmtab[atmi[0]]->flags & ATOMFLAG_MEASURE_ND_RDF)) flag = false;
1334 if (!(atmtab[atmi[1]]->flags & ATOMFLAG_MEASURE_ND_RDF)) flag = false;
1335
1336 if (rdf_eval->count_begin > -0.5) // check the counting window, if needed...
1337 {
1338 if (!(atmtab[atmi[0]]->flags & ATOMFLAG_COUNT_IN_RDF)) flag = false;
1339 if (!(atmtab[atmi[1]]->flags & ATOMFLAG_COUNT_IN_RDF)) flag = false;
1340 }
1341
1342 if (flag) rdf_eval->AddValue(t1c);
1343 }
1344
1345 // f1 = (r/a)^-12 - (r/b)^-6
1346 // df1/dr = -12/a(r/a)^-13 + 6/b(r/b)^-7
1347
1348 f64 t3a = t1c / nbt1_vector[n1].kr;
1349 f64 t3b = t1c / nbt1_vector[n1].kd;
1350
1351 f64 t4a = t3a * t3a * t3a; f64 t4b = t4a * t4a; f64 t4c = t4b * t4b; // ^3 ^6 ^12
1352 f64 t5a = t3b * t3b * t3b; f64 t5b = t5a * t5a; // ^3 ^6
1353
1354 f64 t6a = 1.0 / (t4c) - 1.0 / (t5b);
1355 energy_nbt1a += t6a;
1356
1357 // f2 = Q/r
1358 // df2/dr = -Q/r^2
1359
1360 f64 t6b = nbt1_vector[n1].qq / t1c;
1361 energy_nbt1b += t6b;
1362
1363 // f64 tote = t6a + t6b;
1364
1365 if (ECOMPstore != NULL)
1366 {
1367 const int iA = atmtab[atmi[0]]->ecomp_grp_i;
1368 const int iB = atmtab[atmi[1]]->ecomp_grp_i;
1369
1370 ecomp_AddStore2(iA, iB, ECOMP_DATA_IND_NB_lj, t6a);
1371 ecomp_AddStore2(iA, iB, ECOMP_DATA_IND_NB_es, t6b);
1372 }
1373
1374 if (p1 > 0)
1375 {
1376 f64 t7a = 12.0 / (nbt1_vector[n1].kr * t4c * t3a);
1377 f64 t7b = 6.0 / (nbt1_vector[n1].kd * t5b * t3b);
1378
1379 f64 t8a = nbt1_vector[n1].qq / t1b;
1380
1381 f64 t9a = t7b - t7a - t8a;
1382 for (i32s n2 = 0;n2 < 3;n2++)
1383 {
1384 f64 t9b = (t1a[n2] / t1c) * t9a;
1385
1386 d1[l2g_mm[atmi[0]] * 3 + n2] += t9b;
1387 d1[l2g_mm[atmi[1]] * 3 + n2] -= t9b;
1388 }
1389 }
1390 }
1391 }
1392
1393 /*################################################################################################*/
1394
eng1_mm_default_nbt_mim(setup * p1,i32u p2)1395 eng1_mm_default_nbt_mim::eng1_mm_default_nbt_mim(setup * p1, i32u p2) :
1396 engine(p1, p2),
1397 eng1_mm(p1, p2),
1398 engine_pbc(p1, p2)
1399 {
1400 fGL mindim = box_HALFdim[0];
1401 if (box_HALFdim[1] < mindim) mindim = box_HALFdim[1];
1402 if (box_HALFdim[2] < mindim) mindim = box_HALFdim[2];
1403
1404 sw1 = 0.6; if (sw1 < (mindim - 0.4)) sw1 = mindim - 0.4; // will trigger if boxdim < 2.0 nm!!!
1405 sw2 = shft1 = mindim - 0.2;
1406
1407 limit = mindim;
1408
1409 // calculate the actual values...
1410
1411 sw1 = sw1 * sw1;
1412 sw2 = sw2 * sw2;
1413
1414 swA = 3.0 * sw1;
1415 swB = pow(sw2 - sw1, 3.0);
1416
1417 shft3 = pow(shft1, 3.0);
1418
1419 limit = limit * limit;
1420
1421 nbt1_vector.reserve(250000);
1422
1423 RequestNeighborListUpdate();
1424 }
1425
~eng1_mm_default_nbt_mim(void)1426 eng1_mm_default_nbt_mim::~eng1_mm_default_nbt_mim(void)
1427 {
1428 }
1429
ComputeNBT1(i32u p1)1430 void eng1_mm_default_nbt_mim::ComputeNBT1(i32u p1)
1431 {
1432 energy_nbt1a = 0.0;
1433 energy_nbt1b = 0.0;
1434 energy_nbt1c = 0.0;
1435 energy_nbt1d = 0.0;
1436
1437 atom ** atmtab = GetSetup()->GetMMAtoms();
1438
1439 if (update_neighbor_list) UpdateTerms();
1440
1441 // the nonbonded terms begin...
1442 // the nonbonded terms begin...
1443 // the nonbonded terms begin...
1444
1445 for (i32s n1 = 0;n1 < (i32s) nbt1_vector.size();n1++)
1446 {
1447 i32s * atmi = nbt1_vector[n1].atmi;
1448
1449 f64 t1a[3]; f64 t1b = 0.0;
1450 for (i32s n2 = 0;n2 < 3;n2++)
1451 {
1452 f64 t2a = crd[l2g_mm[atmi[0]] * 3 + n2];
1453 // ^^^^^ either at primary cell OR at 1st neighbor cell (1 of 26).
1454 if (t2a < -box_HALFdim[n2])
1455 {
1456 t2a += 2.0 * box_HALFdim[n2];
1457 if (t2a < -box_HALFdim[n2]) assertion_failed(__FILE__, __LINE__, "PBC failed ; a-");
1458 }
1459 else if (t2a > +box_HALFdim[n2])
1460 {
1461 t2a -= 2.0 * box_HALFdim[n2];
1462 if (t2a > +box_HALFdim[n2]) assertion_failed(__FILE__, __LINE__, "PBC failed ; a+");
1463 }
1464
1465 f64 t2b = crd[l2g_mm[atmi[1]] * 3 + n2];
1466 // ^^^^^ either at primary cell OR at 1st neighbor cell (1 of 26).
1467 if (t2b < -box_HALFdim[n2])
1468 {
1469 t2b += 2.0 * box_HALFdim[n2];
1470 if (t2b < -box_HALFdim[n2]) assertion_failed(__FILE__, __LINE__, "PBC failed ; b-");
1471 }
1472 else if (t2b > +box_HALFdim[n2])
1473 {
1474 t2b -= 2.0 * box_HALFdim[n2];
1475 if (t2b > +box_HALFdim[n2]) assertion_failed(__FILE__, __LINE__, "PBC failed ; b+");
1476 }
1477
1478 t1a[n2] = t2a - t2b;
1479
1480 if (t1a[n2] < -box_HALFdim[n2])
1481 {
1482 t1a[n2] += 2.0 * box_HALFdim[n2];
1483 }
1484 else if (t1a[n2] > +box_HALFdim[n2])
1485 {
1486 t1a[n2] -= 2.0 * box_HALFdim[n2];
1487 }
1488
1489 t1b += t1a[n2] * t1a[n2];
1490 }
1491
1492 f64 t1c = sqrt(t1b);
1493
1494 // f1 = (r/a)^-12 - (r/b)^-6
1495 // df1/dr = -12/a(r/a)^-13 + 6/b(r/b)^-7
1496
1497 f64 t3a = t1c / nbt1_vector[n1].kr;
1498 f64 t3b = t1c / nbt1_vector[n1].kd;
1499
1500 f64 t4a = t3a * t3a * t3a; f64 t4b = t4a * t4a; f64 t4c = t4b * t4b; // ^3 ^6 ^12
1501 f64 t5a = t3b * t3b * t3b; f64 t5b = t5a * t5a; // ^3 ^6
1502
1503 f64 t6a = 1.0 / (t4c) - 1.0 / (t5b);
1504
1505 // s1 = (rE^2 - r^2)^2 * (rE^2 + 2r^2 - 3rB^2) / (rE^2 - rB^2)^3
1506 // ds1/dr = [this will yield 2 terms quite easily...]
1507
1508 f64 t3x; // value
1509 f64 t3y; f64 t3z; // derivative
1510 if (t1b < sw1)
1511 {
1512 t3x = 1.0;
1513 t3y = t3z = 0.0;
1514 }
1515 else if (t1b > sw2)
1516 {
1517 t3x = 0.0;
1518 t3y = t3z = 0.0;
1519 }
1520 else
1521 {
1522 f64 t3c = sw2 - t1b; f64 t3d = t3c * t3c;
1523 f64 t3e = sw2 + 2.0 * t1b - swA;
1524
1525 t3x = t3d * t3e / swB;
1526 t3y = 4.0 * t1c * t3d / swB;
1527 t3z = 4.0 * t1c * t3c * t3e / swB;
1528 }
1529
1530 energy_nbt1a += t6a * t3x;
1531
1532 // f2 = Q/r
1533 // df2/dr = -Q/r^2
1534
1535 f64 t6b = nbt1_vector[n1].qq / t1c;
1536
1537 // s2 = (1 - (r/rE)^3)^2
1538 // ds2/dr = -6r^2 * (1 - (r/rE)^3) / rE^3
1539
1540 f64 t4x; // value
1541 f64 t4y; // derivative
1542 if (t1c > shft1)
1543 {
1544 t4x = 0.0;
1545 t4y = 0.0;
1546 }
1547 else
1548 {
1549 f64 t4d = t1b * t1c / shft3;
1550 f64 t4e = 1.0 - t4d;
1551
1552 t4x = t4e * t4e;
1553 t4y = 6.0 * t1b * t4e / shft3;
1554 }
1555
1556 energy_nbt1b += t6b * t4x;
1557
1558 // f64 tote = t6a * t3x + t6b * t4x;
1559
1560 if (ECOMPstore != NULL)
1561 {
1562 const int iA = atmtab[atmi[0]]->ecomp_grp_i;
1563 const int iB = atmtab[atmi[1]]->ecomp_grp_i;
1564
1565 ecomp_AddStore2(iA, iB, ECOMP_DATA_IND_NB_lj, (t6a * t3x));
1566 ecomp_AddStore2(iA, iB, ECOMP_DATA_IND_NB_es, (t6b * t4x));
1567 }
1568
1569 if (p1 > 0)
1570 {
1571 f64 t7a = 12.0 / (nbt1_vector[n1].kr * t4c * t3a);
1572 f64 t7b = 6.0 / (nbt1_vector[n1].kd * t5b * t3b);
1573
1574 f64 t8a = nbt1_vector[n1].qq / t1b;
1575
1576 f64 t9a = (t7b - t7a) * t3x + t6a * (t3y - t3z);
1577 f64 t9b = t8a * t4x + t6b * t4y;
1578
1579 f64 t9c = t9a - t9b;
1580
1581 for (i32s n2 = 0;n2 < 3;n2++)
1582 {
1583 f64 t9d = (t1a[n2] / t1c) * t9c;
1584
1585 d1[l2g_mm[atmi[0]] * 3 + n2] += t9d;
1586 d1[l2g_mm[atmi[1]] * 3 + n2] -= t9d;
1587
1588 if (do_virial)
1589 {
1590 //f64 crdA = crd[l2g_mm[atmi[0]] * 3 + n2]; eik�y...
1591 //f64 crdB = crd[l2g_mm[atmi[1]] * 3 + n2]; eik�y...
1592 //f64 rAB = crdA - crdB; eik�y...
1593 //////////////////////////////////////////////////////////
1594 f64 rAB = t1a[n2];
1595
1596 virial[n2] -= rAB * t9d; // F = -dE/dr
1597 }
1598 }
1599 }
1600 }
1601 }
1602
UpdateTerms(void)1603 void eng1_mm_default_nbt_mim::UpdateTerms(void)
1604 {
1605 update_neighbor_list = false; // processed...
1606
1607 atom ** atmtab = GetSetup()->GetMMAtoms();
1608 // bond ** bndtab = GetSetup()->GetMMBonds();
1609
1610 ostream * ostr = NULL; // do not print output.
1611 // ostream * ostr = & cout; // print output to cout.
1612
1613 nbt1_vector.resize(0);
1614
1615 if (ostr != NULL) (* ostr) << _("creating nbt1-terms: ");
1616 i32s nbt1_err = 0;
1617
1618 for (i32s ind1 = 0;ind1 < GetSetup()->GetMMAtomCount() - 1;ind1++)
1619 {
1620 for (i32s ind2 = ind1 + 1;ind2 < GetSetup()->GetMMAtomCount();ind2++)
1621 {
1622 i32s test = range_cr1[ind1];
1623 while (test < range_cr1[ind1 + 1])
1624 {
1625 if (cr1[test] == atmtab[ind2]) break;
1626 else test++;
1627 }
1628
1629 // if this is true, then the atoms are not 1-2 or 1-3 related.
1630 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1631 if (test == range_cr1[ind1 + 1])
1632 {
1633 f64 t1a; f64 t1b = 0.0;
1634 for (i32s n1 = 0;n1 < 3;n1++)
1635 {
1636 f64 t2a = crd[l2g_mm[ind1] * 3 + n1];
1637 // ^^^^^ either at primary cell OR at 1st neighbor cell (1 of 26).
1638 if (t2a < -box_HALFdim[n1])
1639 {
1640 t2a += 2.0 * box_HALFdim[n1];
1641 if (t2a < -box_HALFdim[n1]) assertion_failed(__FILE__, __LINE__, "PBC failed ; a-");
1642 }
1643 else if (t2a > +box_HALFdim[n1])
1644 {
1645 t2a -= 2.0 * box_HALFdim[n1];
1646 if (t2a > +box_HALFdim[n1]) assertion_failed(__FILE__, __LINE__, "PBC failed ; a+");
1647 }
1648
1649 f64 t2b = crd[l2g_mm[ind2] * 3 + n1];
1650 // ^^^^^ either at primary cell OR at 1st neighbor cell (1 of 26).
1651 if (t2b < -box_HALFdim[n1])
1652 {
1653 t2b += 2.0 * box_HALFdim[n1];
1654 if (t2b < -box_HALFdim[n1]) assertion_failed(__FILE__, __LINE__, "PBC failed ; b-");
1655 }
1656 else if (t2b > +box_HALFdim[n1])
1657 {
1658 t2b -= 2.0 * box_HALFdim[n1];
1659 if (t2b > +box_HALFdim[n1]) assertion_failed(__FILE__, __LINE__, "PBC failed ; b+");
1660 }
1661
1662 t1a = t2a - t2b;
1663
1664 if (t1a < -box_HALFdim[n1])
1665 {
1666 t1a += 2.0 * box_HALFdim[n1];
1667 }
1668 else if (t1a > +box_HALFdim[n1])
1669 {
1670 t1a -= 2.0 * box_HALFdim[n1];
1671 }
1672
1673 t1b += t1a * t1a;
1674 }
1675
1676 if (t1b > limit) continue;
1677
1678 test = range_cr2[ind1];
1679 while (test < range_cr2[ind1 + 1])
1680 {
1681 if (cr2[test] == atmtab[ind2]) break;
1682 else test++;
1683 }
1684
1685 // if this is true, then the atoms are 1-4 related.
1686 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1687 const bool is14 = (test != range_cr2[ind1 + 1]);
1688 // if (is14) cout << "DEBUG ; is 1-4 : " << ind1 << " " << ind2 << endl;
1689
1690 mm_default_nbt1 newnbt1;
1691 newnbt1.atmi[0] = ind1;
1692 newnbt1.atmi[1] = ind2;
1693
1694 bool success = false;
1695 if (dynamic_cast<setup1_mm *>(GetSetup())->GetExceptions())
1696 {
1697 success = default_tables::GetInstance()->e_Init(this, & newnbt1, is14);
1698 }
1699
1700 if (success != true)
1701 {
1702
1703 // see also eng1_mm_default_nbt_bp ctor ; should be the same!!!
1704 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1705
1706 bool errors = false;
1707 const default_at * at;
1708
1709 f64 r1 = 0.150; f64 e1 = 0.175; // default...
1710 at = default_tables::GetInstance()->GetAtomType(atmtab[ind1]->atmtp);
1711 if (at != NULL) { r1 = at->vdw_R; e1 = at->vdw_E; }
1712 else errors = true;
1713
1714 f64 r2 = 0.150; f64 e2 = 0.175; // default...
1715 at = default_tables::GetInstance()->GetAtomType(atmtab[ind2]->atmtp);
1716 if (at != NULL) { r2 = at->vdw_R; e2 = at->vdw_E; }
1717 else errors = true;
1718
1719 f64 optdist = r1 + r2;
1720 f64 energy = sqrt(e1 * e2);
1721
1722 f64 charge1 = atmtab[ind1]->charge;
1723 f64 charge2 = atmtab[ind2]->charge;
1724 newnbt1.qq = 138.9354518 * charge1 * charge2;
1725
1726 if (is14)
1727 {
1728 energy *= 0.5;
1729 newnbt1.qq *= 0.75;
1730 }
1731
1732 f64 tmp1 = optdist * pow(1.0 * energy, 1.0 / 12.0);
1733 f64 tmp2 = optdist * pow(2.0 * energy, 1.0 / 6.0);
1734
1735 newnbt1.kr = tmp1;
1736 newnbt1.kd = tmp2;
1737
1738 if (!errors) success = true;
1739 }
1740
1741 nbt1_err += !success;
1742 nbt1_vector.push_back(newnbt1);
1743 }
1744 }
1745 }
1746
1747 // report possible errors...
1748
1749 i32s total_err = nbt1_err;
1750 if (total_err && GetSetup()->GetModel()->verbosity >= 2)
1751 {
1752 ostringstream str;
1753 str << _("WARNING : there were ") << total_err << _(" missing parameters in the nonbonded terms.") << endl << ends;
1754 GetSetup()->GetModel()->PrintToLog(str.str().c_str());
1755 }
1756 }
1757
1758 /*################################################################################################*/
1759
eng1_mm_default_bp(setup * p1,i32u p2)1760 eng1_mm_default_bp::eng1_mm_default_bp(setup * p1, i32u p2) :
1761 engine(p1, p2),
1762 eng1_mm(p1, p2),
1763 engine_bp(p1, p2),
1764 eng1_mm_default_bt(p1, p2),
1765 eng1_mm_default_nbt_bp(p1, p2)
1766 {
1767 }
1768
~eng1_mm_default_bp(void)1769 eng1_mm_default_bp::~eng1_mm_default_bp(void)
1770 {
1771 }
1772
1773 /*################################################################################################*/
1774
eng1_mm_default_mim(setup * p1,i32u p2)1775 eng1_mm_default_mim::eng1_mm_default_mim(setup * p1, i32u p2) :
1776 engine(p1, p2),
1777 eng1_mm(p1, p2),
1778 eng1_mm_default_bt(p1, p2),
1779 eng1_mm_default_nbt_mim(p1, p2)
1780 {
1781 }
1782
~eng1_mm_default_mim(void)1783 eng1_mm_default_mim::~eng1_mm_default_mim(void)
1784 {
1785 }
1786
1787 /*################################################################################################*/
1788
1789 // eof
1790