1 // ENG1_SF.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_sf.h" // config.h is here -> we get ENABLE-macros here...
23
24 #include "v3d.h"
25 #include "seqbuild.h"
26
27 #include <fstream>
28 #include <sstream>
29 #include <algorithm>
30 using namespace std;
31
32 /*################################################################################################*/
33
34 // the surface area code apparently contains some bugs, since it sometimes
35 // crashes. another possibility is that the surface area math contains some
36 // bad cases (like arcs/segments with zero length/area ???) which should be
37 // avoided somehow. either way, the numerical and analytical gradients of
38 // surface area seem to match.
39
40 // LOWLIMIT is a cheat to prevent zero-divisions in surface-area calculations.
41 // if a zero-division seems to happen, the values are just changed to LOWLIMIT
42 // in as early stage as possible, thus making minimum effects on results...
43
44 #define LOWLIMIT 0.0000001 // 0.0000001 seems to work quite well...
45
46 /*################################################################################################*/
47
eng1_sf_param_SetDefaultValues(eng1_sf_param & prm)48 void eng1_sf_param_SetDefaultValues(eng1_sf_param & prm)
49 {
50 prm.vdwrad = 1.325;
51 prm.lenjon = 1.225;
52 prm.wescc = 1.000;
53 prm.wescd = 0.500;
54 prm.dipole1 = 2.0;
55 prm.dipole2 = 4.0;
56 prm.epsilon1 = 2.00;
57 prm.epsilon2 = 1.00;
58 prm.epsilon3 = 2.00;
59 prm.epsilon4 = 1.40;
60 prm.epsilon5 = 1.96;
61 prm.epsilon9 = 0.0015;
62 prm.exp_solv_1n = 0.15;
63 prm.exp_solv_1p = 1.00;
64 prm.exp_solv_2 = 1.00;
65 prm.imp_solv_1n = +4.00;
66 prm.imp_solv_1p = +2.00;
67 prm.imp_solv_2 = -12.00;
68 prm.solvrad = 0.15;
69 prm.wang = 0.85;
70 prm.wtor1 = 0.75;
71 prm.wtor2 = 0.50;
72 prm.wrep = 0.035;
73
74 prm.charge_wes[0] = +0.0744;
75 prm.charge_wes[1] = +0.0600;
76 prm.charge_wes[2] = -0.0488;
77 prm.charge_wes[3] = -0.0681;
78
79 prm.charge_sasa1[0] = -9.605;
80 prm.charge_sasa1[1] = -1.505;
81 prm.charge_sasa1[2] = +0.523;
82 prm.charge_sasa1[3] = +2.593;
83
84 prm.charge_sasa2[0] = +2.069;
85 prm.charge_sasa2[1] = -0.629;
86 prm.charge_sasa2[2] = -0.224;
87 prm.charge_sasa2[3] = -0.272;
88
89 prm.charge_pKa[0] = 7.830;
90 prm.charge_pKa[1] = 4.651;
91 prm.charge_pKa[2] = 11.5;
92 prm.charge_pKa[3] = 5.069;
93 prm.charge_pKa[4] = 9.5;
94 prm.charge_pKa[5] = 6.319;
95 prm.charge_pKa[6] = 6.088;
96 prm.charge_pKa[7] = 10.64;
97 prm.charge_pKa[8] = 10.81;
98
99 prm.charge_acid[0] = false;
100 prm.charge_acid[1] = true;
101 prm.charge_acid[2] = false;
102 prm.charge_acid[3] = true;
103 prm.charge_acid[4] = true;
104 prm.charge_acid[5] = true;
105 prm.charge_acid[6] = false;
106 prm.charge_acid[7] = false;
107 prm.charge_acid[8] = true;
108
109 // the pH setting is a default value that can be modified!!!
110 // the pH setting is a default value that can be modified!!!
111 // the pH setting is a default value that can be modified!!!
112
113 prm.pH = 7.0;
114 }
115
eng1_sf_param_MakeCopy(eng1_sf_param & sprm,eng1_sf_param & tprm)116 void eng1_sf_param_MakeCopy(eng1_sf_param & sprm, eng1_sf_param & tprm)
117 {
118 tprm.vdwrad = sprm.vdwrad;
119 tprm.lenjon = sprm.lenjon;
120 tprm.wescc = sprm.wescc;
121 tprm.wescd = sprm.wescd;
122 tprm.dipole1 = sprm.dipole1;
123 tprm.dipole2 = sprm.dipole2;
124 tprm.epsilon1 = sprm.epsilon1;
125 tprm.epsilon2 = sprm.epsilon2;
126 tprm.epsilon3 = sprm.epsilon3;
127 tprm.epsilon4 = sprm.epsilon4;
128 tprm.epsilon5 = sprm.epsilon5;
129 tprm.epsilon9 = sprm.epsilon9;
130 tprm.exp_solv_1n = sprm.exp_solv_1n;
131 tprm.exp_solv_1p = sprm.exp_solv_1p;
132 tprm.exp_solv_2 = sprm.exp_solv_2;
133 tprm.imp_solv_1n = sprm.imp_solv_1n;
134 tprm.imp_solv_1p = sprm.imp_solv_1p;
135 tprm.imp_solv_2 = sprm.imp_solv_2;
136 tprm.solvrad = sprm.solvrad;
137 tprm.wang = sprm.wang;
138 tprm.wtor1 = sprm.wtor1;
139 tprm.wtor2 = sprm.wtor2;
140 tprm.wrep = sprm.wrep;
141
142 for (int n1 = 0;n1 < 4;n1++)
143 {
144 tprm.charge_wes[n1] = sprm.charge_wes[n1];
145 tprm.charge_sasa1[n1] = sprm.charge_sasa1[n1];
146 tprm.charge_sasa2[n1] = sprm.charge_sasa2[n1];
147 }
148
149 for (int n1 = 0;n1 < 9;n1++)
150 {
151 tprm.charge_pKa[n1] = sprm.charge_pKa[n1];
152 tprm.charge_acid[n1] = sprm.charge_acid[n1];
153 }
154
155 tprm.pH = sprm.pH;
156 }
157
158 /*################################################################################################*/
159
sf_chn(void)160 sf_chn::sf_chn(void)
161 {
162 }
163
sf_chn(const sf_chn & p1)164 sf_chn::sf_chn(const sf_chn & p1)
165 {
166 res_vector = p1.res_vector;
167 }
168
~sf_chn(void)169 sf_chn::~sf_chn(void)
170 {
171 }
172
173 /*################################################################################################*/
174
sf_res(void)175 sf_res::sf_res(void)
176 {
177 symbol = '?';
178 for (i32s n1 = 0;n1 < 5;n1++)
179 {
180 peptide[n1] = NULL;
181 }
182
183 natm = 0;
184 atmr[0] = atmr[1] = atmr[2] = NULL;
185 loc_varind[0] = loc_varind[1] = loc_varind[2] = NOT_DEFINED;
186
187 state = SF_STATE_LOOP;
188 }
189
sf_res(char p1,atom * a1,atom * a2,atom * a3,atom * a4,atom * a5,i32s p2,atom * r1,atom * r2,atom * r3,i32s i1,i32s i2,i32s i3)190 sf_res::sf_res(char p1, atom * a1, atom * a2, atom * a3, atom * a4, atom * a5, i32s p2, atom * r1, atom * r2, atom * r3, i32s i1, i32s i2, i32s i3)
191 {
192 symbol = p1;
193 peptide[0] = a1;
194 peptide[1] = a2;
195 peptide[2] = a3;
196 peptide[3] = a4;
197 peptide[4] = a5;
198
199 natm = p2;
200 atmr[0] = r1;
201 atmr[1] = r2;
202 atmr[2] = r3;
203 loc_varind[0] = i1;
204 loc_varind[1] = i2;
205 loc_varind[2] = i3;
206
207 state = SF_STATE_LOOP;
208 }
209
sf_res(const sf_res & p1)210 sf_res::sf_res(const sf_res & p1)
211 {
212 symbol = p1.symbol;
213 for (i32s n1 = 0;n1 < 5;n1++)
214 {
215 peptide[n1] = p1.peptide[n1];
216 }
217
218 natm = p1.natm;
219 for (i32s n1 = 0;n1 < 3;n1++)
220 {
221 atmr[n1] = p1.atmr[n1];
222 loc_varind[n1] = p1.loc_varind[n1];
223 }
224
225 state = p1.state;
226 }
227
~sf_res(void)228 sf_res::~sf_res(void)
229 {
230 }
231
232 /*################################################################################################*/
233
setup1_sf(model * p1,SFmode SFmd,bool convert_crd)234 setup1_sf::setup1_sf(model * p1, SFmode SFmd, bool convert_crd) : setup(p1)
235 {
236 mode = SFmd;
237
238 eng1_sf_param_SetDefaultValues(prm);
239
240 // since the SF is the highest-level model, it must contain all parts of the system.
241 // take copies of all detected chains in the system, whether it was requested or not.
242
243 // initialize...
244
245 GetModel()->UpdateChains();
246 DefineSecondaryStructure(GetModel());
247 vector<chn_info> & ci_vector = (* GetModel()->ref_civ);
248
249 // the main loop; make contiguous chains like this:
250
251 // rA0 -- rB0 -- rC0 -- rD0 -- and_so_on
252 // | |
253 // rB1 rD1
254 // |
255 // rB2
256
257 fGL * crd_sf[3];
258 crd_sf[0] = new fGL[GetModel()->cs_vector.size() * 3];
259 crd_sf[1] = new fGL[GetModel()->cs_vector.size() * 3];
260 crd_sf[2] = new fGL[GetModel()->cs_vector.size() * 3];
261
262 i32s loc_varind_counter = 0;
263 for (i32u n1 = 0;n1 < ci_vector.size();n1++)
264 {
265 sf_chn newchn;
266 chn_vector.push_back(newchn);
267
268 if (ci_vector[n1].GetType() != chn_info::amino_acid) continue; // keep the same number of chains...
269
270 iter_al range1[2];
271 GetModel()->GetRange(1, n1, range1);
272
273 iter_al range2[2];
274 for (i32s n2 = 0;n2 < ci_vector[n1].length;n2++)
275 {
276 GetModel()->GetRange(2, range1, n2, range2);
277 if (range2[0] == range2[1]) { cout << "empty residue!!!" << endl; exit(EXIT_FAILURE); }
278
279 i32s tmp1[2]; tmp1[0] = (* range2[0]).builder_res_id >> 8;
280 for (tmp1[1] = 0;tmp1[1] < (i32s) model::amino_builder->resi_vector.size();tmp1[1]++)
281 {
282 if (model::amino_builder->resi_vector[tmp1[1]].id == tmp1[0]) break;
283 } if (tmp1[1] == (i32s) model::amino_builder->resi_vector.size()) continue;
284
285 char symbol_sf = model::amino_builder->resi_vector[tmp1[1]].symbol1;
286
287 i32s ind = 0; while (ind < 20 && model::sf_symbols[ind] != symbol_sf) ind++;
288 if (ind == 20) { cout << "BUG: unknown residue symbol." << endl; exit(EXIT_FAILURE); }
289
290 i32s natm_sf = NOT_DEFINED;
291 i32s atmi_sf[3] = { NOT_DEFINED, NOT_DEFINED, NOT_DEFINED };
292
293 switch (symbol_sf)
294 {
295 case 'A': natm_sf = 1; break;
296
297 case 'R': natm_sf = 3; break;
298
299 case 'N': natm_sf = 2; break;
300
301 case 'D': natm_sf = 2; break;
302
303 case 'C': natm_sf = 2; break;
304
305 case 'Q': natm_sf = 2; break;
306
307 case 'E': natm_sf = 2; break;
308
309 case 'G': natm_sf = 1; break;
310
311 case 'H': natm_sf = 2; break;
312
313 case 'I': natm_sf = 2; break;
314
315 case 'L': natm_sf = 2; break;
316
317 case 'K': natm_sf = 3; break;
318
319 case 'M': natm_sf = 2; break;
320
321 case 'F': natm_sf = 2; break;
322
323 case 'P': natm_sf = 1; break;
324
325 case 'S': natm_sf = 1; break;
326
327 case 'T': natm_sf = 1; break;
328
329 case 'W': natm_sf = 3; break;
330
331 case 'Y': natm_sf = 2; break;
332
333 case 'V': natm_sf = 1; break;
334
335 default:
336 cout << "BUG: unknown residue " << symbol_sf << " at aa2sf_Convert_sub1()." << endl;
337 exit(EXIT_FAILURE);
338 }
339
340 vector<i32s> idv;
341 idv.push_back(0x01);
342
343 atmi_sf[0] = idv.front();
344 for (i32u n3 = 0;n3 < GetModel()->cs_vector.size();n3++)
345 {
346 GetReducedCRD(range2, idv, & crd_sf[0][n3 * 3], n3);
347 }
348
349 if (natm_sf > 1)
350 {
351 idv.resize(0);
352 switch (symbol_sf)
353 {
354 case 'R':
355 idv.push_back(0x22); idv.push_back(0x21);
356 break;
357
358 case 'N':
359 idv.push_back(0x21);
360 break;
361
362 case 'D':
363 idv.push_back(0x21);
364 break;
365
366 case 'C':
367 idv.push_back(0x21);
368 break;
369
370 case 'Q':
371 idv.push_back(0x22); idv.push_back(0x21);
372 break;
373
374 case 'E':
375 idv.push_back(0x22); idv.push_back(0x21);
376 break;
377
378 case 'H':
379 idv.push_back(0x21); idv.push_back(0x22);
380 idv.push_back(0x23); idv.push_back(0x24); idv.push_back(0x25);
381 break;
382
383 case 'I':
384 idv.push_back(0x22);
385 break;
386
387 case 'L':
388 idv.push_back(0x21);
389 break;
390
391 case 'K':
392 idv.push_back(0x21);
393 break;
394
395 case 'M':
396 idv.push_back(0x22); idv.push_back(0x21); idv.push_back(0x23);
397 break;
398
399 case 'F':
400 idv.push_back(0x21); idv.push_back(0x24);
401 break;
402
403 case 'W':
404 idv.push_back(0x21); idv.push_back(0x23); idv.push_back(0x24);
405 break;
406
407 case 'Y':
408 idv.push_back(0x21); idv.push_back(0x24);
409 break;
410 }
411
412 atmi_sf[1] = idv.front();
413 for (i32u n3 = 0;n3 < GetModel()->cs_vector.size();n3++)
414 {
415 GetReducedCRD(range2, idv, & crd_sf[1][n3 * 3], n3);
416 }
417 }
418
419 if (natm_sf > 2)
420 {
421 idv.resize(0);
422 switch (symbol_sf)
423 {
424 case 'R':
425 idv.push_back(0x24);
426 break;
427
428 case 'K':
429 idv.push_back(0x23);
430 break;
431
432 case 'W':
433 idv.push_back(0x28); idv.push_back(0x25);
434 break;
435 }
436
437 atmi_sf[2] = idv.front();
438 for (i32u n3 = 0;n3 < GetModel()->cs_vector.size();n3++)
439 {
440 GetReducedCRD(range2, idv, & crd_sf[2][n3 * 3], n3);
441 }
442 }
443
444 // hide all atoms of the residue.
445
446 for (iter_al it1 = range2[0];it1 != range2[1];it1++)
447 {
448 (* it1).flags |= ATOMFLAG_IS_HIDDEN;
449 }
450
451 // make the residue and un-hide the simplified atoms.
452
453 atom * peptide_sf[5] = { NULL, NULL, NULL, NULL, NULL };
454
455 atom * atmr_sf[3] = { NULL, NULL, NULL };
456 i32s iloc_sf[3] = { NOT_DEFINED, NOT_DEFINED, NOT_DEFINED };
457
458 bool last_residue = (n2 == ci_vector[n1].length - 1);
459 for (iter_al it1 = range2[0];it1 != range2[1];it1++)
460 {
461 if (((* it1).builder_res_id & 0xFF) == 0x00) peptide_sf[0] = & (* it1); // N
462 if (((* it1).builder_res_id & 0xFF) == 0x01) peptide_sf[1] = & (* it1); // C
463 if (((* it1).builder_res_id & 0xFF) == 0x02) peptide_sf[2] = & (* it1); // C
464 if (((* it1).builder_res_id & 0xFF) == 0x10) peptide_sf[3] = & (* it1); // O
465 if (last_residue && ((* it1).builder_res_id & 0xFF) == 0x11) peptide_sf[4] = & (* it1); // O
466
467 if (atmi_sf[0] != NOT_DEFINED && ((* it1).builder_res_id & 0xFF) == atmi_sf[0]) // 1st
468 {
469 atmr_sf[0] = & (* it1);
470 iloc_sf[0] = loc_varind_counter++;
471 }
472
473 if (atmi_sf[1] != NOT_DEFINED && ((* it1).builder_res_id & 0xFF) == atmi_sf[1]) // 2nd
474 {
475 atmr_sf[1] = & (* it1);
476 iloc_sf[1] = loc_varind_counter++;
477 }
478
479 if (atmi_sf[2] != NOT_DEFINED && ((* it1).builder_res_id & 0xFF) == atmi_sf[2]) // 3rd
480 {
481 atmr_sf[2] = & (* it1);
482 iloc_sf[2] = loc_varind_counter++;
483 }
484 }
485
486 sf_res newres = sf_res(symbol_sf, peptide_sf[0], peptide_sf[1], peptide_sf[2], peptide_sf[3], peptide_sf[4], natm_sf, atmr_sf[0], atmr_sf[1], atmr_sf[2], iloc_sf[0], iloc_sf[1], iloc_sf[2]);
487 chn_vector.back().res_vector.push_back(newres);
488
489 fGL vdwr1 = -1.0; fGL vdwr2 = -1.0; fGL vdwr3 = -1.0;
490 fGL mass1 = -1.0; fGL mass2 = -1.0; fGL mass3 = -1.0;
491 switch (symbol_sf)
492 {
493 case 'A': // 71.1
494 vdwr1 = 0.193391;
495 mass1 = 71.1;
496 break;
497
498 case 'R': // 156.7
499 vdwr1 = 0.222079;
500 mass1 = 53.9; // 1.05
501
502 vdwr2 = 0.195669;
503 mass2 = 51.4; // 1.00
504
505 vdwr3 = 0.189176;
506 mass3 = 51.4; // 1.00
507 break;
508
509 case 'N': // 114.1
510 vdwr1 = 0.224929;
511 mass1 = 58.4; // 1.05
512
513 vdwr2 = 0.187591;
514 mass2 = 55.7; // 1.00
515 break;
516
517 case 'D': // 114.5
518 vdwr1 = 0.225038;
519 mass1 = 58.6; // 1.05
520
521 vdwr2 = 0.187136;
522 mass2 = 55.9; // 1.00
523 break;
524
525 case 'C': // 103.2
526 vdwr1 = 0.229061;
527 mass1 = 51.6; // 1.00
528
529 vdwr2 = 0.179623;
530 mass2 = 51.6; // 1.00
531 break;
532
533 case 'Q': // 128.2
534 vdwr1 = 0.226117;
535 mass1 = 64.1; // 1.00
536
537 vdwr2 = 0.204211;
538 mass2 = 64.1; // 1.00
539 break;
540
541 case 'E': // 128.7
542 vdwr1 = 0.222040;
543 mass1 = 64.3; // 1.00
544
545 vdwr2 = 0.191868;
546 mass2 = 64.3; // 1.00
547 break;
548
549 case 'G': // 57.1
550 vdwr1 = 0.173764;
551 mass1 = 57.1;
552 break;
553
554 case 'H': // 137.7
555 vdwr1 = 0.229728;
556 mass1 = 68.9; // 1.00
557
558 vdwr2 = 0.213733;
559 mass2 = 68.9; // 1.00
560 break;
561
562 case 'I': // 113.1
563 vdwr1 = 0.224550;
564 mass1 = 56.6; // 1.00
565
566 vdwr2 = 0.180242;
567 mass2 = 56.6; // 1.00
568 break;
569
570 case 'L': // 113.1
571 vdwr1 = 0.220761;
572 mass1 = 56.6; // 1.00
573
574 vdwr2 = 0.182951;
575 mass2 = 56.6; // 1.00
576 break;
577
578 case 'K': // 129.2
579 vdwr1 = 0.222174;
580 mass1 = 49.7; // 1.25
581
582 vdwr2 = 0.184924;
583 mass2 = 39.8; // 1.00
584
585 vdwr3 = 0.199069;
586 mass3 = 39.8; // 1.00
587 break;
588
589 case 'M': // 131.2
590 vdwr1 = 0.230562;
591 mass1 = 65.6; // 1.00
592
593 vdwr2 = 0.221798;
594 mass2 = 65.6; // 1.00
595 break;
596
597 case 'F': // 147.2
598 vdwr1 = 0.224414;
599 mass1 = 68.5; // 1.00
600
601 vdwr2 = 0.198277;
602 mass2 = 78.7; // 1.15
603 break;
604
605 case 'P': // 97.1
606 vdwr1 = 0.210783;
607 mass1 = 97.1;
608 break;
609
610 case 'S': // 87.1
611 vdwr1 = 0.204505;
612 mass1 = 87.1;
613 break;
614
615 case 'T': // 101.1
616 vdwr1 = 0.214132;
617 mass1 = 101.1;
618 break;
619
620 case 'W': // 188.2
621 vdwr1 = 0.232409;
622 mass1 = 62.7; // 1.00
623
624 vdwr2 = 0.233762;
625 mass2 = 62.7; // 1.00
626
627 vdwr3 = 0.225098;
628 mass3 = 62.7; // 1.00
629 break;
630
631 case 'Y': // 163.2
632 vdwr1 = 0.224633;
633 mass1 = 72.5; // 1.00
634
635 vdwr2 = 0.198270;
636 mass2 = 90.6; // 1.25
637 break;
638
639 case 'V': // 99.1
640 vdwr1 = 0.213718;
641 mass1 = 99.1;
642 break;
643
644 default:
645 cout << "BUG: unknown residue " << symbol_sf << " at setup1_sf ctor." << endl;
646 exit(EXIT_FAILURE);
647 }
648
649 if (atmr_sf[0] != NULL)
650 {
651 for (i32u n3 = 0;n3 < GetModel()->cs_vector.size();n3++)
652 {
653 fGL x = crd_sf[0][n3 * 3 + 0];
654 fGL y = crd_sf[0][n3 * 3 + 1];
655 fGL z = crd_sf[0][n3 * 3 + 2];
656
657 if (convert_crd) atmr_sf[0]->SetCRD(n3, x, y, z);
658 }
659
660 if (vdwr1 < 0.0 || mass1 < 0.0) { cout << "bad values for atom #1" << endl; exit(EXIT_FAILURE); }
661
662 atmr_sf[0]->vdwr = vdwr1;
663 atmr_sf[0]->mass = mass1;
664 }
665
666 if (atmr_sf[1] != NULL)
667 {
668 for (i32u n3 = 0;n3 < GetModel()->cs_vector.size();n3++)
669 {
670 fGL x = crd_sf[1][n3 * 3 + 0];
671 fGL y = crd_sf[1][n3 * 3 + 1];
672 fGL z = crd_sf[1][n3 * 3 + 2];
673
674 if (convert_crd) atmr_sf[1]->SetCRD(n3, x, y, z);
675 }
676
677 if (vdwr2 < 0.0 || mass2 < 0.0) { cout << "bad values for atom #2" << endl; exit(EXIT_FAILURE); }
678
679 atmr_sf[1]->vdwr = vdwr2;
680 atmr_sf[1]->mass = mass2;
681 }
682
683 if (atmr_sf[2] != NULL)
684 {
685 for (i32u n3 = 0;n3 < GetModel()->cs_vector.size();n3++)
686 {
687 fGL x = crd_sf[2][n3 * 3 + 0];
688 fGL y = crd_sf[2][n3 * 3 + 1];
689 fGL z = crd_sf[2][n3 * 3 + 2];
690
691 if (convert_crd) atmr_sf[2]->SetCRD(n3, x, y, z);
692 }
693
694 if (vdwr3 < 0.0 || mass3 < 0.0) { cout << "bad values for atom #3" << endl; exit(EXIT_FAILURE); }
695
696 atmr_sf[2]->vdwr = vdwr3;
697 atmr_sf[2]->mass = mass3;
698 }
699 }
700 }
701
702 // cout << "loc_varind_counter = " << loc_varind_counter << endl; // this MUST match...
703 // int stop; cin >> stop; // to GetSFAtomCount()!!!
704
705 delete[] crd_sf[0];
706 delete[] crd_sf[1];
707 delete[] crd_sf[2];
708
709 UpdateAtomFlags(); // THE NEW WAY!!!!!
710
711 // disulphide bridges...
712
713 for (i32u n1 = 0;n1 < ci_vector.size();n1++)
714 {
715 iter_al range1a[2];
716 GetModel()->GetRange(1, n1, range1a);
717
718 // intra-chain ones...
719
720 vector<i32s> cys_data1;
721 vector<atom *> cys_ref1;
722
723 for (i32s n2 = 0;n2 < ci_vector[n1].length;n2++)
724 {
725 if (ci_vector[n1].sequence1[n2] == 'C')
726 {
727 bool flag = true;
728
729 iter_al range1b[2];
730 GetModel()->GetRange(2, range1a, n2, range1b);
731 if (range1b[0] == range1b[1]) flag = false;
732
733 iter_al it1 = range1b[0];
734 while (it1 != range1b[1] && ((* it1).builder_res_id & 0xFF) != 0x21) it1++;
735 if (it1 == range1b[1]) flag = false;
736
737 if (flag)
738 {
739 cys_data1.push_back(n2);
740 cys_ref1.push_back(& (* it1));
741 }
742 }
743 }
744
745 for (i32s n2 = 0;n2 < ((i32s) cys_ref1.size()) - 1;n2++)
746 {
747 for (i32s n3 = n2 + 1;n3 < (i32s) cys_ref1.size();n3++)
748 {
749 bond tb = bond(cys_ref1[n2], cys_ref1[n3], bondtype('S'));
750 iter_bl it1 = find(GetModel()->bond_list.begin(), GetModel()->bond_list.end(), tb);
751 if (it1 != GetModel()->bond_list.end())
752 {
753 sf_dsb newdsb;
754 newdsb.chn[0] = n1; newdsb.res[0] = cys_data1[n2];
755 newdsb.chn[1] = n1; newdsb.res[1] = cys_data1[n3];
756 dsb_vector.push_back(newdsb);
757 }
758 }
759 }
760
761 // interchain ones...
762
763 for (i32u n2 = n1 + 1;n2 < ci_vector.size();n2++)
764 {
765 iter_al range2a[2];
766 GetModel()->GetRange(1, n2, range2a);
767
768 vector<i32s> cys_data2;
769 vector<atom *> cys_ref2;
770
771 for (i32s n3 = 0;n3 < ci_vector[n2].length;n3++)
772 {
773 if (ci_vector[n2].sequence1[n3] == 'C')
774 {
775 bool flag = true;
776
777 iter_al range2b[2];
778 GetModel()->GetRange(2, range2a, n3, range2b);
779 if (range2b[0] == range2b[1]) flag = false;
780
781 iter_al it1 = range2b[0];
782 while (it1 != range2b[1] && ((* it1).builder_res_id & 0xFF) != 0x21) it1++;
783 if (it1 == range2b[1]) flag = false;
784
785 if (flag)
786 {
787 cys_data2.push_back(n3);
788 cys_ref2.push_back(& (* it1));
789 }
790 }
791 }
792
793 for (i32u n3 = 0;n3 < cys_ref1.size();n3++)
794 {
795 for (i32u n4 = 0;n4 < cys_ref2.size();n4++)
796 {
797 bond tb = bond(cys_ref1[n3], cys_ref2[n4], bondtype('S'));
798 iter_bl it1 = find(GetModel()->bond_list.begin(), GetModel()->bond_list.end(), tb);
799 if (it1 != GetModel()->bond_list.end())
800 {
801 sf_dsb newdsb;
802 newdsb.chn[0] = n1; newdsb.res[0] = cys_data1[n3];
803 newdsb.chn[1] = n2; newdsb.res[1] = cys_data2[n4];
804 dsb_vector.push_back(newdsb);
805 }
806 }
807 }
808 }
809 }
810
811 // secondary structure...
812
813 // the default state is always loop. strands have the same places as in K&S.
814 // helices are shifted: the smallest helical peptide is "4...." -> "LHHHL"!!!!!
815
816 for (i32u n1 = 0;n1 < chn_vector.size();n1++)
817 {
818 if (ci_vector[n1].GetLength() != (i32s) chn_vector[n1].res_vector.size())
819 {
820 cout << "BUG: the sizes of chn_info and chn_vector differ!!!" << endl;
821 exit(EXIT_FAILURE);
822 }
823
824 const char * state = ci_vector[n1].GetSecStrStates();
825 for (i32s n2 = 0;n2 < (i32s) chn_vector[n1].res_vector.size();n2++)
826 {
827 chn_vector[n1].res_vector[n2].state = SF_STATE_LOOP; // set the default!!!
828
829 switch (state[n2])
830 {
831 case '4':
832 chn_vector[n1].res_vector[n2].state = SF_STATE_HELIX4;
833 break;
834
835 case 'S':
836 chn_vector[n1].res_vector[n2].state = SF_STATE_STRAND;
837 break;
838 }
839 }
840 }
841 }
842
~setup1_sf(void)843 setup1_sf::~setup1_sf(void)
844 {
845 // un-hide all atoms and reset the atomic radii and masses to their default values...
846
847 for (iter_al it1 = GetModel()->GetAtomsBegin();it1 != GetModel()->GetAtomsEnd();it1++)
848 {
849 (* it1).flags &= (~ATOMFLAG_IS_HIDDEN);
850
851 (* it1).vdwr = (* it1).el.GetVDWRadius();
852 (* it1).mass = (* it1).el.GetAtomicMass();
853 }
854 }
855
UpdateAtomFlags(void)856 void setup1_sf::UpdateAtomFlags(void)
857 {
858 // first hide all atoms, then later un-hide what is needed...
859 // first hide all atoms, then later un-hide what is needed...
860 // first hide all atoms, then later un-hide what is needed...
861
862 for (iter_al it1 = GetModel()->GetAtomsBegin();it1 != GetModel()->GetAtomsEnd();it1++)
863 {
864 (* it1).flags |= ATOMFLAG_IS_HIDDEN;
865 }
866
867 // handle chains...
868
869 for (i32u n1 = 0;n1 < chn_vector.size();n1++)
870 {
871 for (i32u n2 = 0;n2 < chn_vector[n1].res_vector.size();n2++)
872 {
873 for (i32s n3 = 0;n3 < chn_vector[n1].res_vector[n2].natm;n3++)
874 {
875 chn_vector[n1].res_vector[n2].atmr[n3]->flags |= ATOMFLAG_IS_SF_ATOM;
876 chn_vector[n1].res_vector[n2].atmr[n3]->flags &= (~ATOMFLAG_IS_HIDDEN);
877 }
878 }
879 }
880
881 // handle solvent molecules...
882
883 for (iter_al it1 = GetModel()->GetAtomsBegin();it1 != GetModel()->GetAtomsEnd();it1++)
884 {
885 if ((* it1).el.GetAtomicNumber() != 8) continue; // assume H2O!!!
886 if (!((* it1).flags & ATOMFLAG_IS_SOLVENT_ATOM)) continue; // assume H2O!!!
887
888 (* it1).flags |= ATOMFLAG_IS_SF_ATOM;
889 (* it1).flags &= (~ATOMFLAG_IS_HIDDEN);
890
891 (* it1).vdwr = 0.155;
892 (* it1).mass = 18.016;
893
894 // TODO : now assumes that oxygen coords = water coords.
895
896 }
897 }
898
GetReducedCRD(iter_al * iter,vector<i32s> & idv,fGL * crd,i32u cset)899 void setup1_sf::GetReducedCRD(iter_al * iter, vector<i32s> & idv, fGL * crd, i32u cset)
900 {
901 vector<atom *> refv;
902 for (i32u n1 = 0;n1 < idv.size();n1++)
903 {
904 iter_al it2 = iter[0];
905 while (it2 != iter[1] && ((* it2).builder_res_id & 0xFF) != idv[n1]) it2++;
906 if (it2 != iter[1]) refv.push_back(& (* it2));
907 }
908
909 if (!refv.size()) { cout << "BUG: no atoms found!" << endl; exit(EXIT_FAILURE); }
910
911 for (i32u n1 = 0;n1 < 3;n1++)
912 {
913 crd[n1] = 0.0;
914 for (i32u n2 = 0;n2 < refv.size();n2++)
915 {
916 const fGL * cdata = refv[n2]->GetCRD(cset);
917 crd[n1] += cdata[n1];
918 }
919
920 crd[n1] /= (f64) refv.size();
921 }
922 }
923
StorePStatesToModel(eng1_sf * eng)924 void setup1_sf::StorePStatesToModel(eng1_sf * eng)
925 {
926 if (!GetModel()->ref_civ) return;
927
928 vector<chn_info> & ci_vector = (* GetModel()->ref_civ);
929 if (chn_vector.size() != ci_vector.size())
930 {
931 cout << "ERROR : chain counts mismatch!" << endl;
932 exit(EXIT_FAILURE);
933 }
934
935 for (i32u n1 = 0;n1 < chn_vector.size();n1++)
936 {
937 if (!chn_vector[n1].res_vector.size()) continue; // empty chain -> nucleic acid.
938 if (chn_vector[n1].res_vector.size() != (i32u) ci_vector[n1].length)
939 {
940 cout << "ERROR : chain lengths mismatch!" << endl;
941 exit(EXIT_FAILURE);
942 }
943
944 // make the tables, if needed...
945 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
946
947 if (ci_vector[n1].p_state == NULL)
948 {
949 ci_vector[n1].p_state = new char[ci_vector[n1].length];
950 }
951
952 for (i32u n2 = 0;n2 < chn_vector[n1].res_vector.size();n2++)
953 {
954 bool flag_tc = false; // check for terminal charges...
955 if (n2 == 0 || ((i32s) n2) == ((i32s) chn_vector[n1].res_vector.size()) - 1)
956 {
957 fGL tmpc = chn_vector[n1].res_vector[n2].atmr[0]->charge;
958 if (tmpc < -0.5 || tmpc > +0.5) flag_tc = true;
959 }
960
961 int charge = 0;
962 for (i32s n3 = 0;n3 < chn_vector[n1].res_vector[n2].natm;n3++)
963 {
964 fGL tmpc = chn_vector[n1].res_vector[n2].atmr[n3]->charge;
965 if (tmpc < -0.5) charge -= 1; if (tmpc > +0.5) charge += 1;
966 }
967
968 // store the information...
969 // ^^^^^^^^^^^^^^^^^^^^^^^^
970
971 char flags = abs(charge);
972 if (charge < 0) flags |= PSTATE_SIGN_NEGATIVE;
973 else flags |= PSTATE_SIGN_POSITIVE;
974
975 if (flag_tc) flags |= PSTATE_CHARGED_TERMINAL;
976
977 ci_vector[n1].p_state[n2] = flags;
978 }
979 }
980 }
981
static_GetEngineCount(void)982 i32u setup1_sf::static_GetEngineCount(void)
983 {
984 return 1;
985 }
986
static_GetEngineIDNumber(i32u eng_index)987 i32u setup1_sf::static_GetEngineIDNumber(i32u eng_index)
988 {
989 if (eng_index < static_GetEngineCount()) return 0x01; // 0x01 is the default so far (we have only a single class).
990 else return (i32u) NOT_DEFINED; // error code...
991 }
992
static_GetEngineName(i32u eng_index)993 const char * setup1_sf::static_GetEngineName(i32u eng_index)
994 {
995 static char name[] = "simplified forcefield";
996
997 if (eng_index < static_GetEngineCount()) return name;
998 else return NULL; // error code...
999 }
1000
static_GetClassName(void)1001 const char * setup1_sf::static_GetClassName(void)
1002 {
1003 static char cn[] = "allsf";
1004 return cn;
1005 }
1006
GetEngineCount(void)1007 i32u setup1_sf::GetEngineCount(void)
1008 {
1009 return static_GetEngineCount();
1010 }
1011
GetEngineIDNumber(i32u eng_index)1012 i32u setup1_sf::GetEngineIDNumber(i32u eng_index)
1013 {
1014 return static_GetEngineIDNumber(eng_index);
1015 }
1016
GetEngineName(i32u eng_index)1017 const char * setup1_sf::GetEngineName(i32u eng_index)
1018 {
1019 return static_GetEngineName(eng_index);
1020 }
1021
GetClassName_lg(void)1022 const char * setup1_sf::GetClassName_lg(void)
1023 {
1024 return static_GetClassName();
1025 }
1026
CreateEngineByIndex(i32u eng_index)1027 engine * setup1_sf::CreateEngineByIndex(i32u eng_index)
1028 {
1029 if (eng_index >= GetEngineCount())
1030 {
1031 cout << "setup1_sf::CreateEngineByIndex() failed!" << endl;
1032 return NULL;
1033 }
1034
1035 GetModel()->use_periodic_boundary_conditions = false;
1036 if (GetModel()->use_boundary_potential) GetModel()->Message("use_boundary_potential = TRUE");
1037
1038 GetModel()->UpdateIndex();
1039 UpdateSetupInfo();
1040
1041 // return new eng1_sf(this, 1, true, false); // explicit...
1042 return new eng1_sf(this, 1, false, true); // implicit...
1043 }
1044
1045 /*################################################################################################*/
1046
1047 // here we are dependent on secondary structure constraints -> if secondary structure is modified,
1048 // the engine class must be discarded and a new one must be created to take effects into account...
1049
1050 // or take a closer look to CopyCRD() ?!?!?!? well this is not that urgent problem.
1051
1052 // IMPORTANT!!! constraints are no longer a part of energy calculations (since they practically affect them,
1053 // due to rounding errors)!!! you have to add constraints to energy to get the total energy!!!
1054
eng1_sf(setup * p1,i32u p2,bool esf,bool isf)1055 eng1_sf::eng1_sf(setup * p1, i32u p2, bool esf, bool isf) :
1056 engine(p1, p2),
1057 engine_bp(p1, p2)
1058 {
1059 use_explicit_solvent = esf;
1060 use_implicit_solvent = isf;
1061 if (use_explicit_solvent) GetSetup()->GetModel()->Message("using EXPLICIT solvent");
1062 if (use_implicit_solvent) GetSetup()->GetModel()->Message("using IMPLICIT solvent");
1063 if (use_explicit_solvent && use_implicit_solvent)
1064 {
1065 cout << "sorry, this is not yet supported!" << endl;
1066 exit(EXIT_FAILURE);
1067 }
1068
1069 if (GetSetup()->GetModel()->GetConstD_count() > 0)
1070 {
1071 GetSetup()->GetModel()->WarningMessage(CONSTRAINTS_NOT_SUPPORTED);
1072 }
1073
1074 // also check engine::bp_center!!!
1075 // also check engine::bp_center!!!
1076 // also check engine::bp_center!!!
1077 bp_fc_solute = 5000.0; // 50 kJ/(mol*�^2) = 5000 kJ/(mol*(nm)^2)
1078 bp_fc_solvent = 12500.0; // 125 kJ/(mol*�^2) = 12500 kJ/(mol*(nm)^2)
1079
1080 if (use_bp)
1081 {
1082 cout << "use_bp ; ";
1083 cout << bp_rad_solute << " " << bp_fc_solute << " ; ";
1084 cout << bp_rad_solvent << " " << bp_fc_solvent << endl;
1085 }
1086
1087 // create the local-to-global lookup table...
1088 // create the local-to-global lookup table...
1089 // create the local-to-global lookup table...
1090
1091 l2g_sf = new i32u[GetSetup()->GetSFAtomCount()];
1092
1093 atom ** atmtab = GetSetup()->GetSFAtoms();
1094 atom ** glob_atmtab = GetSetup()->GetAtoms();
1095 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount();n1++)
1096 {
1097 i32s index = 0;
1098 while (index < GetSetup()->GetAtomCount())
1099 {
1100 if (atmtab[n1] == glob_atmtab[index]) break;
1101 else index++;
1102 }
1103
1104 if (index >= GetSetup()->GetAtomCount())
1105 {
1106 cout << "BUG: eng1_sf ctor failed to create the l2g lookup table." << endl;
1107 exit(EXIT_FAILURE);
1108 }
1109
1110 l2g_sf[n1] = index;
1111 }
1112
1113 // various initialization tasks...
1114 // various initialization tasks...
1115 // various initialization tasks...
1116
1117 setup1_sf * mysu = dynamic_cast<setup1_sf *>(GetSetup());
1118 if (!mysu) { cout << "BUG: cast to setup1_sf failed." << endl; exit(EXIT_FAILURE); }
1119 myprm = & mysu->prm;
1120
1121 index_chn = new i32s[GetSetup()->GetSFAtomCount()];
1122 index_res = new i32s[GetSetup()->GetSFAtomCount()];
1123
1124 num_solvent = 0;
1125
1126 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount();n1++)
1127 {
1128 bool found = false;
1129
1130 for (i32u cc = 0;cc < mysu->chn_vector.size();cc++)
1131 {
1132 for (i32u rc = 0;rc < mysu->chn_vector[cc].res_vector.size();rc++)
1133 {
1134 for (i32s ac = 0;ac < mysu->chn_vector[cc].res_vector[rc].natm;ac++)
1135 {
1136 if (mysu->chn_vector[cc].res_vector[rc].atmr[ac] != atmtab[n1]) continue;
1137
1138 index_chn[n1] = cc;
1139 index_res[n1] = rc;
1140
1141 bool has_valid_varind = true;
1142 if (atmtab[n1]->varind < 0) has_valid_varind = false;
1143 if (atmtab[n1]->varind >= GetSetup()->GetAtomCount()) has_valid_varind = false;
1144 if (!has_valid_varind)
1145 {
1146 cout << "ERROR : invalid varind " << atmtab[n1]->varind << " was found!!!" << endl;
1147 exit(EXIT_FAILURE);
1148 }
1149
1150 if (n1 != mysu->chn_vector[cc].res_vector[rc].loc_varind[ac])
1151 {
1152 cout << "ERROR : an incorrect loc_varind was detected!!!" << endl;
1153 exit(EXIT_FAILURE);
1154 }
1155
1156 if (num_solvent != 0)
1157 {
1158 cout << "ERROR : ordering error : solvent atoms last!" << endl;
1159 exit(EXIT_FAILURE);
1160 }
1161
1162 found = true;
1163 break;
1164 }
1165
1166 if (found) break; // not necessary; an optimization for speed.
1167 }
1168
1169 if (found) break; // not necessary; an optimization for speed.
1170 }
1171
1172 if (!found) // if no match was found, this must be a solvent molecule...
1173 {
1174 if (!(atmtab[n1]->flags & ATOMFLAG_IS_SOLVENT_ATOM))
1175 {
1176 cout << "ERROR : an unknown (non-solvent) atom was found!!!" << endl;
1177 exit(EXIT_FAILURE);
1178 }
1179
1180 bool has_valid_varind = true;
1181 if (atmtab[n1]->varind < 0) has_valid_varind = false;
1182 if (atmtab[n1]->varind >= GetSetup()->GetAtomCount()) has_valid_varind = false;
1183 if (!has_valid_varind)
1184 {
1185 cout << "ERROR : invalid varind " << atmtab[n1]->varind << " was found!!!" << endl;
1186 exit(EXIT_FAILURE);
1187 }
1188
1189 index_chn[n1] = NOT_DEFINED;
1190 index_res[n1] = NOT_DEFINED;
1191 num_solvent++;
1192 }
1193 }
1194
1195 if (use_explicit_solvent && num_solvent == 0)
1196 {
1197 cout << "ERROR : explicit solvent requested but no solvent atoms detected!" << endl;
1198 exit(EXIT_FAILURE);
1199 }
1200
1201 constraints = 0.0;
1202 for (i32u n1 = 0;n1 < mysu->chn_vector.size();n1++)
1203 {
1204 for (i32u n2 = 0;n2 < mysu->chn_vector[n1].res_vector.size();n2++)
1205 {
1206 bool is_helix4 = false;
1207 if ((n2 - 1) >= 0 && mysu->chn_vector[n1].res_vector[n2 - 1].state == SF_STATE_HELIX4) is_helix4 = true;
1208 if ((n2 - 2) >= 0 && mysu->chn_vector[n1].res_vector[n2 - 2].state == SF_STATE_HELIX4) is_helix4 = true;
1209 if ((n2 - 3) >= 0 && mysu->chn_vector[n1].res_vector[n2 - 3].state == SF_STATE_HELIX4) is_helix4 = true;
1210
1211 if (is_helix4)
1212 {
1213 switch (mysu->chn_vector[n1].res_vector[n2].symbol)
1214 {
1215 case 'A': constraints += +1.446e+01; break;
1216 case 'R': constraints += -6.833e+00; break;
1217 case 'N': constraints += +9.452e+00; break;
1218 case 'D': constraints += +5.893e+00; break;
1219 case 'C': constraints += +1.779e+00; break;
1220 case 'Q': constraints += -1.667e+01; break;
1221 case 'E': constraints += +7.172e+00; break;
1222 case 'G': constraints += +1.646e+01; break;
1223 case 'H': constraints += +7.052e+00; break;
1224 case 'I': constraints += -1.236e+01; break;
1225 case 'L': constraints += +1.775e-01; break;
1226 case 'K': constraints += -8.890e-01; break;
1227 case 'M': constraints += +4.644e+00; break;
1228 case 'F': constraints += +5.412e+00; break;
1229 case 'P': constraints += -3.191e+00; break;
1230 case 'S': constraints += +1.971e+01; break;
1231 case 'T': constraints += +4.572e+00; break;
1232 case 'W': constraints += -6.210e+00; break;
1233 case 'Y': constraints += -1.238e+01; break;
1234 case 'V': constraints += -1.009e+01; break;
1235
1236 default:
1237 cout << "BUG: unknown residue (helix) : " << mysu->chn_vector[n1].res_vector[n2].symbol << endl;
1238 exit(EXIT_FAILURE);
1239 }
1240 }
1241
1242 if (mysu->chn_vector[n1].res_vector[n2].state == SF_STATE_STRAND)
1243 {
1244 switch (mysu->chn_vector[n1].res_vector[n2].symbol)
1245 {
1246 case 'A': constraints += -8.465e+00; break;
1247 case 'R': constraints += -7.428e+00; break;
1248 case 'N': constraints += -1.505e+01; break;
1249 case 'D': constraints += -9.706e+00; break;
1250 case 'C': constraints += -9.969e+00; break;
1251 case 'Q': constraints += -8.787e+00; break;
1252 case 'E': constraints += -9.782e+00; break;
1253 case 'G': constraints += -8.186e+00; break;
1254 case 'H': constraints += -1.019e+01; break;
1255 case 'I': constraints += -1.172e+01; break;
1256 case 'L': constraints += -1.168e+01; break;
1257 case 'K': constraints += -6.592e+00; break;
1258 case 'M': constraints += -1.007e+01; break;
1259 case 'F': constraints += -1.414e+01; break;
1260 case 'P': constraints += -9.055e+00; break;
1261 case 'S': constraints += -7.650e+00; break;
1262 case 'T': constraints += -1.405e+01; break;
1263 case 'W': constraints += -9.975e+00; break;
1264 case 'Y': constraints += -1.431e+01; break;
1265 case 'V': constraints += -1.285e+01; break;
1266
1267 default:
1268 cout << "BUG: unknown residue (strand) : " << mysu->chn_vector[n1].res_vector[n2].symbol << endl;
1269 exit(EXIT_FAILURE);
1270 }
1271 }
1272 }
1273 }
1274
1275 tmp_vartab = NULL;
1276 tmp_parames = NULL;
1277 tmp_paramsa1 = NULL;
1278 tmp_paramsa2 = NULL;
1279 tmp_newpKa = NULL;
1280
1281 mass = new f64[GetSetup()->GetSFAtomCount()];
1282 vdwr = new f64[GetSetup()->GetSFAtomCount()];
1283 charge1 = new f64[GetSetup()->GetSFAtomCount()];
1284 charge2 = new f64[GetSetup()->GetSFAtomCount()];
1285
1286 for (i32u n1 = 0;n1 < LAYERS;n1++)
1287 {
1288 vdwr1[n1] = new f64[GetSetup()->GetSFAtomCount() - num_solvent];
1289 vdwr2[n1] = new f64[GetSetup()->GetSFAtomCount() - num_solvent];
1290 sasaE[n1] = new f64[GetSetup()->GetSFAtomCount() - num_solvent];
1291
1292 solv_exp[n1] = new fGL[GetSetup()->GetSFAtomCount() - num_solvent];
1293 }
1294
1295 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount();n1++)
1296 {
1297 mass[n1] = atmtab[n1]->mass;
1298 vdwr[n1] = atmtab[n1]->vdwr * myprm->vdwrad;
1299
1300 charge1[n1] = 0.0; // SetupCharges() will handle this...
1301 charge2[n1] = 0.0; // SetupCharges() will handle this...
1302 }
1303
1304 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount() - num_solvent;n1++)
1305 {
1306 for (i32u n2 = 0;n2 < LAYERS;n2++)
1307 {
1308 vdwr1[n2][n1] = vdwr[n1] + ((f64) (n2 + 1)) * myprm->solvrad;
1309 vdwr2[n2][n1] = vdwr1[n2][n1] * vdwr1[n2][n1];
1310
1311 sasaE[n2][n1] = 0.0; // SetupCharges() will handle this...
1312 }
1313 }
1314 cout << "init ok; natm = " << GetSetup()->GetSFAtomCount() << endl;
1315
1316 /*##############################################*/
1317 /*##############################################*/
1318
1319 // main-chain t1-terms: direction is always from the N- to the C-terminal.
1320 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1321
1322 for (i32u n1 = 0;n1 < mysu->chn_vector.size();n1++)
1323 {
1324 for (i32s n2 = 0;n2 < ((i32s) mysu->chn_vector[n1].res_vector.size()) - 1;n2++)
1325 {
1326 sf_bt1 newbt1;
1327 newbt1.atmi[0] = mysu->chn_vector[n1].res_vector[n2 + 0].loc_varind[0];
1328 newbt1.atmi[1] = mysu->chn_vector[n1].res_vector[n2 + 1].loc_varind[0];
1329
1330 bool proline = (mysu->chn_vector[n1].res_vector[n2 + 1].symbol == 'P');
1331 newbt1.opt = (proline ? 0.352 : 0.38126); // assume 1/3 X-Pro cis...
1332
1333 newbt1.fc = 51.7e+03;
1334
1335 bt1_vector.push_back(newbt1);
1336 }
1337 }
1338 cout << "main bt1 ok; n = " << bt1_vector.size() << endl;
1339
1340 /*##############################################*/
1341 /*##############################################*/
1342
1343 // main-chain t2-terms: add like t1-terms -> directions for the t1-terms
1344 // will always be FALSE for the first one and TRUE for the second one.
1345 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1346
1347 for (i32s n1 = 0;n1 < ((i32s) bt1_vector.size()) - 1;n1++)
1348 {
1349 i32s test1 = bt1_vector[n1].GetIndex(0, false);
1350 i32s test2 = bt1_vector[n1 + 1].GetIndex(0, true);
1351
1352 if (test1 == test2)
1353 {
1354 sf_bt2 newbt2;
1355 newbt2.index1[0] = n1; newbt2.dir1[0] = false;
1356 newbt2.index1[1] = n1 + 1; newbt2.dir1[1] = true;
1357
1358 newbt2.atmi[0] = bt1_vector[newbt2.index1[0]].GetIndex(1, newbt2.dir1[0]);
1359 newbt2.atmi[1] = bt1_vector[newbt2.index1[0]].GetIndex(0, newbt2.dir1[0]);
1360 newbt2.atmi[2] = bt1_vector[newbt2.index1[1]].GetIndex(1, newbt2.dir1[1]);
1361
1362 i32s indc = index_chn[newbt2.atmi[1]];
1363 i32s indr = index_res[newbt2.atmi[1]];
1364
1365 newbt2.ttype = TTYPE_LOOP;
1366
1367 if ((indr - 1) >= 0 && mysu->chn_vector[indc].res_vector[indr - 1].state == SF_STATE_HELIX4) newbt2.ttype = TTYPE_HELIX;
1368 if ((indr - 2) >= 0 && mysu->chn_vector[indc].res_vector[indr - 2].state == SF_STATE_HELIX4) newbt2.ttype = TTYPE_HELIX;
1369 if ((indr - 3) >= 0 && mysu->chn_vector[indc].res_vector[indr - 3].state == SF_STATE_HELIX4) newbt2.ttype = TTYPE_HELIX;
1370
1371 if (mysu->chn_vector[indc].res_vector[indr].state == SF_STATE_STRAND) newbt2.ttype = TTYPE_STRAND;
1372
1373 switch (newbt2.ttype)
1374 {
1375 case TTYPE_HELIX:
1376 newbt2.opt = -0.0712161; // 94.1 deg
1377 newbt2.fc[0] = 181.7 * myprm->wang;
1378 newbt2.fc[1] = mysu->prm.wrep;
1379 break;
1380
1381 case TTYPE_STRAND:
1382 newbt2.opt = -0.543596; // 122.9 deg
1383 newbt2.fc[0] = 40.12 * myprm->wang;
1384 newbt2.fc[1] = mysu->prm.wrep;
1385 break;
1386
1387 default: // this is for TTYPE_LOOP...
1388 newbt2.opt = -0.319972; // 108.7 deg
1389 newbt2.fc[0] = 35.75 * myprm->wang;
1390 newbt2.fc[1] = mysu->prm.wrep;
1391 break;
1392 }
1393
1394 bt2_vector.push_back(newbt2);
1395 }
1396 }
1397 cout << "main bt2 ok; n = " << bt2_vector.size() << endl;
1398
1399 /*##############################################*/
1400 /*##############################################*/
1401
1402 // main-chain t3-terms: again add like t1-terms -> both t2-terms will
1403 // always have directions TRUE+TRUE -> implementation can take care of them.
1404 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1405
1406 for (i32s n1 = 0;n1 < ((i32s) bt2_vector.size()) - 1;n1++)
1407 {
1408 i32s test1i = bt2_vector[n1].index1[1];
1409 bool test1d = bt2_vector[n1].dir1[1];
1410
1411 i32s test2i = bt2_vector[n1 + 1].index1[0];
1412 bool test2d = bt2_vector[n1 + 1].dir1[0];
1413
1414 if (test1i == test2i && test1d != test2d)
1415 {
1416 sf_bt3 newbt3;
1417 newbt3.index2[0] = n1; newbt3.index2[1] = n1 + 1;
1418
1419 newbt3.index1[0] = bt2_vector[newbt3.index2[0]].index1[0];
1420 newbt3.dir1[0] = bt2_vector[newbt3.index2[0]].dir1[0];
1421
1422 newbt3.index1[1] = bt2_vector[newbt3.index2[0]].index1[1];
1423 newbt3.dir1[1] = bt2_vector[newbt3.index2[0]].dir1[1];
1424
1425 newbt3.index1[2] = bt2_vector[newbt3.index2[1]].index1[0];
1426 newbt3.dir1[2] = bt2_vector[newbt3.index2[1]].dir1[0];
1427
1428 newbt3.index1[3] = bt2_vector[newbt3.index2[1]].index1[1];
1429 newbt3.dir1[3] = bt2_vector[newbt3.index2[1]].dir1[1];
1430
1431 newbt3.atmi[0] = bt1_vector[newbt3.index1[0]].GetIndex(1, newbt3.dir1[0]);
1432 newbt3.atmi[1] = bt1_vector[newbt3.index1[0]].GetIndex(0, newbt3.dir1[0]);
1433 newbt3.atmi[2] = bt1_vector[newbt3.index1[3]].GetIndex(0, newbt3.dir1[3]);
1434 newbt3.atmi[3] = bt1_vector[newbt3.index1[3]].GetIndex(1, newbt3.dir1[3]);
1435
1436 ifstream file;
1437
1438 i32s indc1 = index_chn[newbt3.atmi[1]];
1439 i32s indr1 = index_res[newbt3.atmi[1]];
1440
1441 i32s indc2 = index_chn[newbt3.atmi[2]];
1442 i32s indr2 = index_res[newbt3.atmi[2]];
1443
1444 bool h4flag_0 = false; if ((indr1 - 0) >= 0 && mysu->chn_vector[indc1].res_vector[indr1 - 0].state == SF_STATE_HELIX4) h4flag_0 = true;
1445 bool h4flag_1 = false; if ((indr1 - 1) >= 0 && mysu->chn_vector[indc1].res_vector[indr1 - 1].state == SF_STATE_HELIX4) h4flag_1 = true;
1446 bool h4flag_2 = false; if ((indr1 - 2) >= 0 && mysu->chn_vector[indc1].res_vector[indr1 - 2].state == SF_STATE_HELIX4) h4flag_2 = true;
1447 bool h4flag_3 = false; if ((indr1 - 3) >= 0 && mysu->chn_vector[indc1].res_vector[indr1 - 3].state == SF_STATE_HELIX4) h4flag_3 = true;
1448
1449 bool ts1 = (mysu->chn_vector[indc1].res_vector[indr1].state == SF_STATE_STRAND);
1450 bool ts2 = (mysu->chn_vector[indc2].res_vector[indr2].state == SF_STATE_STRAND);
1451
1452 // tor_type... tor_type... tor_type... tor_type... tor_type...
1453 // tor_type... tor_type... tor_type... tor_type... tor_type...
1454 // tor_type... tor_type... tor_type... tor_type... tor_type...
1455
1456 newbt3.tor_ttype = TTYPE_LOOP;
1457 if (h4flag_1 || h4flag_2) newbt3.tor_ttype = TTYPE_HELIX;
1458 if (ts1 && ts2) newbt3.tor_ttype = TTYPE_STRAND;
1459
1460 switch (newbt3.tor_ttype)
1461 {
1462 case TTYPE_HELIX:
1463 newbt3.tors[0] = +0.807732;
1464 newbt3.tors[1] = 31.28; // OBSOLETE!!! * GetModel()->prm.wtor1;
1465 break;
1466
1467 case TTYPE_STRAND:
1468 newbt3.tors[0] = -2.80308;
1469 newbt3.tors[1] = 7.209; // OBSOLETE!!! * GetModel()->prm.wtor1;
1470 break;
1471
1472 default: // this is for TTYPE_LOOP...
1473 model::OpenLibDataFile(file, false, "param_sf/default/looptor.txt");
1474
1475 while (true)
1476 {
1477 if (file.peek() == 'e')
1478 {
1479 cout << "ERROR : end of file was reached at looptor.txt" << endl;
1480 exit(EXIT_FAILURE);
1481 }
1482
1483 char buffer[256]; char tp1; char tp2; file >> tp1 >> tp2;
1484 bool test1 = (tp1 != mysu->chn_vector[indc1].res_vector[indr1].symbol);
1485 bool test2 = (tp2 != mysu->chn_vector[indc2].res_vector[indr2].symbol);
1486 if (test1 || test2) file.getline(buffer, sizeof(buffer));
1487 else
1488 {
1489 f64 value;
1490 file >> value; newbt3.torc[0] = value * myprm->wtor1;
1491 file >> value; newbt3.torc[1] = value * myprm->wtor1;
1492 file >> value; newbt3.torc[2] = value * myprm->wtor1;
1493 file >> value; newbt3.tors[0] = value * myprm->wtor1;
1494 file >> value; newbt3.tors[1] = value * myprm->wtor1;
1495 file >> value; newbt3.tors[2] = value * myprm->wtor1;
1496 break; // exit the loop...
1497 }
1498 }
1499
1500 file.close(); // looptor.txt
1501 break;
1502 }
1503
1504 // dip_type... dip_type... dip_type... dip_type... dip_type...
1505 // dip_type... dip_type... dip_type... dip_type... dip_type...
1506 // dip_type... dip_type... dip_type... dip_type... dip_type...
1507
1508 newbt3.dip_ttype = TTYPE_LOOP;
1509 if (h4flag_1 || h4flag_2) newbt3.dip_ttype = TTYPE_HELIX;
1510 if (ts1 || ts2) newbt3.dip_ttype = TTYPE_STRAND;
1511
1512 switch (newbt3.dip_ttype)
1513 {
1514 case TTYPE_HELIX: break;
1515 case TTYPE_STRAND: break;
1516
1517 default: // this is for TTYPE_LOOP...
1518 model::OpenLibDataFile(file, false, "param_sf/default/loopdip.txt");
1519
1520 while (true)
1521 {
1522 if (file.peek() == 'e')
1523 {
1524 cout << "ERROR : end of file was reached at loopdip.txt" << endl;
1525 exit(EXIT_FAILURE);
1526 }
1527
1528 char buffer[256]; char tp1; char tp2; file >> tp1 >> tp2;
1529 bool test1 = (tp1 != mysu->chn_vector[indc1].res_vector[indr1].symbol);
1530 bool test2 = (tp2 != mysu->chn_vector[indc2].res_vector[indr2].symbol);
1531 if (test1 || test2) file.getline(buffer, sizeof(buffer));
1532 else
1533 {
1534 f64 value;
1535 file >> value; newbt3.dipc[0] = value;
1536 file >> value; newbt3.dipc[1] = value;
1537 file >> value; newbt3.dipc[2] = value;
1538 file >> value; newbt3.dips[0] = value;
1539 file >> value; newbt3.dips[1] = value;
1540 file >> value; newbt3.dips[2] = value;
1541 file >> value; newbt3.dipk[0] = value;
1542 file >> value; newbt3.dipk[1] = value;
1543 break; // exit the loop...
1544 }
1545 }
1546
1547 file.close(); // loopdip.txt
1548 break;
1549 }
1550
1551 // set the skip flag if this is a X-pro case !!!!!!!!!!!!!
1552 // set the skip flag if this is a X-pro case !!!!!!!!!!!!!
1553 // set the skip flag if this is a X-pro case !!!!!!!!!!!!!
1554
1555 newbt3.skip = (mysu->chn_vector[indc2].res_vector[indr2].symbol == 'P');
1556
1557 bt3_vector.push_back(newbt3);
1558 }
1559 }
1560 cout << "main bt3 ok; n = " << bt3_vector.size() << endl;
1561
1562 /*##############################################*/
1563 /*##############################################*/
1564
1565 // side-chain t1-terms: directions are always MAIN -> SIDE1 -> SIDE2.
1566 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1567
1568 i32u bt1_mark_side = bt1_vector.size();
1569 for (i32u n1 = 0;n1 < mysu->chn_vector.size();n1++)
1570 {
1571 for (i32u n2 = 0;n2 < mysu->chn_vector[n1].res_vector.size();n2++)
1572 {
1573 for (i32s n3 = 0;n3 < mysu->chn_vector[n1].res_vector[n2].natm - 1;n3++)
1574 {
1575 sf_bt1 newbt1;
1576 newbt1.atmi[0] = mysu->chn_vector[n1].res_vector[n2].loc_varind[n3];
1577 newbt1.atmi[1] = mysu->chn_vector[n1].res_vector[n2].loc_varind[n3 + 1];
1578
1579 if (!n3)
1580 {
1581 switch (mysu->chn_vector[n1].res_vector[n2].symbol)
1582 {
1583 case 'R':
1584 newbt1.opt = 0.29852; newbt1.fc = 6.24e+03;
1585 break;
1586
1587 case 'N':
1588 newbt1.opt = 0.25417; newbt1.fc = 38.5e+03;
1589 break;
1590
1591 case 'D':
1592 newbt1.opt = 0.25474; newbt1.fc = 38.3e+03;
1593 break;
1594
1595 case 'C':
1596 newbt1.opt = 0.27946; newbt1.fc = 32.6e+03;
1597 break;
1598
1599 case 'Q':
1600 newbt1.opt = 0.29621; newbt1.fc = 5.28e+03;
1601 break;
1602
1603 case 'E':
1604 newbt1.opt = 0.29818; newbt1.fc = 5.48e+03;
1605 break;
1606
1607 case 'H':
1608 newbt1.opt = 0.35533; newbt1.fc = 21.4e+03;
1609 break;
1610
1611 case 'I':
1612 newbt1.opt = 0.25265; newbt1.fc = 49.4e+03;
1613 break;
1614
1615 case 'L':
1616 newbt1.opt = 0.26045; newbt1.fc = 32.9e+03;
1617 break;
1618
1619 case 'K':
1620 newbt1.opt = 0.25549; newbt1.fc = 48.1e+03;
1621 break;
1622
1623 case 'M':
1624 newbt1.opt = 0.35232; newbt1.fc = 1.24e+03;
1625 break;
1626
1627 case 'F':
1628 newbt1.opt = 0.37875; newbt1.fc = 27.7e+03;
1629 break;
1630
1631 case 'W':
1632 newbt1.opt = 0.34120; newbt1.fc = 19.9e+03;
1633 break;
1634
1635 case 'Y':
1636 newbt1.opt = 0.37809; newbt1.fc = 19.6e+03;
1637 break;
1638
1639 default:
1640 cout << "problems: unknown residue (S-T1-A) : " << mysu->chn_vector[n1].res_vector[n2].symbol << endl;
1641 exit(EXIT_FAILURE);
1642 }
1643 }
1644 else
1645 {
1646 switch (mysu->chn_vector[n1].res_vector[n2].symbol)
1647 {
1648 case 'R':
1649 newbt1.opt = 0.28736; newbt1.fc = 9.51e+03;
1650 break;
1651
1652 case 'K':
1653 newbt1.opt = 0.25333; newbt1.fc = 35.4e+03;
1654 break;
1655
1656 case 'W':
1657 newbt1.opt = 0.21051; newbt1.fc = 60.0e+03; // modified fc!!!
1658 break;
1659
1660 default:
1661 cout << "problems: unknown residue (S-T1-B) : " << mysu->chn_vector[n1].res_vector[n2].symbol << endl;
1662 exit(EXIT_FAILURE);
1663 }
1664 }
1665
1666 bt1_vector.push_back(newbt1);
1667 }
1668 }
1669 }
1670 cout << "side bt1 ok; n = " << bt1_vector.size() - bt1_mark_side << endl;
1671
1672 /*##############################################*/
1673 /*##############################################*/
1674
1675 // side-chain t2-terms: directions are always MAIN -> SIDE1 -> SIDE2.
1676 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1677
1678 i32u bt2_mark_side = bt2_vector.size();
1679 for (i32s n1 = (i32s) bt1_mark_side;n1 < ((i32s) bt1_vector.size()) - 1;n1++)
1680 {
1681 i32s test1 = bt1_vector[n1].GetIndex(0, false);
1682 i32s test2 = bt1_vector[n1 + 1].GetIndex(0, true);
1683
1684 i32s indc = index_chn[test1];
1685 i32s indr = index_res[test1];
1686
1687 bool sidechain = (mysu->chn_vector[indc].res_vector[indr].atmr[0] != atmtab[test1]);
1688 if (test1 == test2 && sidechain)
1689 {
1690 sf_bt2 newbt2;
1691 newbt2.index1[0] = n1; newbt2.dir1[0] = false;
1692 newbt2.index1[1] = n1 + 1; newbt2.dir1[1] = true;
1693
1694 newbt2.atmi[0] = bt1_vector[newbt2.index1[0]].GetIndex(1, newbt2.dir1[0]);
1695 newbt2.atmi[1] = bt1_vector[newbt2.index1[0]].GetIndex(0, newbt2.dir1[0]);
1696 newbt2.atmi[2] = bt1_vector[newbt2.index1[1]].GetIndex(1, newbt2.dir1[1]);
1697
1698 newbt2.ttype = TTYPE_SIDE;
1699
1700 switch (mysu->chn_vector[indc].res_vector[indr].symbol)
1701 {
1702 case 'R':
1703 newbt2.opt = -0.62392;
1704 newbt2.fc[0] = 14.7;
1705 newbt2.fc[1] = 0.0;
1706 break;
1707
1708 case 'K':
1709 newbt2.opt = -0.85737;
1710 newbt2.fc[0] = 9.65;
1711 newbt2.fc[1] = 0.0;
1712 break;
1713
1714 case 'W':
1715 newbt2.opt = -0.48517;
1716 newbt2.fc[0] = 68.5;
1717 newbt2.fc[1] = 0.0;
1718 break;
1719
1720 default:
1721 cout << "problems: unknown residue (S-T2) : " << mysu->chn_vector[indc].res_vector[indr].symbol << endl;
1722 exit(EXIT_FAILURE);
1723 }
1724
1725 bt2_vector.push_back(newbt2);
1726 }
1727 }
1728 cout << "side bt2 ok; n = " << bt2_vector.size() - bt2_mark_side << endl;
1729
1730 /*##############################################*/
1731 /*##############################################*/
1732
1733 // t4-terms for all non-terminal residues with side-chains: just add
1734 // t1- and t2-subterms like defined before; all directions are TRUE.
1735 // ^^^^^^^^^^^^^^^^^^^^^^^^
1736
1737 for (i32u n1 = 0;n1 < mysu->chn_vector.size();n1++)
1738 {
1739 for (i32s n2 = 1;n2 < ((i32s) mysu->chn_vector[n1].res_vector.size()) - 1;n2++)
1740 {
1741 if (mysu->chn_vector[n1].res_vector[n2].natm > 1)
1742 {
1743 i32s index[2];
1744 index[0] = mysu->chn_vector[n1].res_vector[n2].loc_varind[0];
1745 index[1] = mysu->chn_vector[n1].res_vector[n2].loc_varind[1];
1746
1747 sf_bt4 newbt4;
1748 newbt4.index1 = 0; // find the first side-chain t1-term:
1749 while (newbt4.index1 < (i32s) bt1_vector.size())
1750 {
1751 bool test1 = (bt1_vector[newbt4.index1].atmi[0] == index[0]);
1752 bool test2 = (bt1_vector[newbt4.index1].atmi[1] == index[1]);
1753 if (test1 && test2) break; else newbt4.index1++;
1754 }
1755
1756 if (newbt4.index1 >= (i32s) bt1_vector.size())
1757 {
1758 cout << "ERROR : bt4/index1 search failed!!!" << endl;
1759 exit(EXIT_FAILURE);
1760 }
1761
1762 newbt4.index2 = 0; // find the first t2-term (main-chain):
1763 while (newbt4.index2 < (i32s) bt2_vector.size())
1764 {
1765 bool test1 = (bt2_vector[newbt4.index2].atmi[1] == index[0]);
1766 if (test1) break; else newbt4.index2++;
1767 }
1768
1769 if (newbt4.index2 >= (i32s) bt2_vector.size())
1770 {
1771 cout << "ERROR : bt4/index2 search failed!!!" << endl;
1772 exit(EXIT_FAILURE);
1773 }
1774
1775 char symbol = mysu->chn_vector[n1].res_vector[n2].symbol;
1776 i32s state = mysu->chn_vector[n1].res_vector[n2].state;
1777
1778 // changed quickly 20050608 ; checkme!!!
1779 if ((n2 - 1) >= 0 && mysu->chn_vector[n1].res_vector[n2 - 1].state == SF_STATE_HELIX4) state = SF_STATE_HELIX4;
1780 if ((n2 - 2) >= 0 && mysu->chn_vector[n1].res_vector[n2 - 2].state == SF_STATE_HELIX4) state = SF_STATE_HELIX4;
1781 if ((n2 - 3) >= 0 && mysu->chn_vector[n1].res_vector[n2 - 3].state == SF_STATE_HELIX4) state = SF_STATE_HELIX4;
1782
1783 switch (symbol)
1784 {
1785 case 'R':
1786 newbt4.opt = 0.82733; newbt4.fc = 78.42;
1787
1788 switch (state)
1789 {
1790 case SF_STATE_HELIX4:
1791 newbt4.fscos[0] = +1.66715; newbt4.fscos[1] = +2.52406; newbt4.fscos[2] = -0.81090;
1792 newbt4.fssin[0] = -0.44044; newbt4.fssin[1] = -0.56621; newbt4.fssin[2] = +0.16812;
1793 break;
1794
1795 case SF_STATE_STRAND:
1796 newbt4.fscos[0] = +2.81092; newbt4.fscos[1] = -2.18907; newbt4.fscos[2] = -1.31969;
1797 newbt4.fssin[0] = +0.80986; newbt4.fssin[1] = +0.21052; newbt4.fssin[2] = -0.07480;
1798 break;
1799
1800 default:
1801 newbt4.fscos[0] = +1.84105; newbt4.fscos[1] = -0.34928; newbt4.fscos[2] = -1.37635;
1802 newbt4.fssin[0] = -1.07776; newbt4.fssin[1] = +0.50658; newbt4.fssin[2] = +0.43117;
1803 }
1804 break;
1805
1806 case 'N':
1807 newbt4.opt = 0.79233; newbt4.fc = 135.9;
1808
1809 switch (state)
1810 {
1811 case SF_STATE_HELIX4:
1812 newbt4.fscos[0] = +2.06084; newbt4.fscos[1] = +2.21810; newbt4.fscos[2] = -2.68616;
1813 newbt4.fssin[0] = -1.72997; newbt4.fssin[1] = +0.13154; newbt4.fssin[2] = +1.59189;
1814 break;
1815
1816 case SF_STATE_STRAND:
1817 newbt4.fscos[0] = +3.12145; newbt4.fscos[1] = -3.27636; newbt4.fscos[2] = -3.00767;
1818 newbt4.fssin[0] = +0.96246; newbt4.fssin[1] = +0.62799; newbt4.fssin[2] = +0.05736;
1819 break;
1820
1821 default:
1822 newbt4.fscos[0] = +1.74298; newbt4.fscos[1] = -0.79241; newbt4.fscos[2] = -2.61712;
1823 newbt4.fssin[0] = -1.03069; newbt4.fssin[1] = +0.19134; newbt4.fssin[2] = +1.25520;
1824 }
1825 break;
1826
1827 case 'D':
1828 newbt4.opt = 0.79856; newbt4.fc = 136.2;
1829
1830 switch (state)
1831 {
1832 case SF_STATE_HELIX4:
1833 newbt4.fscos[0] = +1.85269; newbt4.fscos[1] = +2.31000; newbt4.fscos[2] = -2.26981;
1834 newbt4.fssin[0] = -2.09388; newbt4.fssin[1] = -0.13391; newbt4.fssin[2] = +1.83416;
1835 break;
1836
1837 case SF_STATE_STRAND:
1838 newbt4.fscos[0] = +3.39290; newbt4.fscos[1] = -3.94822; newbt4.fscos[2] = -2.65580;
1839 newbt4.fssin[0] = +1.70965; newbt4.fssin[1] = -0.36173; newbt4.fssin[2] = -0.43294;
1840 break;
1841
1842 default:
1843 newbt4.fscos[0] = +1.59474; newbt4.fscos[1] = -0.74304; newbt4.fscos[2] = -2.43083;
1844 newbt4.fssin[0] = -0.82258; newbt4.fssin[1] = -0.02264; newbt4.fssin[2] = +1.50055;
1845 }
1846 break;
1847
1848 case 'C':
1849 newbt4.opt = 0.77013; newbt4.fc = 134.1;
1850
1851 switch (state)
1852 {
1853 case SF_STATE_HELIX4:
1854 newbt4.fscos[0] = +1.29742; newbt4.fscos[1] = +2.82739; newbt4.fscos[2] = -2.06328;
1855 newbt4.fssin[0] = -1.07293; newbt4.fssin[1] = +0.12694; newbt4.fssin[2] = +1.99976;
1856 break;
1857
1858 case SF_STATE_STRAND:
1859 newbt4.fscos[0] = +2.17623; newbt4.fscos[1] = -2.82029; newbt4.fscos[2] = -3.48759;
1860 newbt4.fssin[0] = +0.45237; newbt4.fssin[1] = +0.60356; newbt4.fssin[2] = -0.52859;
1861 break;
1862
1863 default:
1864 newbt4.fscos[0] = +1.39788; newbt4.fscos[1] = -0.52709; newbt4.fscos[2] = -2.47451;
1865 newbt4.fssin[0] = -0.95942; newbt4.fssin[1] = +0.51342; newbt4.fssin[2] = +0.98760;
1866 }
1867 break;
1868
1869 case 'Q':
1870 newbt4.opt = 0.81477; newbt4.fc = 80.30;
1871
1872 switch (state)
1873 {
1874 case SF_STATE_HELIX4:
1875 newbt4.fscos[0] = +1.71509; newbt4.fscos[1] = +1.93731; newbt4.fscos[2] = -0.76456;
1876 newbt4.fssin[0] = -1.01086; newbt4.fssin[1] = -1.05898; newbt4.fssin[2] = +0.48023;
1877 break;
1878
1879 case SF_STATE_STRAND:
1880 newbt4.fscos[0] = +2.80955; newbt4.fscos[1] = -2.04764; newbt4.fscos[2] = -1.83362;
1881 newbt4.fssin[0] = +0.63332; newbt4.fssin[1] = +0.24795; newbt4.fssin[2] = -0.05034;
1882 break;
1883
1884 default:
1885 newbt4.fscos[0] = +1.85194; newbt4.fscos[1] = -0.29906; newbt4.fscos[2] = -1.32592;
1886 newbt4.fssin[0] = -1.41994; newbt4.fssin[1] = +0.29379; newbt4.fssin[2] = +0.35692;
1887 }
1888 break;
1889
1890 case 'E':
1891 newbt4.opt = 0.81811; newbt4.fc = 97.26;
1892
1893 switch (state)
1894 {
1895 case SF_STATE_HELIX4:
1896 newbt4.fscos[0] = +1.44402; newbt4.fscos[1] = +1.94219; newbt4.fscos[2] = -0.80589;
1897 newbt4.fssin[0] = -1.01587; newbt4.fssin[1] = -1.09656; newbt4.fssin[2] = -0.07502;
1898 break;
1899
1900 case SF_STATE_STRAND:
1901 newbt4.fscos[0] = +2.83636; newbt4.fscos[1] = -2.33477; newbt4.fscos[2] = -1.91387;
1902 newbt4.fssin[0] = +0.41781; newbt4.fssin[1] = -0.03389; newbt4.fssin[2] = +0.16658;
1903 break;
1904
1905 default:
1906 newbt4.fscos[0] = +1.61321; newbt4.fscos[1] = -0.33262; newbt4.fscos[2] = -1.16010;
1907 newbt4.fssin[0] = -1.13465; newbt4.fssin[1] = -0.04564; newbt4.fssin[2] = +0.28664;
1908 }
1909 break;
1910
1911 case 'H':
1912 newbt4.opt = 0.68790; newbt4.fc = 104.7;
1913
1914 switch (state)
1915 {
1916 case SF_STATE_HELIX4:
1917 newbt4.fscos[0] = +2.53650; newbt4.fscos[1] = +2.73410; newbt4.fscos[2] = -3.45313;
1918 newbt4.fssin[0] = +0.20418; newbt4.fssin[1] = +0.76434; newbt4.fssin[2] = +2.08214;
1919 break;
1920
1921 case SF_STATE_STRAND:
1922 newbt4.fscos[0] = +2.26449; newbt4.fscos[1] = -2.53277; newbt4.fscos[2] = -4.92778;
1923 newbt4.fssin[0] = +0.01335; newbt4.fssin[1] = +0.70585; newbt4.fssin[2] = -0.44344;
1924 break;
1925
1926 default:
1927 newbt4.fscos[0] = +1.83546; newbt4.fscos[1] = -0.26818; newbt4.fscos[2] = -3.69531;
1928 newbt4.fssin[0] = -1.45726; newbt4.fssin[1] = +0.46062; newbt4.fssin[2] = +0.96875;
1929 }
1930 break;
1931
1932 case 'I':
1933 newbt4.opt = 0.84384; newbt4.fc = 108.5;
1934
1935 switch (state)
1936 {
1937 case SF_STATE_HELIX4:
1938 newbt4.fscos[0] = +0.36523; newbt4.fscos[1] = +2.54766; newbt4.fscos[2] = -1.29027;
1939 newbt4.fssin[0] = -2.97542; newbt4.fssin[1] = -0.24663; newbt4.fssin[2] = +2.89687;
1940 break;
1941
1942 case SF_STATE_STRAND:
1943 newbt4.fscos[0] = +2.00042; newbt4.fscos[1] = -2.35794; newbt4.fscos[2] = -3.63292;
1944 newbt4.fssin[0] = -0.46656; newbt4.fssin[1] = +2.35817; newbt4.fssin[2] = -0.45473;
1945 break;
1946
1947 default:
1948 newbt4.fscos[0] = +0.67483; newbt4.fscos[1] = -0.72266; newbt4.fscos[2] = -2.26144;
1949 newbt4.fssin[0] = -1.90307; newbt4.fssin[1] = +0.81662; newbt4.fssin[2] = +0.64670;
1950 }
1951 break;
1952
1953 case 'L':
1954 newbt4.opt = 0.82350; newbt4.fc = 175.9;
1955
1956 switch (state)
1957 {
1958 case SF_STATE_HELIX4:
1959 newbt4.fscos[0] = +3.84136; newbt4.fscos[1] = +3.32244; newbt4.fscos[2] = -1.22385;
1960 newbt4.fssin[0] = -1.23669; newbt4.fssin[1] = -0.52039; newbt4.fssin[2] = +1.13885;
1961 break;
1962
1963 case SF_STATE_STRAND:
1964 newbt4.fscos[0] = +4.76904; newbt4.fscos[1] = -2.43156; newbt4.fscos[2] = -3.14965;
1965 newbt4.fssin[0] = +0.65700; newbt4.fssin[1] = +0.38307; newbt4.fssin[2] = -0.30443;
1966 break;
1967
1968 default:
1969 newbt4.fscos[0] = +3.52823; newbt4.fscos[1] = +0.38402; newbt4.fscos[2] = -1.72298;
1970 newbt4.fssin[0] = -2.37699; newbt4.fssin[1] = -0.03600; newbt4.fssin[2] = +0.40297;
1971 }
1972 break;
1973
1974 case 'K':
1975 newbt4.opt = 0.79880; newbt4.fc = 144.6;
1976
1977 switch (state)
1978 {
1979 case SF_STATE_HELIX4:
1980 newbt4.fscos[0] = +1.91050; newbt4.fscos[1] = +2.79925; newbt4.fscos[2] = -1.45299;
1981 newbt4.fssin[0] = -0.47100; newbt4.fssin[1] = -0.31985; newbt4.fssin[2] = +1.06281;
1982 break;
1983
1984 case SF_STATE_STRAND:
1985 newbt4.fscos[0] = +2.98681; newbt4.fscos[1] = -1.90303; newbt4.fscos[2] = -2.31907;
1986 newbt4.fssin[0] = +0.31173; newbt4.fssin[1] = +0.01677; newbt4.fssin[2] = -0.22673;
1987 break;
1988
1989 default:
1990 newbt4.fscos[0] = +1.96145; newbt4.fscos[1] = +0.12819; newbt4.fscos[2] = -1.91454;
1991 newbt4.fssin[0] = -1.33542; newbt4.fssin[1] = +0.29320; newbt4.fssin[2] = +0.55337;
1992 }
1993 break;
1994
1995 case 'M':
1996 newbt4.opt = 0.83147; newbt4.fc = 37.65;
1997
1998 switch (state)
1999 {
2000 case SF_STATE_HELIX4:
2001 newbt4.fscos[0] = +1.88761; newbt4.fscos[1] = +0.43592; newbt4.fscos[2] = +0.43654;
2002 newbt4.fssin[0] = -0.59080; newbt4.fssin[1] = -0.24729; newbt4.fssin[2] = -0.97472;
2003 break;
2004
2005 case SF_STATE_STRAND:
2006 newbt4.fscos[0] = +2.95738; newbt4.fscos[1] = -1.50891; newbt4.fscos[2] = -0.86063;
2007 newbt4.fssin[0] = +0.77162; newbt4.fssin[1] = +0.23999; newbt4.fssin[2] = +0.00057;
2008 break;
2009
2010 default:
2011 newbt4.fscos[0] = +2.17883; newbt4.fscos[1] = -0.87283; newbt4.fscos[2] = -0.23810;
2012 newbt4.fssin[0] = -0.78868; newbt4.fssin[1] = +0.60167; newbt4.fssin[2] = +0.11778;
2013 }
2014 break;
2015
2016 case 'F':
2017 newbt4.opt = 0.67581; newbt4.fc = 114.0;
2018
2019 switch (state)
2020 {
2021 case SF_STATE_HELIX4:
2022 newbt4.fscos[0] = +3.73550; newbt4.fscos[1] = +3.12113; newbt4.fscos[2] = -3.65848;
2023 newbt4.fssin[0] = +0.46830; newbt4.fssin[1] = +1.08106; newbt4.fssin[2] = +1.64272;
2024 break;
2025
2026 case SF_STATE_STRAND:
2027 newbt4.fscos[0] = +2.62637; newbt4.fscos[1] = -3.22908; newbt4.fscos[2] = -5.53782;
2028 newbt4.fssin[0] = -0.34849; newbt4.fssin[1] = +0.99799; newbt4.fssin[2] = -0.06089;
2029 break;
2030
2031 default:
2032 newbt4.fscos[0] = +2.14212; newbt4.fscos[1] = -0.32047; newbt4.fscos[2] = -3.57612;
2033 newbt4.fssin[0] = -1.37620; newbt4.fssin[1] = +0.72569; newbt4.fssin[2] = +0.51568;
2034 }
2035 break;
2036
2037 case 'W':
2038 newbt4.opt = 0.70518; newbt4.fc = 111.5;
2039
2040 switch (state)
2041 {
2042 case SF_STATE_HELIX4:
2043 newbt4.fscos[0] = +2.41235; newbt4.fscos[1] = +2.56594; newbt4.fscos[2] = -3.14962;
2044 newbt4.fssin[0] = +0.41312; newbt4.fssin[1] = +0.98169; newbt4.fssin[2] = +1.69113;
2045 break;
2046
2047 case SF_STATE_STRAND:
2048 newbt4.fscos[0] = +2.23143; newbt4.fscos[1] = -3.10454; newbt4.fscos[2] = -4.07617;
2049 newbt4.fssin[0] = +0.32407; newbt4.fssin[1] = +0.78458; newbt4.fssin[2] = -0.90925;
2050 break;
2051
2052 default:
2053 newbt4.fscos[0] = +1.89729; newbt4.fscos[1] = -0.94277; newbt4.fscos[2] = -3.46953;
2054 newbt4.fssin[0] = -1.08009; newbt4.fssin[1] = +1.02104; newbt4.fssin[2] = +1.06076;
2055 }
2056 break;
2057
2058 case 'Y':
2059 newbt4.opt = 0.67080; newbt4.fc = 106.7;
2060
2061 switch (state)
2062 {
2063 case SF_STATE_HELIX4:
2064 newbt4.fscos[0] = +3.25507; newbt4.fscos[1] = +3.38981; newbt4.fscos[2] = -3.76346;
2065 newbt4.fssin[0] = +0.41418; newbt4.fssin[1] = +0.76818; newbt4.fssin[2] = +1.97564;
2066 break;
2067
2068 case SF_STATE_STRAND:
2069 newbt4.fscos[0] = +2.36446; newbt4.fscos[1] = -2.66255; newbt4.fscos[2] = -5.36649;
2070 newbt4.fssin[0] = -0.19421; newbt4.fssin[1] = +0.99025; newbt4.fssin[2] = -0.25182;
2071 break;
2072
2073 default:
2074 newbt4.fscos[0] = +2.07289; newbt4.fscos[1] = -0.22634; newbt4.fscos[2] = -3.81184;
2075 newbt4.fssin[0] = -1.35353; newbt4.fssin[1] = +0.74664; newbt4.fssin[2] = +0.72466;
2076 }
2077 break;
2078
2079 default:
2080 cout << "problems: unknown residue (S-T4) : " << symbol << endl;
2081 exit(EXIT_FAILURE);
2082 }
2083
2084 newbt4.fc *= myprm->wang;
2085
2086 newbt4.fscos[0] *= myprm->wtor2;
2087 newbt4.fscos[1] *= myprm->wtor2;
2088 newbt4.fscos[2] *= myprm->wtor2;
2089
2090 newbt4.fssin[0] *= myprm->wtor2;
2091 newbt4.fssin[1] *= myprm->wtor2;
2092 newbt4.fssin[2] *= myprm->wtor2;
2093
2094 bt4_vector.push_back(newbt4);
2095 }
2096 }
2097 }
2098 cout << "bt4 ok; n = " << bt4_vector.size() << endl;
2099
2100 /*##############################################*/
2101 /*##############################################*/
2102
2103 // dsb-terms, both bt1 and bt2...
2104
2105 for (i32u n1 = 0;n1 < mysu->dsb_vector.size();n1++)
2106 {
2107 i32s index1[2];
2108 index1[0] = mysu->chn_vector[mysu->dsb_vector[n1].chn[0]].res_vector[mysu->dsb_vector[n1].res[0]].loc_varind[0];
2109 index1[1] = mysu->chn_vector[mysu->dsb_vector[n1].chn[0]].res_vector[mysu->dsb_vector[n1].res[0]].loc_varind[1];
2110
2111 i32s ind1 = 0;
2112 while (ind1 < (i32s) bt1_vector.size())
2113 {
2114 bool test1 = (bt1_vector[ind1].atmi[0] == index1[0]);
2115 bool test2 = (bt1_vector[ind1].atmi[1] == index1[1]);
2116 if (test1 && test2) break;
2117
2118 ind1++;
2119 }
2120
2121 if (ind1 >= (i32s) bt1_vector.size())
2122 {
2123 cout << "DSB error #1" << endl;
2124 exit(EXIT_FAILURE);
2125 }
2126
2127 i32s index2[2];
2128 index2[0] = mysu->chn_vector[mysu->dsb_vector[n1].chn[1]].res_vector[mysu->dsb_vector[n1].res[1]].loc_varind[0];
2129 index2[1] = mysu->chn_vector[mysu->dsb_vector[n1].chn[1]].res_vector[mysu->dsb_vector[n1].res[1]].loc_varind[1];
2130
2131 i32s ind2 = 0;
2132 while (ind2 < (i32s) bt1_vector.size())
2133 {
2134 bool test1 = (bt1_vector[ind2].atmi[0] == index2[0]);
2135 bool test2 = (bt1_vector[ind2].atmi[1] == index2[1]);
2136 if (test1 && test2) break;
2137
2138 ind2++;
2139 }
2140
2141 if (ind2 >= (i32s) bt1_vector.size())
2142 {
2143 cout << "DSB error #2" << endl;
2144 exit(EXIT_FAILURE);
2145 }
2146
2147 sf_bt1 newbt1;
2148 newbt1.atmi[0] = index1[1];
2149 newbt1.atmi[1] = index2[1];
2150
2151 newbt1.opt = 0.203251; newbt1.fc = 60.0e+03; // modified fc!!!
2152
2153 i32s ind3 = bt1_vector.size();
2154 bt1_vector.push_back(newbt1);
2155
2156 sf_bt2 newbt2;
2157 newbt2.ttype = TTYPE_BRIDGE;
2158 newbt2.opt = -0.236514;
2159 newbt2.fc[0] = 26.0;
2160 newbt2.fc[1] = 0.0;
2161
2162 // first angle...
2163
2164 newbt2.index1[0] = ind1; newbt2.dir1[0] = false;
2165 newbt2.index1[1] = ind3; newbt2.dir1[1] = true;
2166
2167 newbt2.atmi[0] = bt1_vector[newbt2.index1[0]].GetIndex(1, newbt2.dir1[0]);
2168 newbt2.atmi[1] = bt1_vector[newbt2.index1[0]].GetIndex(0, newbt2.dir1[0]);
2169 newbt2.atmi[2] = bt1_vector[newbt2.index1[1]].GetIndex(1, newbt2.dir1[1]);
2170
2171 bt2_vector.push_back(newbt2);
2172
2173 // second angle...
2174
2175 newbt2.index1[0] = ind2; newbt2.dir1[0] = false;
2176 newbt2.index1[1] = ind3; newbt2.dir1[1] = false;
2177
2178 newbt2.atmi[0] = bt1_vector[newbt2.index1[0]].GetIndex(1, newbt2.dir1[0]);
2179 newbt2.atmi[1] = bt1_vector[newbt2.index1[0]].GetIndex(0, newbt2.dir1[0]);
2180 newbt2.atmi[2] = bt1_vector[newbt2.index1[1]].GetIndex(1, newbt2.dir1[1]);
2181
2182 bt2_vector.push_back(newbt2);
2183 }
2184 cout << "dsb ok." << endl;
2185
2186 /*##############################################*/
2187 /*##############################################*/
2188
2189 // terminal angles...
2190
2191 for (i32u n1 = 0;n1 < mysu->chn_vector.size();n1++)
2192 {
2193 if (mysu->chn_vector[n1].res_vector[0].natm > 1)
2194 {
2195 i32s index[3];
2196 index[0] = mysu->chn_vector[n1].res_vector[0].loc_varind[1];
2197 index[1] = mysu->chn_vector[n1].res_vector[0].loc_varind[0];
2198 index[2] = mysu->chn_vector[n1].res_vector[1].loc_varind[0];
2199
2200 i32s ind1 = 0;
2201 while (ind1 < (i32s) bt1_vector.size())
2202 {
2203 bool test1 = (bt1_vector[ind1].atmi[0] == index[1]);
2204 bool test2 = (bt1_vector[ind1].atmi[1] == index[0]);
2205 if (test1 && test2) break;
2206
2207 ind1++;
2208 }
2209
2210 if (ind1 >= (i32s) bt1_vector.size())
2211 {
2212 cout << "NT error #1" << endl;
2213 exit(EXIT_FAILURE);
2214 }
2215
2216 i32s ind2 = 0;
2217 while (ind2 < (i32s) bt1_vector.size())
2218 {
2219 bool test1 = (bt1_vector[ind2].atmi[0] == index[1]);
2220 bool test2 = (bt1_vector[ind2].atmi[1] == index[2]);
2221 if (test1 && test2) break;
2222
2223 ind2++;
2224 }
2225
2226 if (ind2 >= (i32s) bt1_vector.size())
2227 {
2228 cout << "NT error #2" << endl;
2229 exit(EXIT_FAILURE);
2230 }
2231
2232 sf_bt2 newbt2;
2233 newbt2.index1[0] = ind1; newbt2.dir1[0] = true;
2234 newbt2.index1[1] = ind2; newbt2.dir1[1] = true;
2235
2236 newbt2.atmi[0] = bt1_vector[newbt2.index1[0]].GetIndex(1, newbt2.dir1[0]);
2237 newbt2.atmi[1] = bt1_vector[newbt2.index1[0]].GetIndex(0, newbt2.dir1[0]);
2238 newbt2.atmi[2] = bt1_vector[newbt2.index1[1]].GetIndex(1, newbt2.dir1[1]);
2239
2240 newbt2.ttype = TTYPE_TERM;
2241 newbt2.opt = -0.328360;
2242 newbt2.fc[0] = 12.0;
2243 newbt2.fc[1] = myprm->wrep;
2244
2245 bt2_vector.push_back(newbt2);
2246 }
2247
2248 i32s sz = mysu->chn_vector[n1].res_vector.size();
2249 if (mysu->chn_vector[n1].res_vector[sz - 1].natm > 1)
2250 {
2251 i32s index[3];
2252 index[0] = mysu->chn_vector[n1].res_vector[sz - 1].loc_varind[1];
2253 index[1] = mysu->chn_vector[n1].res_vector[sz - 1].loc_varind[0];
2254 index[2] = mysu->chn_vector[n1].res_vector[sz - 2].loc_varind[0];
2255
2256 i32s ind1 = 0;
2257 while (ind1 < (i32s) bt1_vector.size())
2258 {
2259 bool test1 = (bt1_vector[ind1].atmi[0] == index[1]);
2260 bool test2 = (bt1_vector[ind1].atmi[1] == index[0]);
2261 if (test1 && test2) break;
2262
2263 ind1++;
2264 }
2265
2266 if (ind1 >= (i32s) bt1_vector.size())
2267 {
2268 cout << "CT error #1" << endl;
2269 exit(EXIT_FAILURE);
2270 }
2271
2272 i32s ind2 = 0;
2273 while (ind2 < (i32s) bt1_vector.size())
2274 {
2275 bool test1 = (bt1_vector[ind2].atmi[0] == index[2]);
2276 bool test2 = (bt1_vector[ind2].atmi[1] == index[1]);
2277 if (test1 && test2) break;
2278
2279 ind2++;
2280 }
2281
2282 if (ind2 >= (i32s) bt1_vector.size())
2283 {
2284 cout << "CT error #2" << endl;
2285 exit(EXIT_FAILURE);
2286 }
2287
2288 sf_bt2 newbt2;
2289 newbt2.index1[0] = ind1; newbt2.dir1[0] = true;
2290 newbt2.index1[1] = ind2; newbt2.dir1[1] = false;
2291
2292 newbt2.atmi[0] = bt1_vector[newbt2.index1[0]].GetIndex(1, newbt2.dir1[0]);
2293 newbt2.atmi[1] = bt1_vector[newbt2.index1[0]].GetIndex(0, newbt2.dir1[0]);
2294 newbt2.atmi[2] = bt1_vector[newbt2.index1[1]].GetIndex(1, newbt2.dir1[1]);
2295
2296 newbt2.ttype = TTYPE_TERM;
2297 newbt2.opt = -0.339739;
2298 newbt2.fc[0] = 20.5;
2299 newbt2.fc[1] = myprm->wrep;
2300
2301 bt2_vector.push_back(newbt2);
2302 }
2303 }
2304 cout << "term ok." << endl;
2305
2306 /*##############################################*/
2307 /*##############################################*/
2308
2309 // nb1-terms...
2310
2311 vector <sf_nonbonded_lookup> lookup;
2312
2313 ifstream file;
2314 model::OpenLibDataFile(file, false, "param_sf/default/nonbonded.txt");
2315
2316 while (file.peek() != 'e')
2317 {
2318 char s1; i32u a1;
2319 char s2; i32u a2;
2320 f64 opte;
2321
2322 file >> s1 >> a1 >> s2 >> a2 >> opte;
2323
2324 char buffer[256];
2325 file.getline(buffer, sizeof(buffer));
2326
2327 sf_nonbonded_lookup newlu;
2328 newlu.s1 = s1; newlu.a1 = a1;
2329 newlu.s2 = s2; newlu.a2 = a2;
2330 newlu.opte = opte;
2331
2332 // cout << "adding a new lookup record : " << s1 << " " << a1 << " " << s2 << " " << a2 << " " << opte << endl;
2333
2334 lookup.push_back(newlu);
2335 }
2336
2337 file.close(); // nonbonded.txt
2338
2339 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount() - 1;n1++)
2340 {
2341 for (i32s n2 = n1 + 1;n2 < GetSetup()->GetSFAtomCount();n2++)
2342 {
2343 sf_nbt1 newnbt1;
2344 newnbt1.atmi[0] = n1; // this is a local index...
2345 newnbt1.atmi[1] = n2; // this is a local index...
2346
2347 if (!InitNBT1(& newnbt1, lookup)) continue;
2348 else nbt1_vector.push_back(newnbt1);
2349 }
2350 }
2351 cout << "nbt1 ok; n = " << nbt1_vector.size() << endl;
2352
2353 /*##############################################*/
2354 /*##############################################*/
2355
2356 // the solvent term... nbt1 are needed for all VA's separated by more than 1 bonds !!!
2357 // the solvent term... nbt1 are needed for all VA's separated by more than 1 bonds !!!
2358 // the solvent term... nbt1 are needed for all VA's separated by more than 1 bonds !!!
2359
2360 for (i32u n1 = 0;n1 < LAYERS;n1++)
2361 {
2362 nbt3_nl[n1] = new sf_nbt3_nl[GetSetup()->GetSFAtomCount() - num_solvent];
2363
2364 for (i32s n2 = 0;n2 < GetSetup()->GetSFAtomCount() - num_solvent;n2++)
2365 {
2366 bool skip = false;
2367
2368 // do not skip any of the layers since solv_exp[] is used and should be up-to-date. see also SetupCharges()
2369 // do not skip any of the layers since solv_exp[] is used and should be up-to-date. see also SetupCharges()
2370 // do not skip any of the layers since solv_exp[] is used and should be up-to-date. see also SetupCharges()
2371
2372 if (skip) /*nbt3_nl[n1][n2].index = NULL*/ exit(EXIT_FAILURE);
2373 else nbt3_nl[n1][n2].index = new i32s[size_nl[n1]];
2374 }
2375 }
2376
2377 dist1 = new i32s[GetSetup()->GetSFAtomCount() - num_solvent];
2378 dist2 = new f64[(GetSetup()->GetSFAtomCount() - num_solvent) * ((GetSetup()->GetSFAtomCount() - num_solvent) - 1) / 2];
2379
2380 i32s n1 = 0; i32s n2 = 0;
2381 while (n2 < GetSetup()->GetSFAtomCount() - num_solvent)
2382 {
2383 dist1[n2++] = n1;
2384 n1 += (GetSetup()->GetSFAtomCount() - num_solvent) - n2;
2385 }
2386
2387 /*##############################################*/
2388 /*##############################################*/
2389
2390 // allocate the arrays for intermediate results...
2391 // allocate the arrays for intermediate results...
2392 // allocate the arrays for intermediate results...
2393
2394 bt1data = new sf_bt1_data[bt1_vector.size()];
2395 bt2data = new sf_bt2_data[bt2_vector.size()];
2396
2397 // finally update the charges...
2398
2399 CopyCRD(GetSetup()->GetModel(), this, 0); // WHICH COORDINATES?!?!? define an "active" crd-set in the model class?!?!?!
2400 SetupCharges();
2401 }
2402
~eng1_sf(void)2403 eng1_sf::~eng1_sf(void)
2404 {
2405 delete[] l2g_sf;
2406 delete[] index_chn;
2407 delete[] index_res;
2408
2409 delete[] mass;
2410 delete[] vdwr;
2411 delete[] charge1;
2412 delete[] charge2;
2413
2414 for (i32u n1 = 0;n1 < LAYERS;n1++)
2415 {
2416 delete[] vdwr1[n1];
2417 delete[] vdwr2[n1];
2418 delete[] sasaE[n1];
2419
2420 delete[] solv_exp[n1];
2421
2422 for (i32s n2 = 0;n2 < GetSetup()->GetSFAtomCount() - num_solvent;n2++)
2423 {
2424 if (!nbt3_nl[n1][n2].index) continue;
2425 else delete[] nbt3_nl[n1][n2].index;
2426 }
2427
2428 delete[] nbt3_nl[n1];
2429 }
2430
2431 delete[] dist1;
2432 delete[] dist2;
2433
2434 delete[] bt1data;
2435 delete[] bt2data;
2436
2437 if (tmp_vartab != NULL)
2438 {
2439 delete[] tmp_vartab; tmp_vartab = NULL;
2440 delete[] tmp_parames; tmp_parames = NULL;
2441 delete[] tmp_paramsa1; tmp_paramsa1 = NULL;
2442 delete[] tmp_paramsa2; tmp_paramsa2 = NULL;
2443 delete[] tmp_newpKa; tmp_newpKa = NULL;
2444 }
2445 }
2446
SetTorsionConstraint(atom *,atom *,atom *,atom *,f64,f64,bool)2447 bool eng1_sf::SetTorsionConstraint(atom *, atom *, atom *, atom *, f64, f64, bool)
2448 {
2449 return false;
2450 }
2451
RemoveTorsionConstraint(atom *,atom *,atom *,atom *)2452 bool eng1_sf::RemoveTorsionConstraint(atom *, atom *, atom *, atom *)
2453 {
2454 return false;
2455 }
2456
2457 // HOW TO PROPERLY DETERMINE THE STATES???? DEPENDENT ON pH AND/OR GEOMETRY?!?!?!?
2458 // HOW TO PROPERLY DETERMINE THE STATES???? DEPENDENT ON pH AND/OR GEOMETRY?!?!?!?
2459 // HOW TO PROPERLY DETERMINE THE STATES???? DEPENDENT ON pH AND/OR GEOMETRY?!?!?!?
2460
SetupCharges(void)2461 void eng1_sf::SetupCharges(void)
2462 {
2463 // first, we calculate energy just to get the sasa values...
2464
2465 Compute(0);
2466
2467 // initialize the charges...
2468
2469 vector<i32s> variables;
2470 atom ** atmtab = GetSetup()->GetSFAtoms();
2471
2472 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount();n1++)
2473 {
2474 charge1[n1] = 0.0;
2475
2476 i32s cgi; i32s cvi;
2477 GetChgGrpVar(n1, cgi, cvi);
2478 if (cvi < 0) continue;
2479
2480 variables.push_back(n1);
2481 }
2482
2483 const i32s smoothing = 10;
2484
2485 if (tmp_vartab != NULL)
2486 {
2487 delete[] tmp_vartab;
2488 delete[] tmp_parames;
2489 delete[] tmp_paramsa1;
2490 delete[] tmp_paramsa2;
2491 delete[] tmp_newpKa;
2492 }
2493
2494 tmp_vartab = new i32s[variables.size() + 1];
2495 tmp_parames = new f64[variables.size()];
2496 tmp_paramsa1 = new f64[variables.size()];
2497 tmp_paramsa2 = new f64[variables.size()];
2498 tmp_newpKa = new f64[variables.size()];
2499
2500 for (i32u n1 = 0;n1 < variables.size();n1++) // copy the contents of variables...
2501 {
2502 tmp_vartab[n1] = variables[n1];
2503 } tmp_vartab[variables.size()] = -1; // terminate by -1!!!
2504
2505 f64 * tmp_charge = new f64[variables.size()];
2506 f64 * tmp_smooth = new f64[variables.size() * smoothing];
2507
2508 // for acid-type groups: HA <--> H+ + A- // product neg. charged; products more charged than reactants.
2509 // for base-type groups: BH+ <--> H+ + B // reactant pos. charged; equal charges on reactants and products.
2510
2511 for (i32u n1 = 0;n1 < variables.size();n1++)
2512 {
2513 i32s cgi; i32s cvi;
2514 GetChgGrpVar(variables[n1], cgi, cvi);
2515 bool is_acid = myprm->charge_acid[cvi];
2516 bool is_base = !myprm->charge_acid[cvi];
2517 f64 pKa = myprm->charge_pKa[cvi];
2518
2519 // pKa = pH + log([HX]/[X-])
2520 // log([X-]/[HX]) = pH - pKa = log(x/(1-x))
2521 // x = 10^(pH - pKa) / (1 + 10^(pH - pKa))
2522
2523 f64 tmp1 = pow(10.0, myprm->pH - pKa);
2524 f64 tmp2 = tmp1 / (1.0 + tmp1);
2525
2526 if (is_acid) charge1[variables[n1]] = -tmp2;
2527 if (is_base) charge1[variables[n1]] = 1.0 - tmp2;
2528
2529 for (i32s n2 = 0;n2 < smoothing;n2++)
2530 {
2531 tmp_smooth[n1 * smoothing + n2] = charge1[variables[n1]];
2532 }
2533 }
2534
2535 // the loop...
2536
2537 //ofstream logf("/tmp/pKa.log", ios::out);
2538 for (i32s loop = 0;loop < 250;loop++)
2539 {
2540 for (i32u n1 = 0;n1 < variables.size();n1++)
2541 {
2542 i32s cgi; i32s cvi;
2543 GetChgGrpVar(variables[n1], cgi, cvi);
2544 bool is_acid = myprm->charge_acid[cvi];
2545 bool is_base = !myprm->charge_acid[cvi];
2546 f64 pKa = myprm->charge_pKa[cvi];
2547
2548 f64 wes = myprm->charge_wes[cgi];
2549 f64 sasa1 = myprm->charge_sasa1[cgi];
2550 f64 sasa2 = myprm->charge_sasa2[cgi];
2551
2552 f64 k1 = pow(10.0, -pKa);
2553 f64 dg = -8.31451 * 300.0 * log(k1);
2554
2555 f64 chargeV = -1.0;
2556 if (is_base) chargeV = +1.0;
2557
2558 f64 parames = 0.0;
2559 f64 * crdV = & crd[l2g_sf[variables[n1]] * 3];
2560 for (i32s n2 = 0;n2 < GetSetup()->GetSFAtomCount();n2++)
2561 {
2562 if (variables[n1] == n2) continue;
2563 if (charge1[n2] == 0.0) continue;
2564
2565 f64 * crdA = & crd[l2g_sf[n2] * 3];
2566
2567 f64 r2 = 0.0;
2568 for (i32s n9 = 0;n9 < 3;n9++)
2569 {
2570 f64 tmp1 = crdV[n9] - crdA[n9];
2571 r2 += tmp1 * tmp1;
2572 }
2573
2574 if (r2 < 0.01) continue; // always skip if r < 1 � !!!
2575 f64 r1 = sqrt(r2);
2576
2577 // if (r1 > 1.5) continue; // also skip large distances; charge neutralization???
2578
2579 // e = 2 + 76 * (( r / A ) ^ n) / (1 + ( r / A ) ^ n) = 2 + 76 * f / g
2580
2581 f64 e1 = solv_exp[0][variables[n1]];
2582 f64 e2 = solv_exp[0][n2];
2583
2584 f64 eps_n = myprm->epsilon1 + r2 * myprm->epsilon2; // assume constant!!!
2585 f64 eps_A = myprm->epsilon3 - log(1.0 + (e1 + e2) * myprm->epsilon4 + e1 * e2 * myprm->epsilon5); // assume constant!!!
2586
2587 f64 t7a = r1 / eps_A;
2588 f64 t7b = pow(t7a, eps_n);
2589 f64 t7c = 1.0 + t7b;
2590 f64 t7d = 2.0 + 76.0 * (t7b / t7c);
2591
2592 // f = Q/(e*r)
2593 // df/dr = -Q * ((1/e*r^2) + (de/dr)/(e^2*r))
2594
2595 f64 qq = 138.9354518 * chargeV * charge1[n2];
2596 parames += qq / (t7d * r1);
2597 }
2598
2599 tmp_parames[n1] = parames;
2600
2601 f64 paramsa1 = solv_exp[0][variables[n1]] + solv_exp[1][variables[n1]] + solv_exp[2][variables[n1]]; // sum : 3*0.5 = 1.5
2602 f64 paramsa2 = 12.0 * solv_exp[0][variables[n1]] * solv_exp[1][variables[n1]] * solv_exp[2][variables[n1]]; // product : 0.5^3 = 0.125
2603
2604 tmp_paramsa1[n1] = paramsa1;
2605 tmp_paramsa2[n1] = paramsa2;
2606
2607 // here "cc" is a coupling constant used to stabilize the fit...
2608 // here "cc" is a coupling constant used to stabilize the fit...
2609 // here "cc" is a coupling constant used to stabilize the fit...
2610
2611 f64 cc = loop / 100.0; // this is for 250 steps!!!
2612 const f64 limit = 1.00; if (cc > limit) cc = limit;
2613
2614 dg += cc * wes * parames * 1000.0;
2615 dg += cc * sasa1 * paramsa1 * 1000.0;
2616 dg += cc * sasa2 * paramsa2 * 1000.0;
2617
2618 f64 k2 = exp(-dg / (8.31451 * 300.0));
2619 pKa = -log10(k2);
2620
2621 f64 tmp1 = pow(10.0, myprm->pH - pKa);
2622 f64 tmp2 = tmp1 / (1.0 + tmp1);
2623
2624 if (is_acid) tmp_charge[n1] = -tmp2;
2625 if (is_base) tmp_charge[n1] = 1.0 - tmp2;
2626
2627 tmp_newpKa[n1] = pKa; // do these oscillate as well???
2628 }
2629
2630 f64 delta = 0.0;
2631
2632 for (i32u n1 = 0;n1 < variables.size();n1++)
2633 {
2634 f64 ch1 = charge1[variables[n1]];
2635 f64 ch2 = tmp_charge[n1];
2636
2637 for (i32s n2 = 0;n2 < smoothing - 1;n2++)
2638 {
2639 f64 tmp1 = tmp_smooth[n1 * smoothing + n2 + 1];
2640 tmp_smooth[n1 * smoothing + n2 + 0] = tmp1;
2641 }
2642
2643 tmp_smooth[n1 * smoothing + (smoothing - 1)] = ch2;
2644
2645 ch2 = 0.0; // now do the smoothing by averaging the array...
2646 for (i32s n2 = 0;n2 < smoothing;n2++) ch2 += tmp_smooth[n1 * smoothing + n2];
2647 ch2 /= (f64) smoothing;
2648
2649 f64 dch = ch2 - ch1;
2650 delta += dch * dch;
2651
2652 charge1[variables[n1]] = ch2;
2653 //logf << ch2 << " ";
2654 }
2655 //logf << endl;
2656
2657 cout << "cycle = " << loop << " delta = " << delta << endl;
2658 // if ((loop % 50) == 49) { int qqq; cin >> qqq; }
2659 }
2660 //logf.close();
2661
2662 delete[] tmp_charge;
2663 delete[] tmp_smooth;
2664
2665 // do a simple charge neutralization...
2666
2667 f64 net_charge = 0.0; i32s nneg = 0; i32s npos = 0;
2668 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount();n1++)
2669 {
2670 charge2[n1] = charge1[n1];
2671
2672 net_charge += charge1[n1];
2673 if (charge1[n1] < -0.05) nneg++;
2674 if (charge1[n1] > +0.05) npos++;
2675 }
2676
2677 if (net_charge < -0.05)
2678 {
2679 cout << "NET NEGATIVE: " << net_charge << endl;
2680 f64 counterion_effect = fabs(net_charge) / (f64) nneg;
2681 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount();n1++)
2682 {
2683 if (charge1[n1] < -0.05) charge2[n1] += counterion_effect;
2684 }
2685 }
2686 else if (net_charge > +0.05)
2687 {
2688 cout << "NET POSITIVE: " << net_charge << endl;
2689 f64 counterion_effect = net_charge / (f64) npos;
2690 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount();n1++)
2691 {
2692 if (charge1[n1] > +0.05) charge2[n1] -= counterion_effect;
2693 }
2694 }
2695
2696 // report the results briefly...
2697
2698 i32u counter = 0;
2699 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount();n1++)
2700 {
2701 atmtab[n1]->charge = 0.0; // update the atom objects as well!!!
2702
2703 i32s cgi; i32s cvi;
2704 GetChgGrpVar(n1, cgi, cvi);
2705 if (cvi < 0) continue;
2706
2707 atmtab[n1]->charge = charge1[n1]; // update the atom objects as well!!! see StorePStatesToModel().
2708
2709 cout << "atmi = " << n1 << " (" << index_chn[n1] << "/" << index_res[n1] << ") ";
2710 cout << "pKa = " << tmp_newpKa[counter] << " (shift = " << (tmp_newpKa[counter] - myprm->charge_pKa[cvi]) << ") ";
2711 cout << "charge1 = " << charge1[n1] << " charge2 = " << charge2[n1] << endl;
2712
2713 counter++;
2714 }
2715
2716 // save the results in mdl::ref_civ if available; mdl::CheckProtonation() uses it!!!
2717 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2718
2719 setup1_sf * susf = dynamic_cast<setup1_sf *>(GetSetup());
2720 if (susf != NULL) susf->StorePStatesToModel(this);
2721
2722 delete[] tmp_vartab; tmp_vartab = NULL; // normally, this is not needed anymore...
2723 delete[] tmp_parames; tmp_parames = NULL; // normally, this is not needed anymore...
2724 delete[] tmp_paramsa1; tmp_paramsa1 = NULL; // normally, this is not needed anymore...
2725 delete[] tmp_paramsa2; tmp_paramsa2 = NULL; // normally, this is not needed anymore...
2726 delete[] tmp_newpKa; tmp_newpKa = NULL; // normally, this is not needed anymore...
2727
2728 // since the charges have changed, update sasa and neighbor lists as well!!!
2729 // since the charges have changed, update sasa and neighbor lists as well!!!
2730 // since the charges have changed, update sasa and neighbor lists as well!!!
2731
2732 setup1_sf * mysu = dynamic_cast<setup1_sf *>(GetSetup());
2733 if (!mysu) { cout << "BUG: cast to setup1_sf failed." << endl; exit(EXIT_FAILURE); }
2734
2735 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount() - num_solvent;n1++)
2736 {
2737 i32s indc = index_chn[n1];
2738 i32s indr = index_res[n1];
2739 f64 tmpc = charge1[n1];
2740
2741 char symbol = mysu->chn_vector[indc].res_vector[indr].symbol;
2742 i32s ind = 0; while (ind < 20 && model::sf_symbols[ind] != symbol) ind++;
2743 if (ind == 20) { cout << "BUG: unknown residue symbol." << endl; exit(EXIT_FAILURE); }
2744
2745 i32s inda = 0;
2746 while (mysu->chn_vector[indc].res_vector[indr].atmr[inda] != atmtab[n1])
2747 {
2748 inda++;
2749 if (inda >= mysu->chn_vector[indc].res_vector[indr].natm)
2750 {
2751 cout << "search of inda failed!!! i = " << n1 << endl;
2752 exit(EXIT_FAILURE);
2753 }
2754 }
2755
2756 i32s tmpt = (ind << 8) | inda;
2757 i32s tmpi = 0; while (tmpi < SF_NUM_TYPES && model::sf_types[tmpi] != tmpt) tmpi++;
2758 if (tmpi >= SF_NUM_TYPES) { cout << "BUG: tmpi overflow."; exit(EXIT_FAILURE); }
2759
2760 f64 tmps = (model::sf_is_polar[tmpi] ? myprm->imp_solv_1p : myprm->imp_solv_1n);
2761 tmps += myprm->imp_solv_2 * tmpc * tmpc;
2762
2763 for (i32u n2 = 0;n2 < LAYERS;n2++)
2764 {
2765 const f64 tmpsasa2[2] = { 0.80, 0.64 }; // NOT_PROPERLY_OPTIMIZED!!!
2766
2767 f64 svalue = tmps;
2768 if (n2 != 0)
2769 {
2770 f64 tmp1 = vdwr2[0][n1] / vdwr2[n2][n1];
2771 svalue *= tmp1 * tmpsasa2[n2 - 1];
2772
2773 // neg: -4 -2 -1 long range???
2774 // pos: +2 0 0 short range???
2775 if (svalue > 0.0) svalue = 0.0; // NOT_PROPERLY_OPTIMIZED!!!
2776 }
2777
2778 sasaE[n2][n1] = svalue * vdwr2[n2][n1];
2779
2780 // now update the neighbor lists...
2781
2782 if (nbt3_nl[n2][n1].index != NULL) delete[] nbt3_nl[n2][n1].index;
2783
2784 bool skip = false;
2785
2786 // do not skip any of the layers since solv_exp[] is used and should be up-to-date!!!
2787 // do not skip any of the layers since solv_exp[] is used and should be up-to-date!!!
2788 // do not skip any of the layers since solv_exp[] is used and should be up-to-date!!!
2789
2790 if (skip) /*nbt3_nl[n2][n1].index = NULL*/ exit(EXIT_FAILURE);
2791 else nbt3_nl[n2][n1].index = new i32s[size_nl[n2]];
2792 }
2793 }
2794 }
2795
GetChgGrpVar(i32s atmi,i32s & cgi,i32s & cvi)2796 void eng1_sf::GetChgGrpVar(i32s atmi, i32s & cgi, i32s & cvi)
2797 {
2798 i32s indc = index_chn[atmi];
2799 i32s indr = index_res[atmi];
2800
2801 if (indc < 0) // detect and handle solvent atom cases...
2802 {
2803 cgi = NOT_DEFINED;
2804 cvi = NOT_DEFINED;
2805 return;
2806 }
2807
2808 atom ** atmtab = GetSetup()->GetSFAtoms();
2809 setup1_sf * mysu = dynamic_cast<setup1_sf *>(GetSetup());
2810 if (!mysu) { cout << "BUG: cast to setup1_sf failed." << endl; exit(EXIT_FAILURE); }
2811 myprm = & mysu->prm;
2812
2813 i32s inda = 0;
2814 while (mysu->chn_vector[indc].res_vector[indr].atmr[inda] != atmtab[atmi])
2815 {
2816 inda++;
2817 if (inda >= mysu->chn_vector[indc].res_vector[indr].natm)
2818 {
2819 cout << "search of inda failed!!! i = " << atmi << endl;
2820 exit(EXIT_FAILURE);
2821 }
2822 }
2823
2824 // now determine the variable index!!!
2825 // now determine the variable index!!!
2826 // now determine the variable index!!!
2827
2828 cvi = NOT_DEFINED;
2829
2830 if (inda == 0 && indr == 0)
2831 {
2832 cvi = 0; // n-terminal
2833 }
2834
2835 if (inda == 0 && indr == ((i32s) mysu->chn_vector[indc].res_vector.size()) - 1)
2836 {
2837 cvi = 1; // c-terminal
2838 }
2839
2840 char symbol = mysu->chn_vector[indc].res_vector[indr].symbol;
2841 if (symbol == 'R' && inda == 2) cvi = 2;
2842 if (symbol == 'D' && inda == 1) cvi = 3;
2843 if (symbol == 'C' && inda == 1) cvi = 4; // cysteine!!!
2844 if (symbol == 'E' && inda == 1) cvi = 5;
2845 if (symbol == 'H' && inda == 1) cvi = 6;
2846 if (symbol == 'K' && inda == 2) cvi = 7;
2847 if (symbol == 'Y' && inda == 1) cvi = 8;
2848
2849 // check the cases where cysteine is disulphide-bridged!!!
2850 // check the cases where cysteine is disulphide-bridged!!!
2851 // check the cases where cysteine is disulphide-bridged!!!
2852
2853 if (cvi == 4) // this is a test for cysteine, see above...
2854 {
2855 bool flag = false;
2856
2857 for (i32u n1 = 0;n1 < mysu->dsb_vector.size();n1++)
2858 {
2859 if (mysu->dsb_vector[n1].chn[0] == indc && mysu->dsb_vector[n1].res[0] == indr) flag = true;
2860 if (mysu->dsb_vector[n1].chn[1] == indc && mysu->dsb_vector[n1].res[1] == indr) flag = true;
2861
2862 if (flag) break;
2863 }
2864
2865 if (flag) cvi = NOT_DEFINED;
2866 }
2867
2868 // now determine the group index!!!
2869 // now determine the group index!!!
2870 // now determine the group index!!!
2871
2872 switch (cvi)
2873 {
2874 case 1: // ct
2875 case 3: // D1
2876 case 5: // E1
2877 cgi = 0; // strong acids
2878 break;
2879
2880 case 4: // C1
2881 case 8: // Y1
2882 cgi = 1; // weak acids
2883 break;
2884
2885 case 0: // nt
2886 case 2: // R2
2887 case 7: // K2
2888 cgi = 2; // strong bases
2889 break;
2890
2891 case 6:
2892 cgi = 3; // weak bases
2893 break;
2894
2895 default:
2896 cgi = NOT_DEFINED;
2897 }
2898 }
2899
InitNBT1(sf_nbt1 * ref,vector<sf_nonbonded_lookup> & lookup)2900 bool eng1_sf::InitNBT1(sf_nbt1 * ref, vector<sf_nonbonded_lookup> & lookup)
2901 {
2902 atom ** atmtab = GetSetup()->GetSFAtoms();
2903
2904 setup1_sf * mysu = dynamic_cast<setup1_sf *>(GetSetup());
2905 if (!mysu) { cout << "BUG: cast to setup1_sf failed." << endl; exit(EXIT_FAILURE); }
2906 myprm = & mysu->prm;
2907
2908 // these are used to save memory when measuring vdw-radii and surface areas... 2.0 * (0.3 + 0.2) = 1.0
2909 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2910 //const fGL * c1 = atmtab[ref->atmi[0]]->GetCRD(0); const fGL * c2 = atmtab[ref->atmi[1]]->GetCRD(0);
2911 //v3d<fGL> v1 = v3d<fGL>(c1, c2); if (v1.len() > 1.5) return false;
2912
2913 i32s indc1 = index_chn[ref->atmi[0]];
2914 i32s indr1 = index_res[ref->atmi[0]];
2915
2916 i32s indc2 = index_chn[ref->atmi[1]];
2917 i32s indr2 = index_res[ref->atmi[1]];
2918
2919 if (indc1 < 0 || indc2 < 0) // solvent-solute or solvent-solvent
2920 {
2921 f64 radius[2];
2922 radius[0] = vdwr[ref->atmi[0]];
2923 radius[1] = vdwr[ref->atmi[1]];
2924
2925 f64 optr = radius[0] + radius[1];
2926 f64 opte;
2927
2928 // nonbonded.txt min/median/max = 0.18/0.83/3.56
2929 // that is multiplied by myprm->lenjon = 1.225 -> average = 0.83 * 1.225 = 1.02
2930
2931 // for H2O, if opte > 3.00 -> more or less frozen state?!?!?!
2932 // difficult to estimate for other cases, since also depends on mass/size...
2933
2934 if (indc1 < 0 && indc2 < 0) opte = 2.00; // opte > 3.00 -> frozen???
2935 else // solvent-solute cases...
2936 {
2937 i32s indc; i32s indr; i32s n1;
2938 if (indc2 < 0) { indc = indc1; indr = indr1; n1 = ref->atmi[0]; }
2939 else { indc = indc2; indr = indr2; n1 = ref->atmi[1]; }
2940
2941 f64 tmpc = charge1[n1];
2942
2943 char symbol = mysu->chn_vector[indc].res_vector[indr].symbol;
2944 i32s ind = 0; while (ind < 20 && model::sf_symbols[ind] != symbol) ind++;
2945 if (ind == 20) { cout << "BUG: unknown residue symbol." << endl; exit(EXIT_FAILURE); }
2946
2947 i32s inda = 0;
2948 while (mysu->chn_vector[indc].res_vector[indr].atmr[inda] != atmtab[n1])
2949 {
2950 inda++;
2951 if (inda >= mysu->chn_vector[indc].res_vector[indr].natm)
2952 {
2953 cout << "search of inda failed!!! i = " << n1 << endl;
2954 exit(EXIT_FAILURE);
2955 }
2956 }
2957
2958 i32s tmpt = (ind << 8) | inda;
2959 i32s tmpi = 0; while (tmpi < SF_NUM_TYPES && model::sf_types[tmpi] != tmpt) tmpi++;
2960 if (tmpi >= SF_NUM_TYPES) { cout << "BUG: tmpi overflow."; exit(EXIT_FAILURE); }
2961
2962 f64 tmps = (model::sf_is_polar[tmpi] ? myprm->exp_solv_1p : myprm->exp_solv_1n);
2963 tmps += myprm->exp_solv_2 * tmpc * tmpc;
2964
2965 // ??? see also SetupCharges()!!!
2966
2967 opte = tmps;
2968 }
2969
2970 InitLenJon(ref, optr, opte);
2971 return true;
2972 }
2973
2974 i32s inda1 = 0;
2975 while (mysu->chn_vector[indc1].res_vector[indr1].atmr[inda1] != atmtab[ref->atmi[0]])
2976 {
2977 inda1++;
2978 if (inda1 >= mysu->chn_vector[indc1].res_vector[indr1].natm)
2979 {
2980 cout << "initNB : search of inda1 failed!!! i = " << ref->atmi[0] << endl;
2981 exit(EXIT_FAILURE);
2982 }
2983 }
2984
2985 i32s inda2 = 0;
2986 while (mysu->chn_vector[indc2].res_vector[indr2].atmr[inda2] != atmtab[ref->atmi[1]])
2987 {
2988 inda2++;
2989 if (inda2 >= mysu->chn_vector[indc2].res_vector[indr2].natm)
2990 {
2991 cout << "initNB : search of inda2 failed!!! i = " << ref->atmi[1] << endl;
2992 exit(EXIT_FAILURE);
2993 }
2994 }
2995
2996 // nbt1 are needed for all VA's separated by more than 1 bonds (solvent term).
2997 // nbt1 are needed for all VA's separated by more than 1 bonds (solvent term).
2998 // nbt1 are needed for all VA's separated by more than 1 bonds (solvent term).
2999
3000 if (inda1 && inda2) // bridge???
3001 {
3002 for (i32u n1 = 0;n1 < mysu->dsb_vector.size();n1++)
3003 {
3004 bool test1a = (indc1 == mysu->dsb_vector[n1].chn[0]);
3005 bool test1b = (indr1 == mysu->dsb_vector[n1].res[0]);
3006 bool test1c = (indc2 == mysu->dsb_vector[n1].chn[1]);
3007 bool test1d = (indr2 == mysu->dsb_vector[n1].res[1]);
3008 if (test1a && test1b && test1c && test1d) return false;
3009
3010 bool test2a = (indc1 == mysu->dsb_vector[n1].chn[1]);
3011 bool test2b = (indr1 == mysu->dsb_vector[n1].res[1]);
3012 bool test2c = (indc2 == mysu->dsb_vector[n1].chn[0]);
3013 bool test2d = (indr2 == mysu->dsb_vector[n1].res[0]);
3014 if (test2a && test2b && test2c && test2d) return false;
3015 }
3016 }
3017
3018 bool same_chn = (indc1 == indc2);
3019
3020 bool rdist_is_0 = same_chn && ((indr2 - indr1) == 0);
3021 bool rdist_is_1 = same_chn && ((indr2 - indr1) == 1);
3022
3023 if (rdist_is_0 && ((inda2 - inda1) < 2)) return false;
3024 if (rdist_is_1 && !inda1 && !inda2) return false;
3025
3026 bool rdist_is_2 = same_chn && ((indr2 - indr1) == 2);
3027 bool rdist_is_3 = same_chn && ((indr2 - indr1) == 3);
3028 bool rdist_is_4 = same_chn && ((indr2 - indr1) == 4);
3029
3030 bool both_in_main_chn = (!inda1 && !inda2);
3031
3032 bool check3 = (indr1 + 1 < (i32s) mysu->chn_vector[indc1].res_vector.size());
3033 bool check4 = (indr1 + 2 < (i32s) mysu->chn_vector[indc1].res_vector.size());
3034 bool check5 = (indr1 + 3 < (i32s) mysu->chn_vector[indc1].res_vector.size());
3035
3036 bool helix3 = check3 && (mysu->chn_vector[indc1].res_vector[indr1 + 1].state == SF_STATE_HELIX4);
3037 bool helix4 = check4 && helix3 && (mysu->chn_vector[indc1].res_vector[indr1 + 2].state == SF_STATE_HELIX4);
3038 bool helix5 = check5 && helix4 && (mysu->chn_vector[indc1].res_vector[indr1 + 3].state == SF_STATE_HELIX4);
3039
3040 f64 radius[2];
3041 radius[0] = vdwr[ref->atmi[0]];
3042 radius[1] = vdwr[ref->atmi[1]];
3043
3044 f64 optr = radius[0] + radius[1];
3045 f64 opte = 0.0;
3046
3047 char symbol[2] =
3048 {
3049 mysu->chn_vector[indc1].res_vector[indr1].symbol,
3050 mysu->chn_vector[indc2].res_vector[indr2].symbol
3051 };
3052
3053 i32u index = 0;
3054 while (true)
3055 {
3056 if (index >= lookup.size()) { cout << "nonbonded lookup failed!!!" << endl; exit(EXIT_FAILURE); }
3057
3058 bool test1 = (lookup[index].s1 == symbol[0] && lookup[index].a1 == (i32u) inda1 && lookup[index].s2 == symbol[1] && lookup[index].a2 == (i32u) inda2);
3059 bool test2 = (lookup[index].s1 == symbol[1] && lookup[index].a1 == (i32u) inda2 && lookup[index].s2 == symbol[0] && lookup[index].a2 == (i32u) inda1);
3060 if (test1 || test2)
3061 {
3062 opte = lookup[index].opte * myprm->lenjon;
3063 break; // exit the loop...
3064 }
3065 else index++;
3066 }
3067
3068 // use modified values in helices!!! these will produce correct helix dimensions.
3069 // use modified values in helices!!! these will produce correct helix dimensions.
3070 // use modified values in helices!!! these will produce correct helix dimensions.
3071
3072 if (both_in_main_chn)
3073 {
3074 bool test3 = (rdist_is_2 && helix3);
3075 bool test4 = (rdist_is_3 && helix4);
3076 bool test5 = (rdist_is_4 && helix5);
3077
3078 if (test3 || test4 || test5)
3079 {
3080 optr = 0.650;
3081 opte = 0.95 * myprm->lenjon;
3082 }
3083 }
3084
3085 InitLenJon(ref, optr, opte);
3086 return true;
3087 }
3088
InitLenJon(sf_nbt1 * ref,f64 opt,f64 fc)3089 void eng1_sf::InitLenJon(sf_nbt1 * ref, f64 opt, f64 fc)
3090 {
3091 if (opt < 0.100)
3092 {
3093 cout << "eng1_sf::InitLenJon() : too small opt : " << opt << endl;
3094 exit(EXIT_FAILURE);
3095 }
3096
3097 if (fc < 0.100)
3098 {
3099 cout << "eng1_sf::InitLenJon() : too small fc : " << fc << endl;
3100 exit(EXIT_FAILURE);
3101 }
3102
3103 // ref->data[0] = opt * pow(2.0 * fc, 1.0 / 9.0); // 6-9
3104 // ref->data[1] = opt * pow(3.0 * fc, 1.0 / 6.0);
3105 ref->data[0] = opt * pow(1.0 * fc, 1.0 / 12.0); // 6-12
3106 ref->data[1] = opt * pow(2.0 * fc, 1.0 / 6.0);
3107 // ref->data[0] = opt * pow(3.0 * fc, 1.0 / 12.0); // 9-12
3108 // ref->data[1] = opt * pow(4.0 * fc, 1.0 / 9.0);
3109 }
3110
ComputeBT1(i32u p1)3111 void eng1_sf::ComputeBT1(i32u p1)
3112 {
3113 energy_bt1 = 0.0;
3114
3115 // data1 -> length of the bond vector, in nanometers [nm].
3116 // data2[0] -> grad[0-2]: for atom 1 when direction = 0, for atom 0 when direction = 1
3117 // data2[1] -> grad[0-2]: for atom 0 when direction = 0, for atom 1 when direction = 1
3118
3119 // mm1_eng_exp1::ComputeBT1() works in a similar way...
3120
3121 for (i32u n1 = 0;n1 < bt1_vector.size();n1++)
3122 {
3123 i32s * atmi = bt1_vector[n1].atmi;
3124
3125 f64 t1a[3]; f64 t1b = 0.0;
3126 for (i32s n2 = 0;n2 < 3;n2++)
3127 {
3128 f64 t9a = crd[l2g_sf[atmi[0]] * 3 + n2];
3129 f64 t9b = crd[l2g_sf[atmi[1]] * 3 + n2];
3130
3131 t1a[n2] = t9a - t9b;
3132 t1b += t1a[n2] * t1a[n2];
3133 }
3134
3135 f64 t1c = sqrt(t1b);
3136 bt1data[n1].data1 = t1c;
3137 for (i32s n2 = 0;n2 < 3;n2++)
3138 {
3139 f64 t9a = t1a[n2] / t1c;
3140 bt1data[n1].data2[0][n2] = +t9a;
3141 bt1data[n1].data2[1][n2] = -t9a;
3142 }
3143
3144 /*##############################################*/
3145 /*##############################################*/
3146 // this is for the surface area term...
3147
3148 bool first = (atmi[0] > atmi[1]);
3149 dist2[dist1[atmi[first]] + (atmi[!first] - atmi[first]) - 1] = t1c;
3150
3151 // 1st layer...
3152 // 1st layer...
3153 // 1st layer...
3154
3155 // none of the SASA terms should be skipped at 1st layer -> no test...
3156
3157 if (t1c < (vdwr1[0][atmi[0]] + vdwr1[0][atmi[1]]))
3158 {
3159 nbt3_nl[0][atmi[0]].index[nbt3_nl[0][atmi[0]].index_count++] = atmi[1];
3160 if (nbt3_nl[0][atmi[0]].index_count >= size_nl[0]) { cout << "BUG: NL overflow 1a!!!" << endl; exit(EXIT_FAILURE); }
3161
3162 nbt3_nl[0][atmi[1]].index[nbt3_nl[0][atmi[1]].index_count++] = atmi[0];
3163 if (nbt3_nl[0][atmi[1]].index_count >= size_nl[0]) { cout << "BUG: NL overflow 1a!!!" << endl; exit(EXIT_FAILURE); }
3164 }
3165
3166 // 2nd layer...
3167 // 2nd layer...
3168 // 2nd layer...
3169
3170 bool test2a = (nbt3_nl[1][atmi[0]].index != NULL);
3171 if (test2a && t1c < (vdwr1[1][atmi[0]] + vdwr1[0][atmi[1]]) && t1c > (vdwr1[1][atmi[0]] - vdwr1[0][atmi[1]]))
3172 {
3173 nbt3_nl[1][atmi[0]].index[nbt3_nl[1][atmi[0]].index_count++] = atmi[1];
3174 if (nbt3_nl[1][atmi[0]].index_count >= size_nl[1]) { cout << "BUG: NL overflow 2a!!!" << endl; exit(EXIT_FAILURE); }
3175 }
3176
3177 bool test2b = (nbt3_nl[1][atmi[1]].index != NULL);
3178 if (test2b && t1c < (vdwr1[0][atmi[0]] + vdwr1[1][atmi[1]]) && t1c > (vdwr1[1][atmi[1]] - vdwr1[0][atmi[0]]))
3179 {
3180 nbt3_nl[1][atmi[1]].index[nbt3_nl[1][atmi[1]].index_count++] = atmi[0];
3181 if (nbt3_nl[1][atmi[1]].index_count >= size_nl[1]) { cout << "BUG: NL overflow 2a!!!" << endl; exit(EXIT_FAILURE); }
3182 }
3183
3184 // 3rd layer...
3185 // 3rd layer...
3186 // 3rd layer...
3187
3188 bool test3a = (nbt3_nl[2][atmi[0]].index != NULL);
3189 if (test3a && t1c < (vdwr1[2][atmi[0]] + vdwr1[0][atmi[1]]) && t1c > (vdwr1[2][atmi[0]] - vdwr1[0][atmi[1]]))
3190 {
3191 nbt3_nl[2][atmi[0]].index[nbt3_nl[2][atmi[0]].index_count++] = atmi[1];
3192 if (nbt3_nl[2][atmi[0]].index_count >= size_nl[2]) { cout << "BUG: NL overflow 3a!!!" << endl; exit(EXIT_FAILURE); }
3193 }
3194
3195 bool test3b = (nbt3_nl[2][atmi[1]].index != NULL);
3196 if (test3b && t1c < (vdwr1[0][atmi[0]] + vdwr1[2][atmi[1]]) && t1c > (vdwr1[2][atmi[1]] - vdwr1[0][atmi[0]]))
3197 {
3198 nbt3_nl[2][atmi[1]].index[nbt3_nl[2][atmi[1]].index_count++] = atmi[0];
3199 if (nbt3_nl[2][atmi[1]].index_count >= size_nl[2]) { cout << "BUG: NL overflow 3a!!!" << endl; exit(EXIT_FAILURE); }
3200 }
3201
3202 // this is for the surface area term...
3203 /*##############################################*/
3204 /*##############################################*/
3205
3206 // f = a(x-b)^2
3207 // df/dx = 2a(x-b)
3208
3209 f64 t2a = t1c - bt1_vector[n1].opt;
3210 energy_bt1 += bt1_vector[n1].fc * t2a * t2a;
3211
3212 if (p1 > 0)
3213 {
3214 f64 t2b = 2.0 * bt1_vector[n1].fc * t2a;
3215 for (i32s n2 = 0;n2 < 3;n2++)
3216 {
3217 f64 t2c = bt1data[n1].data2[0][n2] * t2b;
3218
3219 d1[l2g_sf[atmi[0]] * 3 + n2] += t2c;
3220 d1[l2g_sf[atmi[1]] * 3 + n2] -= t2c;
3221 }
3222 }
3223 }
3224 }
3225
ComputeBT2(i32u p1)3226 void eng1_sf::ComputeBT2(i32u p1)
3227 {
3228 energy_bt2 = 0.0;
3229
3230 // data1 -> cosine of the bond angle, in the usual range [-1.0, +1.0]
3231 // data2[0] -> grad[0-2]: for atom 2 when direction = 0, for atom 0 when direction = 1
3232 // data2[1] -> grad[0-2]: for atom 1 when direction = 0, for atom 1 when direction = 1
3233 // data2[2] -> grad[0-2]: for atom 0 when direction = 0, for atom 2 when direction = 1
3234
3235 // mm1_eng_exp1::ComputeBT2() works in a similar way...
3236
3237 for (i32u n1 = 0;n1 < bt2_vector.size();n1++)
3238 {
3239 i32s * atmi = bt2_vector[n1].atmi;
3240
3241 i32s * index1 = bt2_vector[n1].index1;
3242 bool * dir1 = bt2_vector[n1].dir1;
3243
3244 f64 * t1a = bt1data[index1[0]].data2[dir1[0]];
3245 f64 * t1b = bt1data[index1[1]].data2[dir1[1]];
3246
3247 f64 t1c = t1a[0] * t1b[0] + t1a[1] * t1b[1] + t1a[2] * t1b[2];
3248
3249 if (t1c < -1.0) t1c = -1.0; // domain check...
3250 if (t1c > +1.0) t1c = +1.0; // domain check...
3251
3252 // if the main-chain angles approach 180.0 deg we will have some serious numerical problems in later
3253 // stages. fc[1] must be large enough to prevent such problems. here is how we can monitor this...
3254
3255 bool problem1 = (t1c < -0.999); // is it too close to 180 deg ???
3256 bool problem2 = (bt2_vector[n1].fc[1] > 0.0); // is it really a main-chain angle???
3257 if (problem1 && problem2) { cout << "BUG: BT2 ang -> 180.0 deg." << endl; exit(EXIT_FAILURE); }
3258
3259 bt2data[n1].data1 = t1c;
3260
3261 for (i32s n2 = 0;n2 < 3;n2++)
3262 {
3263 f64 t9a = (t1b[n2] - t1c * t1a[n2]) / bt1data[index1[0]].data1;
3264 f64 t9b = (t1a[n2] - t1c * t1b[n2]) / bt1data[index1[1]].data1;
3265
3266 bt2data[n1].data2[0][n2] = t9a;
3267 bt2data[n1].data2[1][n2] = -(t9a + t9b);
3268 bt2data[n1].data2[2][n2] = t9b;
3269 }
3270
3271 // f = a(x-b)^2
3272 // df/dx = 2a(x-b)
3273
3274 f64 t3b = t1c - bt2_vector[n1].opt;
3275 energy_bt2 += bt2_vector[n1].fc[0] * t3b * t3b;
3276
3277 // f = a/(x+1)^2
3278 // df/dx = -2a/(x+1)^3
3279
3280 f64 t3c = t1c + 1.0;
3281
3282 f64 t3d = t3c * t3c;
3283 energy_bt2 += bt2_vector[n1].fc[1] / t3d;
3284
3285 if (p1 > 0)
3286 {
3287 f64 t2a = 2.0 * bt2_vector[n1].fc[0] * t3b;
3288
3289 f64 t3f = t3d * t3c;
3290 f64 t2b = 2.0 * bt2_vector[n1].fc[1] / t3f;
3291
3292 f64 t2c = t2a - t2b;
3293
3294 for (i32s n2 = 0;n2 < 3;n2++)
3295 {
3296 d1[l2g_sf[atmi[0]] * 3 + n2] += bt2data[n1].data2[0][n2] * t2c;
3297 d1[l2g_sf[atmi[1]] * 3 + n2] += bt2data[n1].data2[1][n2] * t2c;
3298 d1[l2g_sf[atmi[2]] * 3 + n2] += bt2data[n1].data2[2][n2] * t2c;
3299 }
3300 }
3301 }
3302 }
3303
ComputeBT3(i32u p1)3304 void eng1_sf::ComputeBT3(i32u p1)
3305 {
3306 energy_bt3a = 0.0;
3307 energy_bt3b = 0.0;
3308
3309 for (i32u n1 = 0;n1 < bt3_vector.size();n1++)
3310 {
3311 i32s * atmi = bt3_vector[n1].atmi;
3312
3313 i32s * index2 = bt3_vector[n1].index2;
3314 i32s * index1 = bt3_vector[n1].index1;
3315 bool * dir1 = bt3_vector[n1].dir1;
3316
3317 // ang is never > 180.0 deg -> sin(ang) is never < 0.0 -> no need to check signs!!!
3318 // ang is never > 180.0 deg -> sin(ang) is never < 0.0 -> no need to check signs!!!
3319 // ang is never > 180.0 deg -> sin(ang) is never < 0.0 -> no need to check signs!!!
3320
3321 f64 t1a[2] = { bt2data[index2[0]].data1, bt2data[index2[1]].data1 };
3322 f64 t1b[2] = { 1.0 - t1a[0] * t1a[0], 1.0 - t1a[1] * t1a[1] };
3323
3324 f64 t1c[2][3];
3325 t1c[0][0] = bt1data[index1[0]].data2[dir1[0]][0] - t1a[0] * bt1data[index1[1]].data2[dir1[1]][0];
3326 t1c[0][1] = bt1data[index1[0]].data2[dir1[0]][1] - t1a[0] * bt1data[index1[1]].data2[dir1[1]][1];
3327 t1c[0][2] = bt1data[index1[0]].data2[dir1[0]][2] - t1a[0] * bt1data[index1[1]].data2[dir1[1]][2];
3328 t1c[1][0] = bt1data[index1[3]].data2[dir1[3]][0] - t1a[1] * bt1data[index1[2]].data2[dir1[2]][0];
3329 t1c[1][1] = bt1data[index1[3]].data2[dir1[3]][1] - t1a[1] * bt1data[index1[2]].data2[dir1[2]][1];
3330 t1c[1][2] = bt1data[index1[3]].data2[dir1[3]][2] - t1a[1] * bt1data[index1[2]].data2[dir1[2]][2];
3331
3332 f64 t9a = t1c[0][0] * t1c[1][0] + t1c[0][1] * t1c[1][1] + t1c[0][2] * t1c[1][2];
3333 f64 t1d = t9a / sqrt(t1b[0] * t1b[1]);
3334
3335 if (t1d < -1.0) t1d = -1.0; // domain check...
3336 if (t1d > +1.0) t1d = +1.0; // domain check...
3337
3338 f64 t1e = acos(t1d);
3339
3340 // now we still have to determine sign of the result...
3341
3342 v3d<f64> tmpv1 = v3d<f64>(t1c[0]);
3343 v3d<f64> tmpv2 = v3d<f64>(bt1data[index1[2]].data2[dir1[2]]);
3344 v3d<f64> tmpv3 = v3d<f64>(bt1data[index1[3]].data2[dir1[3]]);
3345 v3d<f64> tmpv4 = tmpv2.vpr(tmpv3);
3346
3347 if (tmpv1.spr(tmpv4) < 0.0) t1e = -t1e;
3348
3349 f64 t1f; f64 t1g; // f and df/dx !!!
3350 f64 cs[3]; f64 sn[3];
3351
3352 if (bt3_vector[n1].tor_ttype != TTYPE_LOOP) // helix + strand
3353 {
3354 // Dx = x-b | for the following formulas...
3355
3356 // f = a(2PI-Dx)^2 | if Dx > +PI
3357 // df/fx = -2a(2PI-Dx)
3358
3359 // f = a(2PI+Dx)^2 | if Dx < -PI
3360 // df/fx = +2a(2PI+Dx)
3361
3362 // f = a(Dx)^2 | otherwise
3363 // df/fx = 2a(Dx)
3364
3365 f64 t1h = t1e - bt3_vector[n1].tors[0];
3366 if (t1h > +M_PI)
3367 {
3368 t1h = 2.0 * M_PI - t1h;
3369
3370 t1f = bt3_vector[n1].tors[1] * t1h * t1h;
3371 t1g = -2.0 * bt3_vector[n1].tors[1] * t1h;
3372 }
3373 else if (t1h < -M_PI)
3374 {
3375 t1h = 2.0 * M_PI + t1h;
3376
3377 t1f = bt3_vector[n1].tors[1] * t1h * t1h;
3378 t1g = +2.0 * bt3_vector[n1].tors[1] * t1h;
3379 }
3380 else
3381 {
3382 t1f = bt3_vector[n1].tors[1] * t1h * t1h;
3383 t1g = 2.0 * bt3_vector[n1].tors[1] * t1h;
3384 }
3385
3386 energy_bt3a += t1f;
3387 }
3388 else // loop!!!
3389 {
3390 // f = a*cos(x) + b*cos(2x) + c*cos(3x) + d*sin(x) + e*sin(2x) + f*sin(3x)
3391 // df/dx = d*cos(x) + 2e*cos(2x) + 3f*cos(3x) - a*sin(x) - 2b*sin(2x) - 3c*sin(3x)
3392
3393 f64 t1m = t1e + t1e; f64 t1n = t1m + t1e;
3394 cs[0] = cos(t1e); cs[1] = cos(t1m); cs[2] = cos(t1n);
3395 sn[0] = sin(t1e); sn[1] = sin(t1m); sn[2] = sin(t1n);
3396
3397 f64 * fsc = bt3_vector[n1].torc; f64 * fss = bt3_vector[n1].tors;
3398
3399 t1f = fsc[0] * cs[0] + fsc[1] * cs[1] + fsc[2] * cs[2];
3400 t1f += fss[0] * sn[0] + fss[1] * sn[1] + fss[2] * sn[2];
3401
3402 t1g = fss[0] * cs[0] + 2.0 * fss[1] * cs[1] + 3.0 * fss[2] * cs[2];
3403 t1g -= fsc[0] * sn[0] + 2.0 * fsc[1] * sn[1] + 3.0 * fsc[2] * sn[2];
3404
3405 energy_bt3b += t1f;
3406 }
3407
3408 f64 t4c[3]; f64 t5c[3]; f64 t6a[3]; f64 t7a[3]; // these are needed later...
3409
3410 if (p1 > 0)
3411 {
3412 f64 t2a = bt1data[index1[0]].data1 * t1b[0];
3413 f64 t2b = bt1data[index1[0]].data1 * t1a[0] / bt1data[index1[1]].data1;
3414
3415 f64 t3a = bt1data[index1[3]].data1 * t1b[1];
3416 f64 t3b = bt1data[index1[3]].data1 * t1a[1] / bt1data[index1[2]].data1;
3417
3418 const i32s cp[3][3] = { { 0, 1, 2 }, { 1, 2, 0 }, { 2, 0, 1 } };
3419
3420 for (i32s n2 = 0;n2 < 3;n2++)
3421 {
3422 f64 t4a = bt1data[index1[0]].data2[dir1[0]][cp[n2][1]] * bt1data[index1[1]].data2[dir1[1]][cp[n2][2]];
3423 f64 t4b = bt1data[index1[0]].data2[dir1[0]][cp[n2][2]] * bt1data[index1[1]].data2[dir1[1]][cp[n2][1]];
3424 t4c[n2] = (t4a - t4b) / t2a;
3425
3426 f64 t5a = bt1data[index1[2]].data2[dir1[2]][cp[n2][2]] * bt1data[index1[3]].data2[dir1[3]][cp[n2][1]];
3427 f64 t5b = bt1data[index1[2]].data2[dir1[2]][cp[n2][1]] * bt1data[index1[3]].data2[dir1[3]][cp[n2][2]];
3428 t5c[n2] = (t5a - t5b) / t3a;
3429
3430 d1[l2g_sf[atmi[0]] * 3 + n2] += t1g * t4c[n2];
3431 d1[l2g_sf[atmi[3]] * 3 + n2] += t1g * t5c[n2];
3432
3433 t6a[n2] = (t2b - 1.0) * t4c[n2] - t3b * t5c[n2];
3434 t7a[n2] = (t3b - 1.0) * t5c[n2] - t2b * t4c[n2];
3435
3436 d1[l2g_sf[atmi[1]] * 3 + n2] += t1g * t6a[n2];
3437 d1[l2g_sf[atmi[2]] * 3 + n2] += t1g * t7a[n2];
3438 }
3439 }
3440
3441 // calculate dipole directions (with derivatives if needed)...
3442 // calculate dipole directions (with derivatives if needed)...
3443 // calculate dipole directions (with derivatives if needed)...
3444
3445 // f = a*cos(x) + b*cos(2x) + c*cos(3x) + d*sin(x) + e*sin(2x) + f*sin(3x) + g*x + h
3446 // df/dx = d*cos(x) + 2e*cos(2x) + 3f*cos(3x) - a*sin(x) - 2b*sin(2x) - 3c*sin(3x) + g
3447
3448 f64 loop1; f64 loop2; f64 dpbdd[4][3];
3449 f64 * fsc = bt3_vector[n1].dipc; f64 * fss = bt3_vector[n1].dips;
3450
3451 switch (bt3_vector[n1].dip_ttype)
3452 {
3453 case TTYPE_HELIX:
3454 bt3_vector[n1].pbdd = -0.935131;
3455 break;
3456
3457 case TTYPE_STRAND:
3458 bt3_vector[n1].pbdd = +1.994676;
3459 break;
3460
3461 default: // this is for TTYPE_LOOP...
3462 loop1 = bt3_vector[n1].dipk[0] * t1e + bt3_vector[n1].dipk[1];
3463 loop1 += fsc[0] * cs[0] + fsc[1] * cs[1] + fsc[2] * cs[2];
3464 loop1 += fss[0] * sn[0] + fss[1] * sn[1] + fss[2] * sn[2];
3465 bt3_vector[n1].pbdd = loop1;
3466
3467 if (p1 > 0)
3468 {
3469 loop2 = bt3_vector[n1].dipk[0];
3470 loop2 += fss[0] * cs[0] + 2.0 * fss[1] * cs[1] + 3.0 * fss[2] * cs[2];
3471 loop2 -= fsc[0] * sn[0] + 2.0 * fsc[1] * sn[1] + 3.0 * fsc[2] * sn[2];
3472
3473 for (i32s n2 = 0;n2 < 3;n2++)
3474 {
3475 dpbdd[0][n2] = loop2 * t4c[n2];
3476 dpbdd[3][n2] = loop2 * t5c[n2];
3477
3478 dpbdd[1][n2] = loop2 * t6a[n2];
3479 dpbdd[2][n2] = loop2 * t7a[n2];
3480 }
3481 }
3482 }
3483
3484 // calculate some dipole results already here...
3485 // calculate some dipole results already here...
3486 // calculate some dipole results already here...
3487
3488 if (bt3_vector[n1].skip) continue; // skip all X-pro cases...
3489
3490 f64 t2a[3]; f64 t2b[3][3]; // 1st unit bond vector + derivatives...
3491 f64 t2c[3]; f64 t2d[3][3]; // 2nd unit bond vector + derivatives...
3492
3493 for (i32s n2 = 0;n2 < 3;n2++)
3494 {
3495 t2a[n2] = bt1data[index1[0]].data2[dir1[0]][n2];
3496 t2c[n2] = bt1data[index1[1]].data2[dir1[1]][n2];
3497 }
3498
3499 if (p1 > 0)
3500 {
3501 for (i32s n2 = 0;n2 < 3;n2++)
3502 {
3503 for (i32s n3 = 0;n3 < 3;n3++)
3504 {
3505 if (n2 != n3)
3506 {
3507 t2b[n2][n3] = -t2a[n2] * t2a[n3] / bt1data[index1[0]].data1;
3508 t2d[n2][n3] = -t2c[n2] * t2c[n3] / bt1data[index1[1]].data1;
3509 }
3510 else
3511 {
3512 t2b[n2][n2] = (1.0 - t2a[n2] * t2a[n2]) / bt1data[index1[0]].data1;
3513 t2d[n2][n2] = (1.0 - t2c[n2] * t2c[n2]) / bt1data[index1[1]].data1;
3514 }
3515 }
3516 }
3517 }
3518
3519 f64 t2e = bt2data[index2[0]].data1;
3520 f64 t2f = t2e * t2e; f64 t2g = 1.0 - t2f;
3521 f64 t2h = sqrt(t2g);
3522
3523 // 1st basis vector... 1st basis vector... 1st basis vector... 1st basis vector...
3524 // 1st basis vector... 1st basis vector... 1st basis vector... 1st basis vector...
3525 // 1st basis vector... 1st basis vector... 1st basis vector... 1st basis vector...
3526
3527 bt3_vector[n1].bv[0][0] = (t2a[0] - t2e * t2c[0]) / t2h;
3528 bt3_vector[n1].bv[0][1] = (t2a[1] - t2e * t2c[1]) / t2h;
3529 bt3_vector[n1].bv[0][2] = (t2a[2] - t2e * t2c[2]) / t2h;
3530
3531 if (p1 > 0)
3532 {
3533 for (i32s n2 = 0;n2 < 3;n2++)
3534 {
3535 f64 t3a = -t2e * bt2data[index2[0]].data2[0][n2] / t2h; // D sin(alpha)
3536 f64 t3b = -t2e * bt2data[index2[0]].data2[2][n2] / t2h; // D sin(alpha)
3537
3538 for (i32s n3 = 0;n3 < 3;n3++)
3539 {
3540 f64 t4a = t2h * t2b[n2][n3] - t2a[n3] * t3a;
3541 f64 t4b = t2h * bt2data[index2[0]].data2[0][n2] - t2e * t3a;
3542 f64 t4c = (t4a - t2c[n3] * t4b) / t2g;
3543
3544 f64 t5a = t3b * (t2a[n3] - t2e * t2c[n3]);
3545 f64 t5b = t2h * (t2e * t2d[n2][n3] + t2c[n3] * bt2data[index2[0]].data2[2][n2]);
3546 f64 t5c = -(t5a + t5b) / t2g;
3547
3548 bt3_vector[n1].dbv[0][0][n2][n3] = t4c;
3549 bt3_vector[n1].dbv[0][1][n2][n3] = -(t4c + t5c);
3550 bt3_vector[n1].dbv[0][2][n2][n3] = t5c;
3551 }
3552 }
3553 }
3554
3555 // 2nd basis vector... 2nd basis vector... 2nd basis vector... 2nd basis vector...
3556 // 2nd basis vector... 2nd basis vector... 2nd basis vector... 2nd basis vector...
3557 // 2nd basis vector... 2nd basis vector... 2nd basis vector... 2nd basis vector...
3558
3559 f64 t2i[3];
3560 t2i[0] = bt1data[index1[0]].data2[dir1[0]][1] * bt1data[index1[1]].data2[dir1[1]][2] -
3561 bt1data[index1[0]].data2[dir1[0]][2] * bt1data[index1[1]].data2[dir1[1]][1];
3562 t2i[1] = bt1data[index1[0]].data2[dir1[0]][2] * bt1data[index1[1]].data2[dir1[1]][0] -
3563 bt1data[index1[0]].data2[dir1[0]][0] * bt1data[index1[1]].data2[dir1[1]][2];
3564 t2i[2] = bt1data[index1[0]].data2[dir1[0]][0] * bt1data[index1[1]].data2[dir1[1]][1] -
3565 bt1data[index1[0]].data2[dir1[0]][1] * bt1data[index1[1]].data2[dir1[1]][0];
3566
3567 bt3_vector[n1].bv[1][0] = t2i[0] / t2h;
3568 bt3_vector[n1].bv[1][1] = t2i[1] / t2h;
3569 bt3_vector[n1].bv[1][2] = t2i[2] / t2h;
3570
3571 if (p1 > 0)
3572 {
3573 for (i32s n2 = 0;n2 < 3;n2++)
3574 {
3575 f64 t3a[2]; // ???epsilon i
3576 t3a[0] = bt2data[index2[0]].data2[0][n2] * bt2data[index2[0]].data1 / t2g;
3577 t3a[1] = bt2data[index2[0]].data2[2][n2] * bt2data[index2[0]].data1 / t2g;
3578
3579 f64 t3b[2]; // ???sigma ii / r2X
3580 t3b[0] = (1.0 - bt1data[index1[0]].data2[dir1[0]][n2] * bt1data[index1[0]].data2[dir1[0]][n2]) / bt1data[index1[0]].data1;
3581 t3b[1] = (1.0 - bt1data[index1[1]].data2[dir1[1]][n2] * bt1data[index1[1]].data2[dir1[1]][n2]) / bt1data[index1[1]].data1;
3582
3583 i32s n9[2];
3584 n9[0] = (n2 + 1) % 3; // index j
3585 n9[1] = (n2 + 2) % 3; // index k
3586
3587 f64 t3c[2]; // ???sigma ij / r2X
3588 t3c[0] = -bt1data[index1[0]].data2[dir1[0]][n2] * bt1data[index1[0]].data2[dir1[0]][n9[0]] / bt1data[index1[0]].data1;
3589 t3c[1] = -bt1data[index1[1]].data2[dir1[1]][n2] * bt1data[index1[1]].data2[dir1[1]][n9[0]] / bt1data[index1[1]].data1;
3590
3591 f64 t3d[2]; // ???sigma ik / r2X
3592 t3d[0] = -bt1data[index1[0]].data2[dir1[0]][n2] * bt1data[index1[0]].data2[dir1[0]][n9[1]] / bt1data[index1[0]].data1;
3593 t3d[1] = -bt1data[index1[1]].data2[dir1[1]][n2] * bt1data[index1[1]].data2[dir1[1]][n9[1]] / bt1data[index1[1]].data1;
3594
3595 bt3_vector[n1].dbv[1][0][n2][n2] = (t3c[0] * bt1data[index1[1]].data2[dir1[1]][n9[1]] -
3596 t3d[0] * bt1data[index1[1]].data2[dir1[1]][n9[0]] + t2i[n2] * t3a[0]) / t2h;
3597 bt3_vector[n1].dbv[1][0][n2][n9[0]] = (t3d[0] * bt1data[index1[1]].data2[dir1[1]][n2] -
3598 t3b[0] * bt1data[index1[1]].data2[dir1[1]][n9[1]] + t2i[n9[0]] * t3a[0]) / t2h;
3599 bt3_vector[n1].dbv[1][0][n2][n9[1]] = (t3b[0] * bt1data[index1[1]].data2[dir1[1]][n9[0]] -
3600 t3c[0] * bt1data[index1[1]].data2[dir1[1]][n2] + t2i[n9[1]] * t3a[0]) / t2h;
3601
3602 bt3_vector[n1].dbv[1][2][n2][n2] = (t3d[1] * bt1data[index1[0]].data2[dir1[0]][n9[0]] -
3603 t3c[1] * bt1data[index1[0]].data2[dir1[0]][n9[1]] + t2i[n2] * t3a[1]) / t2h;
3604 bt3_vector[n1].dbv[1][2][n2][n9[0]] = (t3b[1] * bt1data[index1[0]].data2[dir1[0]][n9[1]] -
3605 t3d[1] * bt1data[index1[0]].data2[dir1[0]][n2] + t2i[n9[0]] * t3a[1]) / t2h;
3606 bt3_vector[n1].dbv[1][2][n2][n9[1]] = (t3c[1] * bt1data[index1[0]].data2[dir1[0]][n2] -
3607 t3b[1] * bt1data[index1[0]].data2[dir1[0]][n9[0]] + t2i[n9[1]] * t3a[1]) / t2h;
3608 }
3609
3610 for (i32s n2 = 0;n2 < 3;n2++)
3611 {
3612 bt3_vector[n1].dbv[1][1][n2][0] = -(bt3_vector[n1].dbv[1][0][n2][0] + bt3_vector[n1].dbv[1][2][n2][0]);
3613 bt3_vector[n1].dbv[1][1][n2][1] = -(bt3_vector[n1].dbv[1][0][n2][1] + bt3_vector[n1].dbv[1][2][n2][1]);
3614 bt3_vector[n1].dbv[1][1][n2][2] = -(bt3_vector[n1].dbv[1][0][n2][2] + bt3_vector[n1].dbv[1][2][n2][2]);
3615 }
3616 }
3617
3618 // dipole vector... dipole vector... dipole vector... dipole vector...
3619 // dipole vector... dipole vector... dipole vector... dipole vector...
3620 // dipole vector... dipole vector... dipole vector... dipole vector...
3621
3622 f64 csd = cos(bt3_vector[n1].pbdd);
3623 f64 snd = sin(bt3_vector[n1].pbdd);
3624
3625 bt3_vector[n1].dv[0] = csd * bt3_vector[n1].bv[0][0] - snd * bt3_vector[n1].bv[1][0];
3626 bt3_vector[n1].dv[1] = csd * bt3_vector[n1].bv[0][1] - snd * bt3_vector[n1].bv[1][1];
3627 bt3_vector[n1].dv[2] = csd * bt3_vector[n1].bv[0][2] - snd * bt3_vector[n1].bv[1][2];
3628
3629 ////////// debug!!!
3630 //energy_bt3 += bt3_vector[n1].dv[0] + bt3_vector[n1].dv[1] + bt3_vector[n1].dv[2];
3631 ////////// debug!!!
3632
3633 if (p1 > 0)
3634 {
3635 for (i32s n2 = 0;n2 < 3;n2++)
3636 {
3637 for (i32s n3 = 0;n3 < 3;n3++)
3638 {
3639 bt3_vector[n1].ddv[0][n2][n3] = csd * bt3_vector[n1].dbv[0][0][n2][n3] - snd * bt3_vector[n1].dbv[1][0][n2][n3];
3640 bt3_vector[n1].ddv[1][n2][n3] = csd * bt3_vector[n1].dbv[0][1][n2][n3] - snd * bt3_vector[n1].dbv[1][1][n2][n3];
3641 bt3_vector[n1].ddv[2][n2][n3] = csd * bt3_vector[n1].dbv[0][2][n2][n3] - snd * bt3_vector[n1].dbv[1][2][n2][n3];
3642 bt3_vector[n1].ddv[3][n2][n3] = 0.0;
3643 }
3644 }
3645
3646 if (bt3_vector[n1].dip_ttype == TTYPE_LOOP)
3647 {
3648 for (i32s n2 = 0;n2 < 3;n2++)
3649 {
3650 for (i32s n3 = 0;n3 < 3;n3++)
3651 {
3652 bt3_vector[n1].ddv[0][n2][n3] -= snd * dpbdd[0][n2] * bt3_vector[n1].bv[0][n3] + csd * dpbdd[0][n2] * bt3_vector[n1].bv[1][n3];
3653 bt3_vector[n1].ddv[1][n2][n3] -= snd * dpbdd[1][n2] * bt3_vector[n1].bv[0][n3] + csd * dpbdd[1][n2] * bt3_vector[n1].bv[1][n3];
3654 bt3_vector[n1].ddv[2][n2][n3] -= snd * dpbdd[2][n2] * bt3_vector[n1].bv[0][n3] + csd * dpbdd[2][n2] * bt3_vector[n1].bv[1][n3];
3655 bt3_vector[n1].ddv[3][n2][n3] -= snd * dpbdd[3][n2] * bt3_vector[n1].bv[0][n3] + csd * dpbdd[3][n2] * bt3_vector[n1].bv[1][n3];
3656 }
3657 }
3658 }
3659
3660 ////////// debug!!!
3661 //for (i32s n2 = 0;n2 < 3;n2++) {
3662 //for (i32s n3 = 0;n3 < 3;n3++) {
3663 //d1[l2g_sf[atmi[0]] * 3 + n2] += bt3_vector[n1].ddv[0][n2][n3];
3664 //d1[l2g_sf[atmi[1]] * 3 + n2] += bt3_vector[n1].ddv[1][n2][n3];
3665 //d1[l2g_sf[atmi[2]] * 3 + n2] += bt3_vector[n1].ddv[2][n2][n3];
3666 //d1[l2g_sf[atmi[3]] * 3 + n2] += bt3_vector[n1].ddv[3][n2][n3]; } }
3667 ////////// debug!!!
3668
3669 }
3670 }
3671 }
3672
3673 // it seems that we compute here the same cross product than in BT3... combine ??????????
3674 // it seems that we compute here the same cross product than in BT3... combine ??????????
3675 // it seems that we compute here the same cross product than in BT3... combine ??????????
3676
ComputeBT4(i32u p1)3677 void eng1_sf::ComputeBT4(i32u p1)
3678 {
3679 energy_bt4a = 0.0;
3680 energy_bt4b = 0.0;
3681
3682 for (i32u n1 = 0;n1 < bt4_vector.size();n1++)
3683 {
3684 i32s * atma = bt2_vector[bt4_vector[n1].index2].atmi;
3685 i32s atmb = bt1_vector[bt4_vector[n1].index1].atmi[1];
3686
3687 i32s * index1a = bt2_vector[bt4_vector[n1].index2].index1;
3688 bool * dir1a = bt2_vector[bt4_vector[n1].index2].dir1;
3689
3690 i32s index1b = bt4_vector[n1].index1;
3691 i32s index2 = bt4_vector[n1].index2;
3692
3693 f64 t1a[3];
3694 f64 t1b = 1.0 - bt2data[index2].data1 * bt2data[index2].data1; f64 t1c = sqrt(t1b);
3695 t1a[0] = bt1data[index1a[0]].data2[dir1a[0]][1] * bt1data[index1a[1]].data2[dir1a[1]][2] -
3696 bt1data[index1a[0]].data2[dir1a[0]][2] * bt1data[index1a[1]].data2[dir1a[1]][1];
3697 t1a[1] = bt1data[index1a[0]].data2[dir1a[0]][2] * bt1data[index1a[1]].data2[dir1a[1]][0] -
3698 bt1data[index1a[0]].data2[dir1a[0]][0] * bt1data[index1a[1]].data2[dir1a[1]][2];
3699 t1a[2] = bt1data[index1a[0]].data2[dir1a[0]][0] * bt1data[index1a[1]].data2[dir1a[1]][1] -
3700 bt1data[index1a[0]].data2[dir1a[0]][1] * bt1data[index1a[1]].data2[dir1a[1]][0];
3701
3702 f64 t2a[3];
3703 f64 t2b = 2.0 * (bt2data[index2].data1 + 1.0); f64 t2c = sqrt(t2b);
3704 t2a[0] = bt1data[index1a[0]].data2[dir1a[0]][0] + bt1data[index1a[1]].data2[dir1a[1]][0];
3705 t2a[1] = bt1data[index1a[0]].data2[dir1a[0]][1] + bt1data[index1a[1]].data2[dir1a[1]][1];
3706 t2a[2] = bt1data[index1a[0]].data2[dir1a[0]][2] + bt1data[index1a[1]].data2[dir1a[1]][2];
3707
3708 static const f64 sb = sin(M_PI * BETA / 180.0); // sin(beta) = cos(gamma)
3709 static const f64 cb = cos(M_PI * BETA / 180.0); // cos(beta) = sin(gamma)
3710
3711 f64 t3a[3]; // vector e
3712 t3a[0] = cb * t1a[0] / t1c - sb * t2a[0] / t2c;
3713 t3a[1] = cb * t1a[1] / t1c - sb * t2a[1] / t2c;
3714 t3a[2] = cb * t1a[2] / t1c - sb * t2a[2] / t2c;
3715
3716 f64 t3b[3]; // cos(kappa)
3717 t3b[0] = t3a[0] * bt1data[index1b].data2[1][0];
3718 t3b[1] = t3a[1] * bt1data[index1b].data2[1][1];
3719 t3b[2] = t3a[2] * bt1data[index1b].data2[1][2];
3720 f64 t3c = t3b[0] + t3b[1] + t3b[2];
3721
3722 if (t3c < -1.0) t3c = -1.0; // domain check...
3723 if (t3c > +1.0) t3c = +1.0; // domain check...
3724
3725 // here like in BT2 if the kappa-angles approach 180.0 deg we will have some
3726 // serious numerical problems in later stages. here is how we can monitor this...
3727
3728 bool problem1 = (t3c < -0.999); // is it too close to 180 deg ???
3729 if (problem1) { cout << "BUG: BT4 ang -> 180.0 deg." << endl; exit(EXIT_FAILURE); }
3730
3731 // f = a(x-b)^2
3732 // df/dx = 2a(x-b)
3733
3734 f64 t3g = t3c - bt4_vector[n1].opt;
3735 f64 t9a = bt4_vector[n1].fc * t3g * t3g;
3736 f64 t9b = 2.0 * bt4_vector[n1].fc * t3g;
3737
3738 energy_bt4a += t9a;
3739
3740 // f = a/(x-1)^2
3741 // df/dx = -2a/(x-1)^3
3742
3743 // rep-fc is here (45->0) smaller than in BT2 (120->180).
3744 // rep-fc is here (45->0) smaller than in BT2 (120->180).
3745 // rep-fc is here (45->0) smaller than in BT2 (120->180).
3746 const f64 rep_weight = 0.030;
3747
3748 f64 t9x = t3c - 1.0;
3749
3750 f64 t9y = t9x * t9x;
3751 energy_bt4a += rep_weight * myprm->wrep / t9y;
3752
3753 f64 t9z = t9y * t9x;
3754 t9b -= 2.0 * rep_weight * myprm->wrep / t9z;
3755
3756 f64 t3d = 1.0 - t3c * t3c; // sin(kappa)^2
3757 f64 t3e = sqrt(t3d); // sin(kappa)
3758
3759 // here we will calculate cos(lambda) with the derivatives. however, we also have terms which
3760 // depend on sin(lambda), and there we will run into problems at 0 deg and 180 deg.
3761
3762 // this problem is mainly due to bad design, but the sine-terms are still quite important...
3763
3764 // easiest way out is to cheat a bit and use the LOWLIMIT-idea also here... this is certainly
3765 // NOT an ideal fix, since it will affect also accuracy, but still better than dividing by zero.
3766
3767 f64 t4a = 0.0; // cos(lambda)
3768 for (i32s n2 = 0;n2 < 3;n2++)
3769 {
3770 t4a += sb / t3e * bt1data[index1b].data2[1][n2] * t1a[n2] / t1c;
3771 t4a += cb / t3e * bt1data[index1b].data2[1][n2] * t2a[n2] / t2c;
3772 }
3773
3774 if (t4a < -1.0) t4a = -1.0; // domain check...
3775 if (t4a > +1.0) t4a = +1.0; // domain check...
3776
3777 f64 t4b[3]; // vector h
3778 t4b[0] = (bt1data[index1b].data2[1][0] - t3a[0] * t3c) / t3e;
3779 t4b[1] = (bt1data[index1b].data2[1][1] - t3a[1] * t3c) / t3e;
3780 t4b[2] = (bt1data[index1b].data2[1][2] - t3a[2] * t3c) / t3e;
3781
3782 f64 t4c[3]; // vector k
3783 t4c[0] = (t2a[1] * t1a[2] - t2a[2] * t1a[1]) / (t2c * t1c);
3784 t4c[1] = (t2a[2] * t1a[0] - t2a[0] * t1a[2]) / (t2c * t1c);
3785 t4c[2] = (t2a[0] * t1a[1] - t2a[1] * t1a[0]) / (t2c * t1c);
3786
3787 f64 t4d = t4b[0] * t4c[0] + t4b[1] * t4c[1] + t4b[2] * t4c[2];
3788
3789 f64 t4e = sqrt(1.0 - t4a * t4a); // sin(lambda)
3790
3791 if (t4e > +1.0) t4e = +1.0; // domain check...
3792 if (t4e < LOWLIMIT) t4e = LOWLIMIT; // apply the LOWLIMIT cheat...
3793 if (t4d < 0.0) t4e = -t4e; // sign check...
3794
3795 f64 t9c = t4a * t4a; // x*x
3796 f64 t9d = 4.0 * t9c - 1.0; // 4*x*x-1
3797
3798 // f1 = a*x + b*(2*x*x-1) + c*x*(4*x*x-3)
3799 // df1/dx = a + 4*b*x + 3*c*(4*x*x-1)
3800
3801 f64 t9e = bt4_vector[n1].fscos[0] * t4a;
3802 t9e += bt4_vector[n1].fscos[1] * (2.0 * t9c - 1.0);
3803 t9e += bt4_vector[n1].fscos[2] * t4a * (4.0 * t9c - 3.0);
3804
3805 f64 t9f = bt4_vector[n1].fscos[0];
3806 t9f += 4.0 * bt4_vector[n1].fscos[1] * t4a;
3807 t9f += 3.0 * bt4_vector[n1].fscos[2] * t9d;
3808
3809 // f2 = d*Sy + e*2*Sy*x + f*Sy*(4*x*x-1)
3810 // df2/dx = -d*x/Sy + 2e*(-x*x/Sy+Sy) + f*(-x/Sy*(4*x*x-1)+8*x*Sy)
3811
3812 t9e += bt4_vector[n1].fssin[0] * t4e;
3813 t9e += 2.0 * bt4_vector[n1].fssin[1] * t4e * t4a;
3814 t9e += bt4_vector[n1].fssin[2] * t4e * t9d;
3815
3816 t9f -= bt4_vector[n1].fssin[0] * t4a / t4e;
3817 t9f += 2.0 * bt4_vector[n1].fssin[1] * (t4e - t9c / t4e);
3818 t9f += bt4_vector[n1].fssin[2] * (8.0 * t4a * t4e - t4a * t9d / t4e);
3819
3820 energy_bt4b += t9e;
3821
3822 if (p1 > 0)
3823 {
3824 for (i32s n2 = 0;n2 < 3;n2++)
3825 {
3826 f64 t6a[2]; // epsilon i
3827 t6a[0] = bt2data[index2].data2[0][n2] * bt2data[index2].data1 / t1b;
3828 t6a[1] = bt2data[index2].data2[2][n2] * bt2data[index2].data1 / t1b;
3829
3830 f64 t6b[2]; // sigma ii / r2X
3831 t6b[0] = (1.0 - bt1data[index1a[0]].data2[dir1a[0]][n2] * bt1data[index1a[0]].data2[dir1a[0]][n2]) / bt1data[index1a[0]].data1;
3832 t6b[1] = (1.0 - bt1data[index1a[1]].data2[dir1a[1]][n2] * bt1data[index1a[1]].data2[dir1a[1]][n2]) / bt1data[index1a[1]].data1;
3833
3834 i32s n3[2];
3835 n3[0] = (n2 + 1) % 3; // index j
3836 n3[1] = (n2 + 2) % 3; // index k
3837
3838 f64 t6c[2]; // sigma ij / r2X
3839 t6c[0] = -bt1data[index1a[0]].data2[dir1a[0]][n2] * bt1data[index1a[0]].data2[dir1a[0]][n3[0]] / bt1data[index1a[0]].data1;
3840 t6c[1] = -bt1data[index1a[1]].data2[dir1a[1]][n2] * bt1data[index1a[1]].data2[dir1a[1]][n3[0]] / bt1data[index1a[1]].data1;
3841
3842 f64 t6d[2]; // sigma ik / r2X
3843 t6d[0] = -bt1data[index1a[0]].data2[dir1a[0]][n2] * bt1data[index1a[0]].data2[dir1a[0]][n3[1]] / bt1data[index1a[0]].data1;
3844 t6d[1] = -bt1data[index1a[1]].data2[dir1a[1]][n2] * bt1data[index1a[1]].data2[dir1a[1]][n3[1]] / bt1data[index1a[1]].data1;
3845
3846 f64 t5a[2][3]; // components of dc/di
3847 t5a[0][n2] = (t6c[0] * bt1data[index1a[1]].data2[dir1a[1]][n3[1]] -
3848 t6d[0] * bt1data[index1a[1]].data2[dir1a[1]][n3[0]] + t1a[n2] * t6a[0]) / t1c;
3849 t5a[0][n3[0]] = (t6d[0] * bt1data[index1a[1]].data2[dir1a[1]][n2] -
3850 t6b[0] * bt1data[index1a[1]].data2[dir1a[1]][n3[1]] + t1a[n3[0]] * t6a[0]) / t1c;
3851 t5a[0][n3[1]] = (t6b[0] * bt1data[index1a[1]].data2[dir1a[1]][n3[0]] -
3852 t6c[0] * bt1data[index1a[1]].data2[dir1a[1]][n2] + t1a[n3[1]] * t6a[0]) / t1c;
3853 t5a[1][n2] = (t6d[1] * bt1data[index1a[0]].data2[dir1a[0]][n3[0]] -
3854 t6c[1] * bt1data[index1a[0]].data2[dir1a[0]][n3[1]] + t1a[n2] * t6a[1]) / t1c;
3855 t5a[1][n3[0]] = (t6b[1] * bt1data[index1a[0]].data2[dir1a[0]][n3[1]] -
3856 t6d[1] * bt1data[index1a[0]].data2[dir1a[0]][n2] + t1a[n3[0]] * t6a[1]) / t1c;
3857 t5a[1][n3[1]] = (t6c[1] * bt1data[index1a[0]].data2[dir1a[0]][n2] -
3858 t6b[1] * bt1data[index1a[0]].data2[dir1a[0]][n3[0]] + t1a[n3[1]] * t6a[1]) / t1c;
3859
3860 f64 t5b[2][3]; // components of dd/di
3861 f64 t6e = (bt1data[index1a[0]].data2[dir1a[0]][n2] + bt1data[index1a[1]].data2[dir1a[1]][n2]) / t2b;
3862 t5b[0][n2] = (t6b[0] - bt2data[index2].data2[0][n2] * t6e) / t2c;
3863 t5b[1][n2] = (t6b[1] - bt2data[index2].data2[2][n2] * t6e) / t2c;
3864 t6e = (bt1data[index1a[0]].data2[dir1a[0]][n3[0]] + bt1data[index1a[1]].data2[dir1a[1]][n3[0]]) / t2b;
3865 t5b[0][n3[0]] = (t6c[0] - bt2data[index2].data2[0][n2] * t6e) / t2c;
3866 t5b[1][n3[0]] = (t6c[1] - bt2data[index2].data2[2][n2] * t6e) / t2c;
3867 t6e = (bt1data[index1a[0]].data2[dir1a[0]][n3[1]] + bt1data[index1a[1]].data2[dir1a[1]][n3[1]]) / t2b;
3868 t5b[0][n3[1]] = (t6d[0] - bt2data[index2].data2[0][n2] * t6e) / t2c;
3869 t5b[1][n3[1]] = (t6d[1] - bt2data[index2].data2[2][n2] * t6e) / t2c;
3870
3871 f64 t5c[4] = { 0.0, 0.0, 0.0, 0.0 }; // dcos(kappa)/di
3872 for (i32s n3 = 0;n3 < 3;n3++)
3873 {
3874 t5c[0] += bt1data[index1b].data2[1][n3] * (cb * t5a[0][n3] - sb * t5b[0][n3]);
3875 t5c[2] += bt1data[index1b].data2[1][n3] * (cb * t5a[1][n3] - sb * t5b[1][n3]);
3876
3877 if (n2 != n3) t6e = -bt1data[index1b].data2[1][n2] * bt1data[index1b].data2[1][n3];
3878 else t6e = 1.0 - bt1data[index1b].data2[1][n2] * bt1data[index1b].data2[1][n2];
3879 t5c[3] += t3a[n3] * t6e / bt1data[index1b].data1;
3880 }
3881
3882 t5c[1] = -(t5c[0] + t5c[2] + t5c[3]);
3883
3884 d1[l2g_sf[atma[0]] * 3 + n2] += t9b * t5c[0];
3885 d1[l2g_sf[atma[1]] * 3 + n2] += t9b * t5c[1];
3886 d1[l2g_sf[atma[2]] * 3 + n2] += t9b * t5c[2];
3887 d1[l2g_sf[atmb] * 3 + n2] += t9b * t5c[3];
3888
3889 f64 t7e[3];
3890 t7e[0] = t3c * t5c[0] / t3d;
3891 t7e[1] = t3c * t5c[2] / t3d;
3892 t7e[2] = t3c * t5c[3] / t3d;
3893
3894 f64 t5d[4] = { 0.0, 0.0, 0.0, 0.0 }; // dcos(lambda)/di
3895 for (i32s n3 = 0;n3 < 3;n3++)
3896 {
3897 f64 t7a = t1a[n3] / t1c; // f64 t7b = t7a * t7a;
3898 f64 t7c = t2a[n3] / t2c; // f64 t7d = t7c * t7c;
3899
3900 t5d[0] += sb * bt1data[index1b].data2[1][n3] * (t7e[0] * t7a + t5a[0][n3]) / t3e;
3901 t5d[0] += cb * bt1data[index1b].data2[1][n3] * (t7e[0] * t7c + t5b[0][n3]) / t3e;
3902
3903 t5d[2] += sb * bt1data[index1b].data2[1][n3] * (t7e[1] * t7a + t5a[1][n3]) / t3e;
3904 t5d[2] += cb * bt1data[index1b].data2[1][n3] * (t7e[1] * t7c + t5b[1][n3]) / t3e;
3905
3906 if (n2 != n3) t6e = -bt1data[index1b].data2[1][n2] * bt1data[index1b].data2[1][n3];
3907 else t6e = 1.0 - bt1data[index1b].data2[1][n2] * bt1data[index1b].data2[1][n2];
3908
3909 f64 t7f = t7e[2] * bt1data[index1b].data2[1][n3] + t6e / bt1data[index1b].data1;
3910 t5d[3] += (sb * t7a + cb * t7c) * t7f / t3e;
3911 }
3912
3913 t5d[1] = -(t5d[0] + t5d[2] + t5d[3]);
3914
3915 d1[l2g_sf[atma[0]] * 3 + n2] += t9f * t5d[0];
3916 d1[l2g_sf[atma[1]] * 3 + n2] += t9f * t5d[1];
3917 d1[l2g_sf[atma[2]] * 3 + n2] += t9f * t5d[2];
3918 d1[l2g_sf[atmb] * 3 + n2] += t9f * t5d[3];
3919 }
3920 }
3921 }
3922 }
3923
ComputeNBT1(i32u p1)3924 void eng1_sf::ComputeNBT1(i32u p1)
3925 {
3926 energy_nbt1a = 0.0;
3927 energy_nbt1b = 0.0;
3928 energy_nbt1c = 0.0; // not needed anymore???
3929
3930 atom ** atmtab = GetSetup()->GetSFAtoms();
3931
3932 // an additional pass for the boundary potential (optional).
3933 // an additional pass for the boundary potential (optional).
3934 // an additional pass for the boundary potential (optional).
3935
3936 if (use_bp)
3937 {
3938 if (nd_eval != NULL) nd_eval->AddCycle();
3939 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount();n1++)
3940 {
3941 f64 radius = bp_rad_solute;
3942 f64 fc = bp_fc_solute;
3943
3944 if (atmtab[n1]->flags & ATOMFLAG_IS_SOLVENT_ATOM)
3945 {
3946 radius = bp_rad_solvent;
3947 fc = bp_fc_solvent;
3948 }
3949
3950 f64 t1a[3]; f64 t1b = 0.0;
3951 for (i32s n2 = 0;n2 < 3;n2++)
3952 {
3953 f64 t9a = 0.0;//bp_center[n2];
3954 f64 t9b = crd[l2g_sf[n1] * 3 + n2];
3955
3956 t1a[n2] = t9a - t9b;
3957 t1b += t1a[n2] * t1a[n2];
3958 }
3959
3960 f64 t1c = sqrt(t1b);
3961
3962 if (nd_eval != NULL && (atmtab[n1]->flags & ATOMFLAG_MEASURE_ND_RDF)) nd_eval->AddValue(t1c);
3963
3964 if (rdf_eval != NULL && rdf_eval->count_begin > -0.5)
3965 {
3966 bool is_in_rdf_counting_window = true;
3967 if (t1c < rdf_eval->count_begin) is_in_rdf_counting_window = false;
3968 if (t1c >= rdf_eval->count_end) is_in_rdf_counting_window = false;
3969
3970 if (is_in_rdf_counting_window) atmtab[n1]->flags |= ATOMFLAG_COUNT_IN_RDF;
3971 else atmtab[n1]->flags &= (~ATOMFLAG_COUNT_IN_RDF);
3972 }
3973
3974 if (t1c < radius) continue;
3975
3976 // f = a(x-b)^2
3977 // df/dx = 2a(x-b)
3978
3979 f64 t2a = t1c - radius;
3980 f64 t2b = fc * t2a * t2a;
3981
3982 energy_bt1 += t2b;
3983
3984 if (p1 > 0)
3985 {
3986 f64 t2c = 2.0 * fc * t2a;
3987
3988 for (i32s n2 = 0;n2 < 3;n2++)
3989 {
3990 f64 t2d = (t1a[n2] / t1c) * t2c;
3991
3992 d1[l2g_sf[n1] * 3 + n2] -= t2d;
3993 }
3994 }
3995 }
3996 }
3997
3998 i32s limit = GetSetup()->GetSFAtomCount() - num_solvent;
3999 for (i32u n1 = 0;n1 < nbt1_vector.size();n1++)
4000 {
4001 i32s * atmi = nbt1_vector[n1].atmi;
4002
4003 f64 t1a[3]; f64 t1b = 0.0;
4004 for (i32s n2 = 0;n2 < 3;n2++)
4005 {
4006 f64 t2a = crd[l2g_sf[atmi[0]] * 3 + n2];
4007 f64 t2b = crd[l2g_sf[atmi[1]] * 3 + n2];
4008
4009 t1a[n2] = t2a - t2b;
4010 t1b += t1a[n2] * t1a[n2];
4011 }
4012
4013 f64 t1c = sqrt(t1b);
4014
4015 // f = ((x+s)/a)^-9 - ((x+s)/b)^-6
4016 // df/dx = -9/a((x+s)/a)^-10 + 6/b((x+s)/b)^-7
4017
4018 const f64 buff = 0.005;
4019 f64 tmp1 = t1c + buff;
4020
4021 f64 t3a = tmp1 / nbt1_vector[n1].data[0];
4022 f64 t3b = tmp1 / nbt1_vector[n1].data[1];
4023
4024 f64 t4a = t3a * t3a; f64 t4b = t4a * t3a; f64 t4c = t4a * t4b; // ^2 ^3 ^5
4025 f64 t5a = t3b * t3b; f64 t5b = t5a * t3b; //f64 t5c = t5a * t5b; // ^2 ^3 ^5
4026
4027 // energy_nbt1a += 1.0 / (t4b * t4b * t4b) - 1.0 / (t5b * t5b); // 6-9 << also see InitLenJon()!!!
4028 energy_nbt1a += 1.0 / (t4c * t4c * t4a) - 1.0 / (t5b * t5b); // 6-12 << also see InitLenJon()!!!
4029 // energy_nbt1a += 1.0 / (t4c * t4c * t4a) - 1.0 / (t5b * t5b * t5b); // 9-12 << also see InitLenJon()!!!
4030
4031 f64 t6c = 0.0;
4032 if (atmi[0] < limit && atmi[1] < limit)
4033 {
4034 // e = 2 + 76 * (( r / A ) ^ n) / (1 + ( r / A ) ^ n) = 2 + 76 * f / g
4035 // de/dr = 76 * (( g * Df - f * Dg ) / g ^ 2 )
4036
4037 f64 e1 = solv_exp[0][atmi[0]];
4038 f64 e2 = solv_exp[0][atmi[1]];
4039
4040 f64 eps_n = myprm->epsilon1 + t1b * myprm->epsilon2; // assume constant!!!
4041 f64 eps_A = myprm->epsilon3 - log(1.0 + (e1 + e2) * myprm->epsilon4 + e1 * e2 * myprm->epsilon5); // assume constant!!!
4042 if (eps_A < 0.2) { cout << "BUGGER!!! " << eps_A << endl; exit(EXIT_FAILURE); }
4043
4044 f64 t7a = t1c / eps_A;
4045 f64 t7b = pow(t7a, eps_n);
4046 f64 t7c = 1.0 + t7b;
4047 f64 t7d = 2.0 + 76.0 * (t7b / t7c);
4048
4049 f64 t7x = eps_n * pow(t7a, eps_n - 1.0) / eps_A;
4050 f64 t7y = 76.0 * ((t7c * t7x - t7b * t7x) / (t7c * t7c));
4051
4052 f64 qq = 138.9354518 * charge2[atmi[0]] * charge2[atmi[1]] * myprm->wescc;
4053
4054 // f = Q/(e*r)
4055 // df/dr = -Q * ((1/e*r^2) + (de/dr)/(e^2*r))
4056
4057 energy_nbt1b += qq / (t7d * t1c);
4058
4059 t6c = -qq * (1.0 / (t7d * t1b) + t7y / (t7d * t7d * t1c));
4060 }
4061
4062 if (p1 > 0)
4063 {
4064 // f64 t6a = 9.0 / (nbt1_vector[n1].data[0] * t4c * t4c); // 6-9
4065 // f64 t6b = 6.0 / (nbt1_vector[n1].data[1] * t5a * t5a * t5b);
4066 f64 t6a = 12.0 / (nbt1_vector[n1].data[0] * t4c * t4c * t4b); // 6-12
4067 f64 t6b = 6.0 / (nbt1_vector[n1].data[1] * t5a * t5a * t5b);
4068 // f64 t6a = 12.0 / (nbt1_vector[n1].data[0] * t4c * t4c * t4b); // 9-12
4069 // f64 t6b = 9.0 / (nbt1_vector[n1].data[1] * t5c * t5c);
4070
4071 f64 t2a = (t6b - t6a) + t6c;
4072
4073 for (i32s n2 = 0;n2 < 3;n2++)
4074 {
4075 f64 t1e = (t1a[n2] / t1c) * t2a;
4076
4077 d1[l2g_sf[atmi[0]] * 3 + n2] += t1e;
4078 d1[l2g_sf[atmi[1]] * 3 + n2] -= t1e;
4079 }
4080 }
4081 }
4082 }
4083
ComputeNBT2(i32u p1)4084 void eng1_sf::ComputeNBT2(i32u p1)
4085 {
4086 energy_nbt2a = 0.0;
4087 energy_nbt2b = 0.0; // not needed anymore???
4088 energy_nbt2c = 0.0;
4089
4090 // dipole moment for the peptide unit:
4091 // CRC handbook of Chem & Phys, 1st Student Edition, 1988, CRC Press Inc. (page E-53)
4092 // "acetyl methylamine" 3.73 debyes -> 3.73 * 3.338e-30 / 1.6021773e-19 * 1.0e+09 = 0.07771137
4093
4094 /*##############################################*/
4095 /*##############################################*/
4096
4097 for (i32u n1 = 0;n1 < bt3_vector.size();n1++)
4098 {
4099 if (bt3_vector[n1].skip) continue; // skip all X-pro cases...
4100 i32s * atmi = bt3_vector[n1].atmi;
4101
4102 for (i32s n2 = 0;n2 < GetSetup()->GetSFAtomCount() - num_solvent;n2++)
4103 {
4104 if ((charge2[n2] == 0.0)) continue;
4105 if (atmi[1] == (i32s) n2 || atmi[2] == (i32s) n2) continue;
4106
4107 f64 t1a[3];
4108 t1a[0] = crd[l2g_sf[n2] * 3 + 0] - 0.5 * (crd[l2g_sf[atmi[1]] * 3 + 0] + crd[l2g_sf[atmi[2]] * 3 + 0]);
4109 t1a[1] = crd[l2g_sf[n2] * 3 + 1] - 0.5 * (crd[l2g_sf[atmi[1]] * 3 + 1] + crd[l2g_sf[atmi[2]] * 3 + 1]);
4110 t1a[2] = crd[l2g_sf[n2] * 3 + 2] - 0.5 * (crd[l2g_sf[atmi[1]] * 3 + 2] + crd[l2g_sf[atmi[2]] * 3 + 2]);
4111
4112 f64 t1b = t1a[0] * t1a[0] + t1a[1] * t1a[1] + t1a[2] * t1a[2];
4113 f64 t1c = sqrt(t1b); f64 t1d[2][3]; f64 t1e[2][3][3];
4114
4115 // t1d = dr/dx, and therefore t1d[0] is a unit vector...
4116 // t1d = dr/dx, and therefore t1d[0] is a unit vector...
4117 // t1d = dr/dx, and therefore t1d[0] is a unit vector...
4118
4119 for (i32s n3 = 0;n3 < 3;n3++)
4120 {
4121 t1d[0][n3] = t1a[n3] / t1c;
4122 t1d[1][n3] = -0.5 * t1d[0][n3];
4123 }
4124
4125 f64 t1f[3]; f64 t1g[5][3]; // theta
4126 t1f[0] = t1d[0][0] * bt3_vector[n1].dv[0];
4127 t1f[1] = t1d[0][1] * bt3_vector[n1].dv[1];
4128 t1f[2] = t1d[0][2] * bt3_vector[n1].dv[2];
4129 f64 t1h = t1f[0] + t1f[1] + t1f[2];
4130
4131 ////////// debug!!!
4132 //energy_nbt2 += t1h;
4133 ////////// debug!!!
4134
4135 if (p1 > 0)
4136 {
4137 for (i32s n3 = 0;n3 < 3;n3++)
4138 {
4139 for (i32s n4 = 0;n4 < 3;n4++)
4140 {
4141 if (n3 != n4) t1e[0][n3][n4] = -t1d[0][n3] * t1d[0][n4] / t1c;
4142 else t1e[0][n3][n4] = (1.0 - t1d[0][n3] * t1d[0][n3]) / t1c;
4143 t1e[1][n3][n4] = -0.5 * t1e[0][n3][n4];
4144 }
4145 }
4146
4147 for (i32s n3 = 0;n3 < 3;n3++)
4148 {
4149 t1g[0][n3] = t1g[1][n3] = t1g[2][n3] = t1g[3][n3] = t1g[4][n3] = 0.0;
4150
4151 for (i32s n4 = 0;n4 < 3;n4++)
4152 {
4153 t1g[0][n3] += t1d[0][n4] * bt3_vector[n1].ddv[0][n3][n4];
4154 t1g[1][n3] += t1d[0][n4] * bt3_vector[n1].ddv[1][n3][n4] + t1e[1][n3][n4] * bt3_vector[n1].dv[n4];
4155 t1g[2][n3] += t1d[0][n4] * bt3_vector[n1].ddv[2][n3][n4] + t1e[1][n3][n4] * bt3_vector[n1].dv[n4];
4156 t1g[3][n3] += t1d[0][n4] * bt3_vector[n1].ddv[3][n3][n4];
4157 t1g[4][n3] += t1e[0][n3][n4] * bt3_vector[n1].dv[n4];
4158 }
4159
4160 ////////// debug!!!
4161 //d1[l2g_sf[atmi[0]] * 3 + n3] += t1g[0][n3];
4162 //d1[l2g_sf[atmi[1]] * 3 + n3] += t1g[1][n3];
4163 //d1[l2g_sf[atmi[2]] * 3 + n3] += t1g[2][n3];
4164 //d1[l2g_sf[atmi[3]] * 3 + n3] += t1g[3][n3];
4165 //d1[l2g_sf[n2] * 3 + n3] += t1g[4][n3];
4166 ////////// debug!!!
4167
4168 }
4169 }
4170
4171 f64 t9x = t1b * t1b; // this is r^4 !!!!!!!
4172 f64 t9y = t9x * t9x; // this is r^8 !!!!!!!
4173
4174 // e = 2 + 76 * (( r / A ) ^ n) / (1 + ( r / A ) ^ n) + Z/r^9 = 2 + 76 * f / g + Z/r^9
4175 // de/dr = 76 * (( g * Df - f * Dg ) / g ^ 2 ) - 9Z/r^10
4176
4177 f64 eps_n = myprm->epsilon1 + t1b * myprm->epsilon2; // assume constant!!!
4178 f64 eps_A = 0.75; // most chg-dip interactions at surface -> use a small value!
4179
4180 f64 t7a = t1c / eps_A;
4181 f64 t7b = pow(t7a, eps_n);
4182 f64 t7c = 1.0 + t7b;
4183 f64 t7d = 2.0 + 76.0 * (t7b / t7c) + myprm->epsilon9 / (t9y * t1c);
4184
4185 f64 t7x = eps_n * pow(t7a, eps_n - 1.0) / eps_A;
4186 f64 t7y = 76.0 * ((t7c * t7x - t7b * t7x) / (t7c * t7c)) - 9.0 * myprm->epsilon9 / (t9y * t1b);
4187
4188 f64 qd = 138.9354518 * charge2[n2] * 0.077711 * myprm->dipole1 * myprm->wescd;
4189
4190 // f = Q/(e*r^2)
4191 // df/dr = -Q * ((2/e*r^3) + (de/dr)/(e^2*r^2))
4192
4193 energy_nbt2a += qd * t1h / (t7d * t1b);
4194
4195 if (p1 > 0)
4196 {
4197 f64 t3a = -qd * t1h * (2.0 / (t7d * t1b * t1c) + t7y / (t7d * t7d * t1b));
4198 f64 t3b = qd / (t7d * t1b);
4199
4200 for (i32s n3 = 0;n3 < 3;n3++)
4201 {
4202 d1[l2g_sf[atmi[1]] * 3 + n3] += t3a * t1d[1][n3];
4203 d1[l2g_sf[atmi[2]] * 3 + n3] += t3a * t1d[1][n3];
4204 d1[l2g_sf[n2] * 3 + n3] += t3a * t1d[0][n3];
4205
4206 d1[l2g_sf[atmi[0]] * 3 + n3] += t3b * t1g[0][n3];
4207 d1[l2g_sf[atmi[1]] * 3 + n3] += t3b * t1g[1][n3];
4208 d1[l2g_sf[atmi[2]] * 3 + n3] += t3b * t1g[2][n3];
4209 d1[l2g_sf[atmi[3]] * 3 + n3] += t3b * t1g[3][n3];
4210 d1[l2g_sf[n2] * 3 + n3] += t3b * t1g[4][n3];
4211 }
4212 }
4213 }
4214 }
4215
4216 /*##############################################*/
4217 /*##############################################*/
4218
4219 for (i32s n1 = 0;n1 < ((i32s) bt3_vector.size()) - 1;n1++)
4220 {
4221 if (bt3_vector[n1].skip) continue; // skip all X-pro cases...
4222 i32s * atmi1 = bt3_vector[n1].atmi;
4223
4224 for (i32s n2 = n1 + 1;n2 < (i32s) bt3_vector.size();n2++)
4225 {
4226 if (bt3_vector[n2].skip) continue; // skip all X-pro cases...
4227 i32s * atmi2 = bt3_vector[n2].atmi;
4228
4229 f64 t1a[3];
4230 t1a[0] = 0.5 * (crd[l2g_sf[atmi2[1]] * 3 + 0] + crd[l2g_sf[atmi2[2]] * 3 + 0] - crd[l2g_sf[atmi1[1]] * 3 + 0] - crd[l2g_sf[atmi1[2]] * 3 + 0]);
4231 t1a[1] = 0.5 * (crd[l2g_sf[atmi2[1]] * 3 + 1] + crd[l2g_sf[atmi2[2]] * 3 + 1] - crd[l2g_sf[atmi1[1]] * 3 + 1] - crd[l2g_sf[atmi1[2]] * 3 + 1]);
4232 t1a[2] = 0.5 * (crd[l2g_sf[atmi2[1]] * 3 + 2] + crd[l2g_sf[atmi2[2]] * 3 + 2] - crd[l2g_sf[atmi1[1]] * 3 + 2] - crd[l2g_sf[atmi1[2]] * 3 + 2]);
4233
4234 f64 t1b = t1a[0] * t1a[0] + t1a[1] * t1a[1] + t1a[2] * t1a[2];
4235 f64 t1c = sqrt(t1b); f64 t1d[3]; f64 t1e[3][3];
4236
4237 // t1d = dr/dx, and therefore is a vector of constant length 1/2...
4238 // t1d = dr/dx, and therefore is a vector of constant length 1/2...
4239 // t1d = dr/dx, and therefore is a vector of constant length 1/2...
4240
4241 t1d[0] = 0.5 * t1a[0] / t1c;
4242 t1d[1] = 0.5 * t1a[1] / t1c;
4243 t1d[2] = 0.5 * t1a[2] / t1c;
4244
4245 f64 t1f[3]; f64 t1g[5][3]; // theta1 * 1/2
4246 t1f[0] = t1d[0] * bt3_vector[n1].dv[0];
4247 t1f[1] = t1d[1] * bt3_vector[n1].dv[1];
4248 t1f[2] = t1d[2] * bt3_vector[n1].dv[2];
4249 f64 t1h = t1f[0] + t1f[1] + t1f[2];
4250
4251 f64 t1i[3]; f64 t1j[5][3]; // theta2 * 1/2
4252 t1i[0] = t1d[0] * bt3_vector[n2].dv[0];
4253 t1i[1] = t1d[1] * bt3_vector[n2].dv[1];
4254 t1i[2] = t1d[2] * bt3_vector[n2].dv[2];
4255 f64 t1k = t1i[0] + t1i[1] + t1i[2];
4256
4257 f64 t1l[3]; f64 t1m[8][3];
4258 t1l[0] = bt3_vector[n1].dv[0] * bt3_vector[n2].dv[0];
4259 t1l[1] = bt3_vector[n1].dv[1] * bt3_vector[n2].dv[1];
4260 t1l[2] = bt3_vector[n1].dv[2] * bt3_vector[n2].dv[2];
4261 f64 t1n = t1l[0] + t1l[1] + t1l[2];
4262
4263 ////////// debug!!!
4264 //energy_nbt2 += t1h;
4265 //energy_nbt2 += t1k;
4266 //energy_nbt2 += t1n;
4267 ////////// debug!!!
4268
4269 if (p1 > 0)
4270 {
4271 for (i32s n3 = 0;n3 < 3;n3++)
4272 {
4273 for (i32s n4 = 0;n4 < 3;n4++)
4274 {
4275 if (n3 != n4) t1e[n3][n4] = -t1d[n3] * t1d[n4] / t1c;
4276 else t1e[n3][n4] = (0.25 - t1d[n3] * t1d[n3]) / t1c;
4277 }
4278 }
4279
4280 for (i32s n3 = 0;n3 < 3;n3++)
4281 {
4282 t1g[0][n3] = t1g[1][n3] = t1g[2][n3] = t1g[3][n3] = t1g[4][n3] = 0.0;
4283 t1j[0][n3] = t1j[1][n3] = t1j[2][n3] = t1j[3][n3] = t1j[4][n3] = 0.0;
4284
4285 t1m[0][n3] = t1m[1][n3] = t1m[2][n3] = t1m[3][n3] = 0.0;
4286 t1m[4][n3] = t1m[5][n3] = t1m[6][n3] = t1m[7][n3] = 0.0;
4287
4288 for (i32s n4 = 0;n4 < 3;n4++)
4289 {
4290 t1g[0][n3] += t1d[n4] * bt3_vector[n1].ddv[0][n3][n4];
4291 t1g[1][n3] += t1d[n4] * bt3_vector[n1].ddv[1][n3][n4] - t1e[n3][n4] * bt3_vector[n1].dv[n4];
4292 t1g[2][n3] += t1d[n4] * bt3_vector[n1].ddv[2][n3][n4] - t1e[n3][n4] * bt3_vector[n1].dv[n4];
4293 t1g[3][n3] += t1d[n4] * bt3_vector[n1].ddv[3][n3][n4];
4294 t1g[4][n3] += t1e[n3][n4] * bt3_vector[n1].dv[n4];
4295
4296 t1j[0][n3] += t1d[n4] * bt3_vector[n2].ddv[0][n3][n4];
4297 t1j[1][n3] += t1d[n4] * bt3_vector[n2].ddv[1][n3][n4] + t1e[n3][n4] * bt3_vector[n2].dv[n4];
4298 t1j[2][n3] += t1d[n4] * bt3_vector[n2].ddv[2][n3][n4] + t1e[n3][n4] * bt3_vector[n2].dv[n4];
4299 t1j[3][n3] += t1d[n4] * bt3_vector[n2].ddv[3][n3][n4];
4300 t1j[4][n3] -= t1e[n3][n4] * bt3_vector[n2].dv[n4];
4301
4302 t1m[0][n3] += bt3_vector[n1].ddv[0][n3][n4] * bt3_vector[n2].dv[n4];
4303 t1m[1][n3] += bt3_vector[n1].ddv[1][n3][n4] * bt3_vector[n2].dv[n4];
4304 t1m[2][n3] += bt3_vector[n1].ddv[2][n3][n4] * bt3_vector[n2].dv[n4];
4305 t1m[3][n3] += bt3_vector[n1].ddv[3][n3][n4] * bt3_vector[n2].dv[n4];
4306 t1m[4][n3] += bt3_vector[n1].dv[n4] * bt3_vector[n2].ddv[0][n3][n4];
4307 t1m[5][n3] += bt3_vector[n1].dv[n4] * bt3_vector[n2].ddv[1][n3][n4];
4308 t1m[6][n3] += bt3_vector[n1].dv[n4] * bt3_vector[n2].ddv[2][n3][n4];
4309 t1m[7][n3] += bt3_vector[n1].dv[n4] * bt3_vector[n2].ddv[3][n3][n4];
4310 }
4311
4312 ////////// debug!!!
4313 //d1[l2g_sf[atmi1[0]] * 3 + n3] += t1g[0][n3];
4314 //d1[l2g_sf[atmi1[1]] * 3 + n3] += t1g[1][n3];
4315 //d1[l2g_sf[atmi1[2]] * 3 + n3] += t1g[2][n3];
4316 //d1[l2g_sf[atmi1[3]] * 3 + n3] += t1g[3][n3];
4317 //d1[l2g_sf[atmi2[1]] * 3 + n3] += t1g[4][n3];
4318 //d1[l2g_sf[atmi2[2]] * 3 + n3] += t1g[4][n3];
4319 //
4320 //d1[l2g_sf[atmi2[0]] * 3 + n3] += t1j[0][n3];
4321 //d1[l2g_sf[atmi2[1]] * 3 + n3] += t1j[1][n3];
4322 //d1[l2g_sf[atmi2[2]] * 3 + n3] += t1j[2][n3];
4323 //d1[l2g_sf[atmi2[3]] * 3 + n3] += t1j[3][n3];
4324 //d1[l2g_sf[atmi1[1]] * 3 + n3] += t1j[4][n3];
4325 //d1[l2g_sf[atmi1[2]] * 3 + n3] += t1j[4][n3];
4326 //
4327 //d1[l2g_sf[atmi1[0]] * 3 + n3] += t1m[0][n3];
4328 //d1[l2g_sf[atmi1[1]] * 3 + n3] += t1m[1][n3];
4329 //d1[l2g_sf[atmi1[2]] * 3 + n3] += t1m[2][n3];
4330 //d1[l2g_sf[atmi1[3]] * 3 + n3] += t1m[3][n3];
4331 //d1[l2g_sf[atmi2[0]] * 3 + n3] += t1m[4][n3];
4332 //d1[l2g_sf[atmi2[1]] * 3 + n3] += t1m[5][n3];
4333 //d1[l2g_sf[atmi2[2]] * 3 + n3] += t1m[6][n3];
4334 //d1[l2g_sf[atmi2[3]] * 3 + n3] += t1m[7][n3];
4335 ////////// debug!!!
4336
4337 }
4338 }
4339
4340 f64 t9x = t1b * t1b; // this is r^4 !!!!!!!
4341 f64 t9y = t9x * t9x; // this is r^8 !!!!!!!
4342
4343 // e = 2 + 76 * (( r / A ) ^ n) / (1 + ( r / A ) ^ n) + Z/r^9 = 2 + 76 * f / g + Z/r^9
4344 // de/dr = 76 * (( g * Df - f * Dg ) / g ^ 2 ) - 9Z/r^10
4345
4346 f64 eps_n = myprm->epsilon1 + t1b * myprm->epsilon2; // assume constant!!!
4347 f64 eps_A = 1.75; // most dip-dip interactions at interior -> use a large value!
4348
4349 f64 t7a = t1c / eps_A;
4350 f64 t7b = pow(t7a, eps_n);
4351 f64 t7c = 1.0 + t7b;
4352 f64 t7d = 2.0 + 76.0 * (t7b / t7c) + myprm->epsilon9 / (t9y * t1c);
4353
4354 f64 t7x = eps_n * pow(t7a, eps_n - 1.0) / eps_A;
4355 f64 t7y = 76.0 * ((t7c * t7x - t7b * t7x) / (t7c * t7c)) - 9.0 * myprm->epsilon9 / (t9y * t1b);
4356
4357 // f = Q/(e*r^3)
4358 // df/dr = -Q * ((3/e*r^4) + (de/dr)/(e^2*r^3))
4359
4360 f64 dd = 138.9354518 * 0.077711 * 0.077711 * myprm->dipole1 * myprm->dipole1;
4361
4362 // enhance alpha-helices...
4363 // enhance alpha-helices...
4364 // enhance alpha-helices...
4365
4366 if ((n2 - n1) == 3)
4367 {
4368 bool at1 = (bt3_vector[n1 + 0].dip_ttype == TTYPE_HELIX);
4369 bool at2 = (bt3_vector[n1 + 1].dip_ttype == TTYPE_HELIX);
4370 bool at3 = (bt3_vector[n1 + 2].dip_ttype == TTYPE_HELIX);
4371 bool at4 = (bt3_vector[n1 + 3].dip_ttype == TTYPE_HELIX);
4372
4373 bool same_chn = (bt3_vector[n1].atmi[3] == bt3_vector[n2].atmi[0]);
4374
4375 if (at1 && at2 && at3 && at4 && same_chn) dd *= myprm->dipole2;
4376 }
4377
4378 // enhance beta-sheets...
4379 // enhance beta-sheets...
4380 // enhance beta-sheets...
4381
4382 bool bt1 = (bt3_vector[n1].dip_ttype == TTYPE_STRAND);
4383 bool bt2 = (bt3_vector[n2].dip_ttype == TTYPE_STRAND);
4384 if (bt1 && bt2) dd *= myprm->dipole2;
4385
4386 f64 t1x = dd / (t7d * t1b * t1c); // radial part...
4387 f64 t1y = t1n - 12.0 * t1h * t1k; // angular part...
4388 energy_nbt2c += t1x * t1y;
4389
4390 if (p1 > 0)
4391 {
4392 f64 t3a = -dd * t1y * (3.0 / (t7d * t9x) + t7y / (t7d * t7d * t1b * t1c));
4393 f64 t3b = -12.0 * t1k * t1x;
4394 f64 t3c = -12.0 * t1h * t1x;
4395
4396 for (i32s n3 = 0;n3 < 3;n3++)
4397 {
4398 d1[l2g_sf[atmi1[1]] * 3 + n3] -= t3a * t1d[n3];
4399 d1[l2g_sf[atmi1[2]] * 3 + n3] -= t3a * t1d[n3];
4400 d1[l2g_sf[atmi2[1]] * 3 + n3] += t3a * t1d[n3];
4401 d1[l2g_sf[atmi2[2]] * 3 + n3] += t3a * t1d[n3];
4402
4403 d1[l2g_sf[atmi1[0]] * 3 + n3] += t3b * t1g[0][n3];
4404 d1[l2g_sf[atmi1[1]] * 3 + n3] += t3b * t1g[1][n3];
4405 d1[l2g_sf[atmi1[2]] * 3 + n3] += t3b * t1g[2][n3];
4406 d1[l2g_sf[atmi1[3]] * 3 + n3] += t3b * t1g[3][n3];
4407 d1[l2g_sf[atmi2[1]] * 3 + n3] += t3b * t1g[4][n3];
4408 d1[l2g_sf[atmi2[2]] * 3 + n3] += t3b * t1g[4][n3];
4409
4410 d1[l2g_sf[atmi2[0]] * 3 + n3] += t3c * t1j[0][n3];
4411 d1[l2g_sf[atmi2[1]] * 3 + n3] += t3c * t1j[1][n3];
4412 d1[l2g_sf[atmi2[2]] * 3 + n3] += t3c * t1j[2][n3];
4413 d1[l2g_sf[atmi2[3]] * 3 + n3] += t3c * t1j[3][n3];
4414 d1[l2g_sf[atmi1[1]] * 3 + n3] += t3c * t1j[4][n3];
4415 d1[l2g_sf[atmi1[2]] * 3 + n3] += t3c * t1j[4][n3];
4416
4417 d1[l2g_sf[atmi1[0]] * 3 + n3] += t1x * t1m[0][n3];
4418 d1[l2g_sf[atmi1[1]] * 3 + n3] += t1x * t1m[1][n3];
4419 d1[l2g_sf[atmi1[2]] * 3 + n3] += t1x * t1m[2][n3];
4420 d1[l2g_sf[atmi1[3]] * 3 + n3] += t1x * t1m[3][n3];
4421 d1[l2g_sf[atmi2[0]] * 3 + n3] += t1x * t1m[4][n3];
4422 d1[l2g_sf[atmi2[1]] * 3 + n3] += t1x * t1m[5][n3];
4423 d1[l2g_sf[atmi2[2]] * 3 + n3] += t1x * t1m[6][n3];
4424 d1[l2g_sf[atmi2[3]] * 3 + n3] += t1x * t1m[7][n3];
4425 }
4426 }
4427 }
4428 }
4429 }
4430
4431 // Richmond TJ : "Solvent Accessible Surface Area and Extended Volume in Proteins"
4432 // J. Mol. Biol 178, 63-89 (1984)
4433
4434 // Fraczkiewicz R, Braun W : "Exact and Efficient Analytical Calculation of the Accessible
4435 // Surface Areas and Their Gradients for Macromolecules" J. Comp. Chem 19, 319-333, (1998)
4436
4437 // Weiser J, Weiser AA, Shenkin PS, Still WC : "Neigbor-List Reduction: Optimization for
4438 // Computation of Molecular van der Waals and Solvent-Accessible Surface Areas"
4439 // J. Comp. Chem. 19, 797-808, (1998)
4440
4441 // Do Carmo MP : "Differential Geometry of Curves and Surfaces" Prentice Hall Inc., 1976
4442
ComputeNBT3(i32u p1)4443 void eng1_sf::ComputeNBT3(i32u p1)
4444 {
4445 i32s limit = GetSetup()->GetSFAtomCount() - num_solvent;
4446 for (i32u n1 = 0;n1 < nbt1_vector.size();n1++)
4447 {
4448 i32s * atmi = nbt1_vector[n1].atmi;
4449
4450 if (atmi[0] < limit && atmi[1] < limit)
4451 {
4452 f64 t1a[3]; f64 t1b = 0.0;
4453 for (i32s n2 = 0;n2 < 3;n2++)
4454 {
4455 f64 t2a = crd[l2g_sf[atmi[0]] * 3 + n2];
4456 f64 t2b = crd[l2g_sf[atmi[1]] * 3 + n2];
4457
4458 t1a[n2] = t2a - t2b;
4459 t1b += t1a[n2] * t1a[n2];
4460 }
4461
4462 f64 t1c = sqrt(t1b);
4463
4464 bool first = (atmi[0] > atmi[1]);
4465 dist2[dist1[atmi[first]] + (atmi[!first] - atmi[first]) - 1] = t1c;
4466
4467 // 1st layer...
4468 // 1st layer...
4469 // 1st layer...
4470
4471 // none of the SASA terms should be skipped at 1st layer -> no test...
4472
4473 if (t1c < (vdwr1[0][atmi[0]] + vdwr1[0][atmi[1]]))
4474 {
4475 nbt3_nl[0][atmi[0]].index[nbt3_nl[0][atmi[0]].index_count++] = atmi[1];
4476 if (nbt3_nl[0][atmi[0]].index_count >= size_nl[0]) { cout << "BUG: NL overflow 1b!!!" << endl; exit(EXIT_FAILURE); }
4477
4478 nbt3_nl[0][atmi[1]].index[nbt3_nl[0][atmi[1]].index_count++] = atmi[0];
4479 if (nbt3_nl[0][atmi[1]].index_count >= size_nl[0]) { cout << "BUG: NL overflow 1b!!!" << endl; exit(EXIT_FAILURE); }
4480 }
4481
4482 // 2nd layer...
4483 // 2nd layer...
4484 // 2nd layer...
4485
4486 bool test2a = (nbt3_nl[1][atmi[0]].index != NULL);
4487 if (test2a && t1c < (vdwr1[1][atmi[0]] + vdwr1[0][atmi[1]]) && t1c > (vdwr1[1][atmi[0]] - vdwr1[0][atmi[1]]))
4488 {
4489 nbt3_nl[1][atmi[0]].index[nbt3_nl[1][atmi[0]].index_count++] = atmi[1];
4490 if (nbt3_nl[1][atmi[0]].index_count >= size_nl[1]) { cout << "BUG: NL overflow 2b!!!" << endl; exit(EXIT_FAILURE); }
4491 }
4492
4493 bool test2b = (nbt3_nl[1][atmi[1]].index != NULL);
4494 if (test2b && t1c < (vdwr1[0][atmi[0]] + vdwr1[1][atmi[1]]) && t1c > (vdwr1[1][atmi[1]] - vdwr1[0][atmi[0]]))
4495 {
4496 nbt3_nl[1][atmi[1]].index[nbt3_nl[1][atmi[1]].index_count++] = atmi[0];
4497 if (nbt3_nl[1][atmi[1]].index_count >= size_nl[1]) { cout << "BUG: NL overflow 2b!!!" << endl; exit(EXIT_FAILURE); }
4498 }
4499
4500 // 3rd layer...
4501 // 3rd layer...
4502 // 3rd layer...
4503
4504 bool test3a = (nbt3_nl[2][atmi[0]].index != NULL);
4505 if (test3a && t1c < (vdwr1[2][atmi[0]] + vdwr1[0][atmi[1]]) && t1c > (vdwr1[2][atmi[0]] - vdwr1[0][atmi[1]]))
4506 {
4507 nbt3_nl[2][atmi[0]].index[nbt3_nl[2][atmi[0]].index_count++] = atmi[1];
4508 if (nbt3_nl[2][atmi[0]].index_count >= size_nl[2]) { cout << "BUG: NL overflow 3b!!!" << endl; exit(EXIT_FAILURE); }
4509 }
4510
4511 bool test3b = (nbt3_nl[2][atmi[1]].index != NULL);
4512 if (test3b && t1c < (vdwr1[0][atmi[0]] + vdwr1[2][atmi[1]]) && t1c > (vdwr1[2][atmi[1]] - vdwr1[0][atmi[0]]))
4513 {
4514 nbt3_nl[2][atmi[1]].index[nbt3_nl[2][atmi[1]].index_count++] = atmi[0];
4515 if (nbt3_nl[2][atmi[1]].index_count >= size_nl[2]) { cout << "BUG: NL overflow 3b!!!" << endl; exit(EXIT_FAILURE); }
4516 }
4517 }
4518
4519 }
4520
4521 energy_nbt3a = 0.0;
4522 energy_nbt3b = 0.0;
4523
4524 for (i32u layer = 0;layer < LAYERS;layer++)
4525 {
4526 sf_nbt3_nl * nlist = nbt3_nl[layer];
4527 // if (!nlist) continue; // skip if no neighbor list; should exist for all layers now!!!
4528
4529 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount() - num_solvent;n1++)
4530 {
4531 sf_nbt3_nd ndt[MAX_SIZE_NL];
4532 for (i32s n2 = 0;n2 < nlist[n1].index_count;n2++)
4533 {
4534 ndt[n2].index = nlist[n1].index[n2];
4535 i32s atmi[2] = { n1, ndt[n2].index }; bool first = (atmi[0] > atmi[1]);
4536 ndt[n2].distance = dist2[dist1[atmi[first]] + (atmi[!first] - atmi[first]) - 1];
4537 }
4538
4539 sort(ndt, ndt + nlist[n1].index_count);
4540 i32s n_count = 0; i32s nt[SIZE_NT];
4541
4542 // neighbor-list reduction... THIS WON'T WORK IF SOME BT1/NBT1-TERMS ARE LEFT OUT!!!
4543 // neighbor-list reduction... THIS WON'T WORK IF SOME BT1/NBT1-TERMS ARE LEFT OUT!!!
4544 // neighbor-list reduction... THIS WON'T WORK IF SOME BT1/NBT1-TERMS ARE LEFT OUT!!!
4545
4546 // test this against a slow-but-simple implementation?!?!?!
4547
4548 for (i32s n2 = 0;n2 < nlist[n1].index_count;n2++)
4549 {
4550 i32s ind1 = ndt[n2].index;
4551 f64 dij = ndt[n2].distance;
4552
4553 bool flag = true;
4554
4555 for (i32s n3 = n2 + 1;n3 < nlist[n1].index_count;n3++)
4556 {
4557 i32s ind2 = ndt[n3].index;
4558
4559 i32s atmi[2] = { ind1, ind2 }; bool first = (atmi[0] > atmi[1]);
4560 f64 djk = dist2[dist1[atmi[first]] + (atmi[!first] - atmi[first]) - 1];
4561
4562 if (djk > dij) continue;
4563
4564 f64 dij2 = dij * dij; f64 djk2 = djk * djk;
4565 f64 dik = ndt[n3].distance; f64 dik2 = dik * dik;
4566
4567 // here dij and dik both represent distances which should never be
4568 // very close to zero (if LJ-terms work as they should) -> no checking
4569
4570 f64 ca = (vdwr2[layer][n1] + dij2 - vdwr2[0][ind1]) / (2.0 * vdwr1[layer][n1] * dij);
4571 f64 cb = (vdwr2[layer][n1] + dik2 - vdwr2[0][ind2]) / (2.0 * vdwr1[layer][n1] * dik);
4572 f64 cg = (dij2 + dik2 - djk2) / (2.0 * dij * dik);
4573
4574 f64 sa2 = 1.0 - ca * ca;
4575 f64 sg2 = 1.0 - cg * cg;
4576
4577 f64 dc = sa2 * sg2;
4578 if (dc < 0.0) dc = 0.0; // domain check...
4579
4580 if (cb < ca * cg - sqrt(dc))
4581 {
4582 flag = false;
4583 break;
4584 }
4585 }
4586
4587 if (flag)
4588 {
4589 nt[n_count++] = ind1;
4590 if (n_count >= SIZE_NT) { cout << "BUG: NT overflow!!!" << endl; exit(EXIT_FAILURE); }
4591 }
4592 }
4593
4594 i32s coi_count = 0; sf_nbt3_coi coit[SIZE_COI];
4595
4596 // next we will create the coi-table...
4597 // next we will create the coi-table...
4598 // next we will create the coi-table...
4599
4600 for (i32s n2 = 0;n2 < n_count;n2++)
4601 {
4602 coit[coi_count].index = nt[n2]; coit[coi_count].flag = false;
4603
4604 f64 t1a[3]; f64 t1b = 0.0;
4605 for (i32s n3 = 0;n3 < 3;n3++)
4606 {
4607 f64 t9a = crd[l2g_sf[n1] * 3 + n3];
4608 f64 t9b = crd[l2g_sf[coit[coi_count].index] * 3 + n3];
4609
4610 t1a[n3] = t9b - t9a;
4611 t1b += t1a[n3] * t1a[n3];
4612 }
4613
4614 f64 t1c = sqrt(t1b);
4615 coit[coi_count].dist = t1c;
4616
4617 // also t1c is a distance which should never be close to zero -> no checking
4618
4619 f64 t2a[3];
4620 for (i32s n3 = 0;n3 < 3;n3++)
4621 {
4622 t2a[n3] = t1a[n3] / t1c;
4623 coit[coi_count].dv[n3] = t2a[n3];
4624 }
4625
4626 coit[coi_count].g = (t1b + vdwr2[layer][n1] - vdwr2[0][coit[coi_count].index]) / (2.0 * t1c);
4627 coit[coi_count].ct = coit[coi_count].g / vdwr1[layer][n1];
4628
4629 if (p1 > 0 && use_implicit_solvent)
4630 {
4631 for (i32s n3 = 0;n3 < 3;n3++)
4632 {
4633 for (i32s n4 = 0;n4 < 3;n4++)
4634 {
4635 f64 t9a = t2a[n3] * t2a[n4]; f64 t9b;
4636 if (n3 != n4) t9b = -t9a; else t9b = 1.0 - t9a;
4637 coit[coi_count].ddv[n3][n4] = t9b / t1c;
4638 }
4639 }
4640
4641 f64 t3a = (t1c - coit[coi_count].g) / t1c;
4642 coit[coi_count].dg[0] = t3a * coit[coi_count].dv[0];
4643 coit[coi_count].dg[1] = t3a * coit[coi_count].dv[1];
4644 coit[coi_count].dg[2] = t3a * coit[coi_count].dv[2];
4645
4646 coit[coi_count].dct[0] = coit[coi_count].dg[0] / vdwr1[layer][n1];
4647 coit[coi_count].dct[1] = coit[coi_count].dg[1] / vdwr1[layer][n1];
4648 coit[coi_count].dct[2] = coit[coi_count].dg[2] / vdwr1[layer][n1];
4649 }
4650
4651 coit[coi_count++].ipd_count = 0;
4652 if (coi_count >= SIZE_COI) { cout << "BUG: COI overflow!!!" << endl; exit(EXIT_FAILURE); }
4653 }
4654
4655 i32s ips_total_count = 0;
4656 i32s ips_count = 0; sf_nbt3_ips ipst[SIZE_IPS];
4657
4658 // next we will create the ips-table...
4659 // next we will create the ips-table...
4660 // next we will create the ips-table...
4661
4662 for (i32s n2 = 0;n2 < coi_count - 1;n2++)
4663 {
4664 for (i32s n3 = n2 + 1;n3 < coi_count;n3++)
4665 {
4666 f64 t1a[3];
4667 t1a[0] = coit[n2].dv[0] * coit[n3].dv[0];
4668 t1a[1] = coit[n2].dv[1] * coit[n3].dv[1];
4669 t1a[2] = coit[n2].dv[2] * coit[n3].dv[2];
4670
4671 f64 t1b = t1a[0] + t1a[1] + t1a[2]; // cos phi
4672
4673 if (t1b < -1.0) t1b = -1.0; // domain check...
4674 if (t1b > +1.0) t1b = +1.0; // domain check...
4675
4676 f64 t1c = 1.0 - t1b * t1b; // sin^2 phi
4677 if (t1c < LOWLIMIT) t1c = LOWLIMIT;
4678
4679 f64 t2a = (coit[n2].g - coit[n3].g * t1b) / t1c; // tau_kj
4680 f64 t2b = (coit[n3].g - coit[n2].g * t1b) / t1c; // tau_jk
4681
4682 f64 t2c = vdwr2[layer][n1] - coit[n2].g * t2a - coit[n3].g * t2b; // gamma^2
4683 if (t2c < LOWLIMIT) continue; // these will not intercept...
4684
4685 ips_total_count++;
4686 coit[n2].flag = true;
4687 coit[n3].flag = true;
4688
4689 f64 t3a[3]; // eta
4690 t3a[0] = coit[n2].dv[0] * t2a + coit[n3].dv[0] * t2b;
4691 t3a[1] = coit[n2].dv[1] * t2a + coit[n3].dv[1] * t2b;
4692 t3a[2] = coit[n2].dv[2] * t2a + coit[n3].dv[2] * t2b;
4693
4694 f64 t1d = sqrt(t1c); // sin phi
4695
4696 f64 t3b[3]; // omega
4697 t3b[0] = (coit[n2].dv[1] * coit[n3].dv[2] - coit[n2].dv[2] * coit[n3].dv[1]) / t1d;
4698 t3b[1] = (coit[n2].dv[2] * coit[n3].dv[0] - coit[n2].dv[0] * coit[n3].dv[2]) / t1d;
4699 t3b[2] = (coit[n2].dv[0] * coit[n3].dv[1] - coit[n2].dv[1] * coit[n3].dv[0]) / t1d;
4700
4701 f64 t2d = sqrt(t2c); // gamma
4702
4703 for (i32s n4 = 0;n4 < 3;n4++)
4704 {
4705 f64 t9a = t3b[n4] * t2d;
4706 ipst[ips_count].ipv[0][n4] = t3a[n4] - t9a;
4707 ipst[ips_count].ipv[1][n4] = t3a[n4] + t9a;
4708 }
4709
4710 // skip those intersection points that fall inside any other sphere...
4711 // SKIP ALSO IF EQUAL DISTANCE??? i.e. compare using "<" or "<=" ???
4712
4713 bool skip_both = false;
4714 bool skip[2] = { false, false };
4715 for (i32s n4 = 0;n4 < n_count;n4++)
4716 {
4717 i32s n5 = nt[n4];
4718 if (n5 == coit[n2].index || n5 == coit[n3].index) continue;
4719
4720 f64 t9a[3];
4721 t9a[0] = (crd[l2g_sf[n1] * 3 + 0] + ipst[ips_count].ipv[0][0]) - crd[l2g_sf[n5] * 3 + 0];
4722 t9a[1] = (crd[l2g_sf[n1] * 3 + 1] + ipst[ips_count].ipv[0][1]) - crd[l2g_sf[n5] * 3 + 1];
4723 t9a[2] = (crd[l2g_sf[n1] * 3 + 2] + ipst[ips_count].ipv[0][2]) - crd[l2g_sf[n5] * 3 + 2];
4724 f64 t9b = t9a[0] * t9a[0] + t9a[1] * t9a[1] + t9a[2] * t9a[2];
4725 if (t9b < vdwr2[0][n5]) skip[0] = true;
4726
4727 f64 t9c[3];
4728 t9c[0] = (crd[l2g_sf[n1] * 3 + 0] + ipst[ips_count].ipv[1][0]) - crd[l2g_sf[n5] * 3 + 0];
4729 t9c[1] = (crd[l2g_sf[n1] * 3 + 1] + ipst[ips_count].ipv[1][1]) - crd[l2g_sf[n5] * 3 + 1];
4730 t9c[2] = (crd[l2g_sf[n1] * 3 + 2] + ipst[ips_count].ipv[1][2]) - crd[l2g_sf[n5] * 3 + 2];
4731 f64 t9d = t9c[0] * t9c[0] + t9c[1] * t9c[1] + t9c[2] * t9c[2];
4732 if (t9d < vdwr2[0][n5]) skip[1] = true;
4733
4734 skip_both = (skip[0] && skip[1]);
4735 if (skip_both) break;
4736 }
4737
4738 if (skip_both) continue; // overwrite this one...
4739
4740 ipst[ips_count].coi[0] = n2;
4741 ipst[ips_count].coi[1] = n3;
4742
4743 if (!skip[0])
4744 {
4745 coit[n2].AddIPD(ipst[ips_count].ipv[0], ips_count);
4746 coit[n3].AddIPD(ipst[ips_count].ipv[0], ips_count | ORDER_FLAG);
4747 }
4748
4749 if (!skip[1])
4750 {
4751 coit[n2].AddIPD(ipst[ips_count].ipv[1], ips_count | INDEX_FLAG | ORDER_FLAG);
4752 coit[n3].AddIPD(ipst[ips_count].ipv[1], ips_count | INDEX_FLAG);
4753 }
4754
4755 if (p1 > 0 && use_implicit_solvent)
4756 {
4757 f64 t1f[3]; // d(cos phi) / dXk
4758 t1f[0] = (coit[n3].dv[0] - t1b * coit[n2].dv[0]) / coit[n2].dist;
4759 t1f[1] = (coit[n3].dv[1] - t1b * coit[n2].dv[1]) / coit[n2].dist;
4760 t1f[2] = (coit[n3].dv[2] - t1b * coit[n2].dv[2]) / coit[n2].dist;
4761
4762 f64 t1g[3]; // d(cos phi) / dXj
4763 t1g[0] = (coit[n2].dv[0] - t1b * coit[n3].dv[0]) / coit[n3].dist;
4764 t1g[1] = (coit[n2].dv[1] - t1b * coit[n3].dv[1]) / coit[n3].dist;
4765 t1g[2] = (coit[n2].dv[2] - t1b * coit[n3].dv[2]) / coit[n3].dist;
4766
4767 f64 t2e[3]; // d(tau_kj) / dXk
4768 t2e[0] = (t1f[0] * (2.0 * t2a * t1b - coit[n3].g) + coit[n2].dg[0]) / t1c;
4769 t2e[1] = (t1f[1] * (2.0 * t2a * t1b - coit[n3].g) + coit[n2].dg[1]) / t1c;
4770 t2e[2] = (t1f[2] * (2.0 * t2a * t1b - coit[n3].g) + coit[n2].dg[2]) / t1c;
4771
4772 f64 t2f[3]; // d(tau_kj) / dXj
4773 t2f[0] = (t1g[0] * (2.0 * t2a * t1b - coit[n3].g) - t1b * coit[n3].dg[0]) / t1c;
4774 t2f[1] = (t1g[1] * (2.0 * t2a * t1b - coit[n3].g) - t1b * coit[n3].dg[1]) / t1c;
4775 t2f[2] = (t1g[2] * (2.0 * t2a * t1b - coit[n3].g) - t1b * coit[n3].dg[2]) / t1c;
4776
4777 f64 t2g[3]; // d(tau_jk) / dXk
4778 t2g[0] = (t1f[0] * (2.0 * t2b * t1b - coit[n2].g) - t1b * coit[n2].dg[0]) / t1c;
4779 t2g[1] = (t1f[1] * (2.0 * t2b * t1b - coit[n2].g) - t1b * coit[n2].dg[1]) / t1c;
4780 t2g[2] = (t1f[2] * (2.0 * t2b * t1b - coit[n2].g) - t1b * coit[n2].dg[2]) / t1c;
4781
4782 f64 t2h[3]; // d(tau_jk) / dXj
4783 t2h[0] = (t1g[0] * (2.0 * t2b * t1b - coit[n2].g) + coit[n3].dg[0]) / t1c;
4784 t2h[1] = (t1g[1] * (2.0 * t2b * t1b - coit[n2].g) + coit[n3].dg[1]) / t1c;
4785 t2h[2] = (t1g[2] * (2.0 * t2b * t1b - coit[n2].g) + coit[n3].dg[2]) / t1c;
4786
4787 f64 t3c[3][3]; // d(eta) / dXk
4788 f64 t3d[3][3]; // d(eta) / dXj
4789
4790 for (i32s n4 = 0;n4 < 3;n4++)
4791 {
4792 for (i32s n5 = 0;n5 < 3;n5++)
4793 {
4794 f64 t9a = coit[n2].dv[n5]; f64 t9b = coit[n3].dv[n5];
4795 t3c[n4][n5] = t9a * t2e[n4] + t9b * t2g[n4] + t2a * coit[n2].ddv[n4][n5];
4796 t3d[n4][n5] = t9a * t2f[n4] + t9b * t2h[n4] + t2b * coit[n3].ddv[n4][n5];
4797 }
4798 }
4799
4800 f64 t3e[3][3]; // d(omega) / dXk
4801 f64 t3f[3][3]; // d(omega) / dXj
4802
4803 for (i32s n4 = 0;n4 < 3;n4++)
4804 {
4805 for (i32s n5 = 0;n5 < 3;n5++)
4806 {
4807 t3e[n4][n5] = t1b * t3b[n5] * t1f[n4] / t1c;
4808 t3f[n4][n5] = t1b * t3b[n5] * t1g[n4] / t1c;
4809 }
4810
4811 t3e[n4][0] += (coit[n2].ddv[n4][1] * coit[n3].dv[2] - coit[n2].ddv[n4][2] * coit[n3].dv[1]) / t1d;
4812 t3e[n4][1] += (coit[n2].ddv[n4][2] * coit[n3].dv[0] - coit[n2].ddv[n4][0] * coit[n3].dv[2]) / t1d;
4813 t3e[n4][2] += (coit[n2].ddv[n4][0] * coit[n3].dv[1] - coit[n2].ddv[n4][1] * coit[n3].dv[0]) / t1d;
4814
4815 t3f[n4][0] += (coit[n2].dv[1] * coit[n3].ddv[n4][2] - coit[n2].dv[2] * coit[n3].ddv[n4][1]) / t1d;
4816 t3f[n4][1] += (coit[n2].dv[2] * coit[n3].ddv[n4][0] - coit[n2].dv[0] * coit[n3].ddv[n4][2]) / t1d;
4817 t3f[n4][2] += (coit[n2].dv[0] * coit[n3].ddv[n4][1] - coit[n2].dv[1] * coit[n3].ddv[n4][0]) / t1d;
4818 }
4819
4820 f64 t2i[3]; // d(gamma) / dXk
4821 t2i[0] = -(coit[n2].g * t2e[0] + t2a * coit[n2].dg[0] + coit[n3].g * t2g[0]) / (2.0 * t2d);
4822 t2i[1] = -(coit[n2].g * t2e[1] + t2a * coit[n2].dg[1] + coit[n3].g * t2g[1]) / (2.0 * t2d);
4823 t2i[2] = -(coit[n2].g * t2e[2] + t2a * coit[n2].dg[2] + coit[n3].g * t2g[2]) / (2.0 * t2d);
4824
4825 f64 t2j[3]; // d(gamma) / dXj
4826 t2j[0] = -(coit[n2].g * t2f[0] + coit[n3].g * t2h[0] + t2b * coit[n3].dg[0]) / (2.0 * t2d);
4827 t2j[1] = -(coit[n2].g * t2f[1] + coit[n3].g * t2h[1] + t2b * coit[n3].dg[1]) / (2.0 * t2d);
4828 t2j[2] = -(coit[n2].g * t2f[2] + coit[n3].g * t2h[2] + t2b * coit[n3].dg[2]) / (2.0 * t2d);
4829
4830 // the final result is derivatives for points dipv[2][2][3][3].
4831 // indexes are as follows: [point][atom][variable][xyz].
4832
4833 for (i32s n4 = 0;n4 < 3;n4++)
4834 {
4835 for (i32s n5 = 0;n5 < 3;n5++)
4836 {
4837 ipst[ips_count].dipv[0][0][n4][n5] = t3c[n4][n5];
4838 ipst[ips_count].dipv[0][1][n4][n5] = t3d[n4][n5];
4839 ipst[ips_count].dipv[1][0][n4][n5] = t3c[n4][n5];
4840 ipst[ips_count].dipv[1][1][n4][n5] = t3d[n4][n5];
4841 }
4842
4843 for (i32s n5 = 0;n5 < 3;n5++)
4844 {
4845 f64 t9a = t3b[n5] * t2i[n4] + t2d * t3e[n4][n5];
4846 f64 t9b = t3b[n5] * t2j[n4] + t2d * t3f[n4][n5];
4847
4848 ipst[ips_count].dipv[0][0][n4][n5] -= t9a;
4849 ipst[ips_count].dipv[0][1][n4][n5] -= t9b;
4850 ipst[ips_count].dipv[1][0][n4][n5] += t9a;
4851 ipst[ips_count].dipv[1][1][n4][n5] += t9b;
4852 }
4853 }
4854 }
4855
4856 ips_count++;
4857 if (ips_count >= SIZE_IPS) { cout << "BUG: IPS overflow!!!" << endl; exit(EXIT_FAILURE); }
4858 }
4859 }
4860
4861 i32s arc_count = 0; sf_nbt3_arc arct[SIZE_ARC];
4862
4863 // next we will create the arc-table...
4864 // next we will create the arc-table...
4865 // next we will create the arc-table...
4866
4867 for (i32s n2 = 0;n2 < coi_count;n2++)
4868 {
4869 f64 t1z = vdwr2[layer][n1] - coit[n2].g * coit[n2].g;
4870 if (t1z < 0.0) t1z = 0.0; // domain check...
4871
4872 f64 t1a = sqrt(t1z);
4873 if (t1a < LOWLIMIT) t1a = LOWLIMIT;
4874
4875 sort(coit[n2].ipdt, coit[n2].ipdt + coit[n2].ipd_count);
4876
4877 for (i32s n3 = 0;n3 < coit[n2].ipd_count;n3++)
4878 {
4879 if (coit[n2].ipdt[n3].ipdata & ORDER_FLAG) continue;
4880 i32s n4 = n3 + 1; if (n4 == coit[n2].ipd_count) n4 = 0;
4881 if (!(coit[n2].ipdt[n4].ipdata & ORDER_FLAG)) continue;
4882
4883 arct[arc_count].coi = n2; arct[arc_count].flag = false;
4884
4885 arct[arc_count].ipdata[0] = (coit[n2].ipdt[n3].ipdata & ~ORDER_FLAG);
4886 arct[arc_count].ipdata[1] = (coit[n2].ipdt[n4].ipdata & ~ORDER_FLAG);
4887
4888 i32s i1a = (arct[arc_count].ipdata[0] & FLAG_MASK);
4889 bool i1b = (arct[arc_count].ipdata[0] & INDEX_FLAG ? 1 : 0);
4890
4891 i32s i2a = (arct[arc_count].ipdata[1] & FLAG_MASK);
4892 bool i2b = (arct[arc_count].ipdata[1] & INDEX_FLAG ? 1 : 0);
4893
4894 arct[arc_count].index[0][0] = coit[ipst[i1a].coi[i1b]].index;
4895 arct[arc_count].index[0][1] = coit[ipst[i1a].coi[!i1b]].index;
4896
4897 arct[arc_count].index[1][0] = coit[ipst[i2a].coi[!i2b]].index;
4898 arct[arc_count].index[1][1] = coit[ipst[i2a].coi[i2b]].index;
4899
4900 // let's compute the tangent vectors...
4901
4902 f64 * ref1 = ipst[i1a].ipv[i1b];
4903 arct[arc_count].tv[0][0] = (ref1[1] * coit[n2].dv[2] - ref1[2] * coit[n2].dv[1]) / t1a;
4904 arct[arc_count].tv[0][1] = (ref1[2] * coit[n2].dv[0] - ref1[0] * coit[n2].dv[2]) / t1a;
4905 arct[arc_count].tv[0][2] = (ref1[0] * coit[n2].dv[1] - ref1[1] * coit[n2].dv[0]) / t1a;
4906
4907 f64 * ref2 = ipst[i2a].ipv[i2b];
4908 arct[arc_count].tv[1][0] = (ref2[1] * coit[n2].dv[2] - ref2[2] * coit[n2].dv[1]) / t1a;
4909 arct[arc_count].tv[1][1] = (ref2[2] * coit[n2].dv[0] - ref2[0] * coit[n2].dv[2]) / t1a;
4910 arct[arc_count].tv[1][2] = (ref2[0] * coit[n2].dv[1] - ref2[1] * coit[n2].dv[0]) / t1a;
4911
4912 if (p1 > 0 && use_implicit_solvent)
4913 {
4914 for (i32s n4 = 0;n4 < 3;n4++)
4915 {
4916 f64 t9a = coit[n2].g * coit[n2].dg[n4] / t1a;
4917 for (i32s n5 = 0;n5 < 3;n5++)
4918 {
4919 arct[arc_count].dtv[0][0][n4][n5] = t9a * arct[arc_count].tv[0][n5];
4920 arct[arc_count].dtv[1][0][n4][n5] = t9a * arct[arc_count].tv[1][n5];
4921 }
4922
4923 f64 * ref1a = ipst[i1a].dipv[i1b][i1b][n4]; // d(P1) / dXk
4924 arct[arc_count].dtv[0][0][n4][0] += ref1a[1] * coit[n2].dv[2] - ref1a[2] * coit[n2].dv[1];
4925 arct[arc_count].dtv[0][0][n4][1] += ref1a[2] * coit[n2].dv[0] - ref1a[0] * coit[n2].dv[2];
4926 arct[arc_count].dtv[0][0][n4][2] += ref1a[0] * coit[n2].dv[1] - ref1a[1] * coit[n2].dv[0];
4927
4928 f64 * ref1b = ipst[i2a].dipv[i2b][!i2b][n4]; // d(P2) / dXk
4929 arct[arc_count].dtv[1][0][n4][0] += ref1b[1] * coit[n2].dv[2] - ref1b[2] * coit[n2].dv[1];
4930 arct[arc_count].dtv[1][0][n4][1] += ref1b[2] * coit[n2].dv[0] - ref1b[0] * coit[n2].dv[2];
4931 arct[arc_count].dtv[1][0][n4][2] += ref1b[0] * coit[n2].dv[1] - ref1b[1] * coit[n2].dv[0];
4932
4933 f64 * ref2a = ipst[i1a].ipv[i1b];
4934 arct[arc_count].dtv[0][0][n4][0] += ref2a[1] * coit[n2].ddv[n4][2] - ref2a[2] * coit[n2].ddv[n4][1];
4935 arct[arc_count].dtv[0][0][n4][1] += ref2a[2] * coit[n2].ddv[n4][0] - ref2a[0] * coit[n2].ddv[n4][2];
4936 arct[arc_count].dtv[0][0][n4][2] += ref2a[0] * coit[n2].ddv[n4][1] - ref2a[1] * coit[n2].ddv[n4][0];
4937
4938 f64 * ref2b = ipst[i2a].ipv[i2b];
4939 arct[arc_count].dtv[1][0][n4][0] += ref2b[1] * coit[n2].ddv[n4][2] - ref2b[2] * coit[n2].ddv[n4][1];
4940 arct[arc_count].dtv[1][0][n4][1] += ref2b[2] * coit[n2].ddv[n4][0] - ref2b[0] * coit[n2].ddv[n4][2];
4941 arct[arc_count].dtv[1][0][n4][2] += ref2b[0] * coit[n2].ddv[n4][1] - ref2b[1] * coit[n2].ddv[n4][0];
4942
4943 for (i32s n5 = 0;n5 < 3;n5++)
4944 {
4945 arct[arc_count].dtv[0][0][n4][n5] /= t1a;
4946 arct[arc_count].dtv[1][0][n4][n5] /= t1a;
4947 }
4948
4949 f64 * ref3a = ipst[i1a].dipv[i1b][!i1b][n4]; // d(P1) / dXj
4950 arct[arc_count].dtv[0][1][n4][0] = (ref3a[1] * coit[n2].dv[2] - ref3a[2] * coit[n2].dv[1]) / t1a;
4951 arct[arc_count].dtv[0][1][n4][1] = (ref3a[2] * coit[n2].dv[0] - ref3a[0] * coit[n2].dv[2]) / t1a;
4952 arct[arc_count].dtv[0][1][n4][2] = (ref3a[0] * coit[n2].dv[1] - ref3a[1] * coit[n2].dv[0]) / t1a;
4953
4954 f64 * ref3b = ipst[i2a].dipv[i2b][i2b][n4]; // d(P2) / dXj
4955 arct[arc_count].dtv[1][1][n4][0] = (ref3b[1] * coit[n2].dv[2] - ref3b[2] * coit[n2].dv[1]) / t1a;
4956 arct[arc_count].dtv[1][1][n4][1] = (ref3b[2] * coit[n2].dv[0] - ref3b[0] * coit[n2].dv[2]) / t1a;
4957 arct[arc_count].dtv[1][1][n4][2] = (ref3b[0] * coit[n2].dv[1] - ref3b[1] * coit[n2].dv[0]) / t1a;
4958 }
4959 }
4960
4961 arc_count++;
4962 if (arc_count >= SIZE_ARC) { cout << "BUG: ARC overflow!!!" << endl; exit(EXIT_FAILURE); }
4963 }
4964 }
4965
4966 // all cases will pass through this point!!!
4967 // all cases will pass through this point!!!
4968 // all cases will pass through this point!!!
4969
4970 f64 area;
4971 if (!arc_count)
4972 {
4973 if (ips_total_count)
4974 {
4975 solv_exp[layer][n1] = 0.0; // fully buried...
4976 continue; // fully buried...
4977 }
4978 else area = 4.0 * M_PI;
4979 }
4980 else
4981 {
4982 area = 0.0;
4983 i32s arc_counter = 0;
4984
4985 do
4986 {
4987 i32s prev; i32s curr = 0;
4988 while (arct[curr].flag)
4989 {
4990 curr++;
4991 if (curr == arc_count)
4992 {
4993 cout << "area_panic: can't find the first arc!!!" << endl;
4994 goto area_panic;
4995 }
4996 }
4997
4998 i32s first = curr;
4999
5000 f64 sum1 = 0.0;
5001 f64 sum2 = 0.0;
5002
5003 while (true)
5004 {
5005 i32s coi = arct[curr].coi;
5006
5007 f64 t1a[3];
5008 t1a[0] = arct[curr].tv[1][1] * arct[curr].tv[0][2] - arct[curr].tv[1][2] * arct[curr].tv[0][1];
5009 t1a[1] = arct[curr].tv[1][2] * arct[curr].tv[0][0] - arct[curr].tv[1][0] * arct[curr].tv[0][2];
5010 t1a[2] = arct[curr].tv[1][0] * arct[curr].tv[0][1] - arct[curr].tv[1][1] * arct[curr].tv[0][0];
5011
5012 f64 t1b[3];
5013 t1b[0] = coit[coi].dv[0] * t1a[0];
5014 t1b[1] = coit[coi].dv[1] * t1a[1];
5015 t1b[2] = coit[coi].dv[2] * t1a[2];
5016
5017 f64 t1c = (t1b[0] + t1b[1] + t1b[2] < 0.0 ? -1.0 : +1.0);
5018
5019 f64 t2a[3];
5020 t2a[0] = arct[curr].tv[0][0] * arct[curr].tv[1][0];
5021 t2a[1] = arct[curr].tv[0][1] * arct[curr].tv[1][1];
5022 t2a[2] = arct[curr].tv[0][2] * arct[curr].tv[1][2];
5023
5024 f64 t2b = t2a[0] + t2a[1] + t2a[2];
5025
5026 if (t2b < -1.0) t2b = -1.0; // domain check...
5027 if (t2b > +1.0) t2b = +1.0; // domain check...
5028
5029 f64 t2c = (1.0 - t1c) * M_PI + t1c * acos(t2b);
5030 sum1 += t2c * coit[coi].ct;
5031
5032 if (p1 > 0 && use_implicit_solvent)
5033 {
5034 f64 t2x = fabs(sin(t2c));
5035 if (t2x < LOWLIMIT) t2x = LOWLIMIT;
5036
5037 f64 t2y = -sasaE[layer][n1] * coit[coi].ct * t1c / t2x;
5038
5039 // 1st are same points and 2nd are different ones...
5040 // 1st are same points and 2nd are different ones...
5041 // 1st are same points and 2nd are different ones...
5042
5043 for (i32s n2 = 0;n2 < 3;n2++)
5044 {
5045 f64 t3a[3];
5046 t3a[0] = arct[curr].dtv[0][0][n2][0] * arct[curr].tv[1][0];
5047 t3a[1] = arct[curr].dtv[0][0][n2][1] * arct[curr].tv[1][1];
5048 t3a[2] = arct[curr].dtv[0][0][n2][2] * arct[curr].tv[1][2];
5049 f64 t3b = t3a[0] + t3a[1] + t3a[2];
5050
5051 f64 t3c[3];
5052 t3c[0] = arct[curr].tv[0][0] * arct[curr].dtv[1][0][n2][0];
5053 t3c[1] = arct[curr].tv[0][1] * arct[curr].dtv[1][0][n2][1];
5054 t3c[2] = arct[curr].tv[0][2] * arct[curr].dtv[1][0][n2][2];
5055 f64 t3d = t3c[0] + t3c[1] + t3c[2];
5056
5057 f64 t4a[3];
5058 t4a[0] = arct[curr].dtv[0][1][n2][0] * arct[curr].tv[1][0];
5059 t4a[1] = arct[curr].dtv[0][1][n2][1] * arct[curr].tv[1][1];
5060 t4a[2] = arct[curr].dtv[0][1][n2][2] * arct[curr].tv[1][2];
5061 f64 t4b = t4a[0] + t4a[1] + t4a[2];
5062
5063 f64 t4c[3];
5064 t4c[0] = arct[curr].tv[0][0] * arct[curr].dtv[1][1][n2][0];
5065 t4c[1] = arct[curr].tv[0][1] * arct[curr].dtv[1][1][n2][1];
5066 t4c[2] = arct[curr].tv[0][2] * arct[curr].dtv[1][1][n2][2];
5067 f64 t4d = t4c[0] + t4c[1] + t4c[2];
5068
5069 f64 t3e = t2y * (t3b + t3d) + sasaE[layer][n1] * t2c * coit[coi].dct[n2];
5070 f64 t5a = t2y * t4b; f64 t5b = t2y * t4d;
5071
5072 d1[l2g_sf[arct[curr].index[0][0]] * 3 + n2] += t3e;
5073 d1[l2g_sf[arct[curr].index[0][1]] * 3 + n2] += t5a;
5074 d1[l2g_sf[arct[curr].index[1][1]] * 3 + n2] += t5b;
5075 d1[l2g_sf[n1] * 3 + n2] -= t3e + t5a + t5b;
5076 }
5077 }
5078
5079 prev = curr; curr = 0;
5080 i32u ipd = arct[prev].ipdata[1];
5081 while (true)
5082 {
5083 if (arct[curr].ipdata[0] != ipd) curr++;
5084 else break;
5085
5086 if (curr == arc_count)
5087 {
5088 cout << "area_panic: incomplete set of arcs!!!" << endl;
5089 goto area_panic;
5090 }
5091 }
5092
5093 arc_counter++;
5094 arct[curr].flag = true;
5095
5096 f64 t2d[3];
5097 t2d[0] = arct[prev].tv[1][0] * arct[curr].tv[0][0];
5098 t2d[1] = arct[prev].tv[1][1] * arct[curr].tv[0][1];
5099 t2d[2] = arct[prev].tv[1][2] * arct[curr].tv[0][2];
5100
5101 f64 t2e = t2d[0] + t2d[1] + t2d[2];
5102
5103 if (t2e < -1.0) t2e = -1.0; // domain check...
5104 if (t2e > +1.0) t2e = +1.0; // domain check...
5105
5106 f64 t2f = -acos(t2e); sum2 += t2f;
5107
5108 if (p1 > 0 && use_implicit_solvent)
5109 {
5110 f64 t2x = fabs(sin(t2f));
5111 if (t2x < LOWLIMIT) t2x = LOWLIMIT;
5112
5113 f64 t2y = sasaE[layer][n1] / t2x;
5114
5115 // prev_k = curr_j and prev_j = curr_k !!!
5116 // prev_k = curr_j and prev_j = curr_k !!!
5117 // prev_k = curr_j and prev_j = curr_k !!!
5118
5119 for (i32s n2 = 0;n2 < 3;n2++)
5120 {
5121 f64 t3a[3];
5122 t3a[0] = arct[prev].dtv[1][0][n2][0] * arct[curr].tv[0][0];
5123 t3a[1] = arct[prev].dtv[1][0][n2][1] * arct[curr].tv[0][1];
5124 t3a[2] = arct[prev].dtv[1][0][n2][2] * arct[curr].tv[0][2];
5125 f64 t3b = t3a[0] + t3a[1] + t3a[2];
5126
5127 f64 t3c[3];
5128 t3c[0] = arct[prev].tv[1][0] * arct[curr].dtv[0][1][n2][0];
5129 t3c[1] = arct[prev].tv[1][1] * arct[curr].dtv[0][1][n2][1];
5130 t3c[2] = arct[prev].tv[1][2] * arct[curr].dtv[0][1][n2][2];
5131 f64 t3d = t3c[0] + t3c[1] + t3c[2];
5132
5133 f64 t4a[3];
5134 t4a[0] = arct[prev].dtv[1][1][n2][0] * arct[curr].tv[0][0];
5135 t4a[1] = arct[prev].dtv[1][1][n2][1] * arct[curr].tv[0][1];
5136 t4a[2] = arct[prev].dtv[1][1][n2][2] * arct[curr].tv[0][2];
5137 f64 t4b = t4a[0] + t4a[1] + t4a[2];
5138
5139 f64 t4c[3];
5140 t4c[0] = arct[prev].tv[1][0] * arct[curr].dtv[0][0][n2][0];
5141 t4c[1] = arct[prev].tv[1][1] * arct[curr].dtv[0][0][n2][1];
5142 t4c[2] = arct[prev].tv[1][2] * arct[curr].dtv[0][0][n2][2];
5143 f64 t4d = t4c[0] + t4c[1] + t4c[2];
5144
5145 f64 t3e = t2y * (t3b + t3d);
5146 f64 t4e = t2y * (t4b + t4d);
5147
5148 d1[l2g_sf[arct[prev].index[1][0]] * 3 + n2] += t3e;
5149 d1[l2g_sf[arct[prev].index[1][1]] * 3 + n2] += t4e;
5150 d1[l2g_sf[n1] * 3 + n2] -= t3e + t4e;
5151 }
5152 }
5153
5154 if (curr == first) break;
5155 }
5156
5157 area += 2.0 * M_PI + sum1 + sum2;
5158 } while (arc_counter < arc_count);
5159
5160 // when we have some problems somewhere above (for example, if we have
5161 // an incomplete set of arcs or no arcs at all; these things are possible
5162 // in rare special cases; for example we might have to reject some arcs
5163 // if they contained some singular intermediate values) we will truncate
5164 // the sum and jump right here.
5165
5166 // in this case we will calculate incorrect value for the area, but the
5167 // good news is that the value and the gradient will still be consistent.
5168
5169 // since these cases are very rare, this probably won't make big problems
5170 // in any applications...
5171
5172 area_panic: // we will jump here in all problematic cases...
5173
5174 while (area > 4.0 * M_PI) area -= 4.0 * M_PI;
5175 }
5176
5177 // finally here we will handle the single separate patches...
5178 // finally here we will handle the single separate patches...
5179 // finally here we will handle the single separate patches...
5180
5181 for (i32s n2 = 0;n2 < coi_count;n2++)
5182 {
5183 if (coit[n2].flag) continue;
5184
5185 f64 t1a = 2.0 * M_PI / vdwr1[layer][n1];
5186 area -= t1a * (vdwr1[layer][n1] - coit[n2].g);
5187
5188 if (p1 > 0 && use_implicit_solvent)
5189 {
5190 for (i32s n3 = 0;n3 < 3;n3++)
5191 {
5192 f64 t1b = sasaE[layer][n1] * t1a * coit[n2].dg[n3];
5193
5194 d1[l2g_sf[coit[n2].index] * 3 + n3] += t1b;
5195 d1[l2g_sf[n1] * 3 + n3] -= t1b;
5196 }
5197 }
5198 }
5199
5200 f64 max_area; // this is 4pi precisely, but the values below are practical max values.
5201 switch (layer)
5202 {
5203 case 0: max_area = 8.7; break;
5204 case 1: max_area = 10.0; break;
5205 case 2: max_area = 10.8; break;
5206 default:
5207 cout << "unknown layer!" << endl;
5208 exit(EXIT_FAILURE);
5209 }
5210
5211 solv_exp[layer][n1] = area / max_area;
5212 if (solv_exp[layer][n1] < 0.0) solv_exp[layer][n1] = 0.0; // domain check!!!
5213 if (solv_exp[layer][n1] > 1.0) solv_exp[layer][n1] = 1.0; // domain check!!!
5214
5215 if (!use_implicit_solvent) continue;
5216
5217 f64 value = sasaE[layer][n1] * area;
5218 if (sasaE[layer][n1] > 0.0) energy_nbt3a += value; // positive values...
5219 else energy_nbt3b += value; // negative values...
5220 }
5221 //cout << "layer = " << layer << " energy_nbt3 = " << energy_nbt3 << endl;
5222 }
5223 }
5224
Compute(i32u p1,bool)5225 void eng1_sf::Compute(i32u p1, bool)
5226 {
5227 if (p1 > 0)
5228 {
5229 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount();n1++)
5230 {
5231 d1[l2g_sf[n1] * 3 + 0] = 0.0;
5232 d1[l2g_sf[n1] * 3 + 1] = 0.0;
5233 d1[l2g_sf[n1] * 3 + 2] = 0.0;
5234 }
5235 }
5236
5237 // this is for the surface area term...
5238 // this is for the surface area term...
5239 // this is for the surface area term...
5240
5241 for (i32u n1 = 0;n1 < LAYERS;n1++)
5242 {
5243 for (i32s n2 = 0;n2 < GetSetup()->GetSFAtomCount() - num_solvent;n2++)
5244 {
5245 nbt3_nl[n1][n2].index_count = 0;
5246 }
5247 }
5248
5249 ComputeBT1(p1); // we need this also for the surface area term...
5250 ComputeBT2(p1);
5251 ComputeBT3(p1);
5252 ComputeBT4(p1);
5253 //energy_bt1 = 0.0;
5254 //energy_bt2 = 0.0;
5255 //energy_bt3a = 0.0; energy_bt3b = 0.0;
5256 //energy_bt4a = 0.0; energy_bt4b = 0.0;
5257
5258 ComputeNBT3(p1); // we need this also for determining epsilon in electrostatics...
5259 ComputeNBT2(p1);
5260 ComputeNBT1(p1);
5261 //energy_nbt3a = 0.0; energy_nbt3b = 0.0;
5262 //energy_nbt2a = 0.0; energy_nbt2b = 0.0; energy_nbt2c = 0.0;
5263 //energy_nbt1a = 0.0; energy_nbt1b = 0.0; energy_nbt1c = 0.0;
5264
5265 energy = energy_bt1 + energy_bt2;
5266 energy += energy_bt3a + energy_bt3b;
5267 energy += energy_bt4a + energy_bt4b;
5268 energy += energy_nbt1a + energy_nbt1b + energy_nbt1c;
5269 energy += energy_nbt2a + energy_nbt2b + energy_nbt2c;
5270 energy += energy_nbt3a + energy_nbt3b;
5271
5272 /////////////////////////////////////////////
5273 /////////////////////////////////////////////
5274 // this will print also the components...
5275 // this will print also the components...
5276 // this will print also the components...
5277 /*ostringstream str;
5278 str.setf(ios::fixed); str.precision(8);
5279
5280 str << "B: ";
5281 str << energy_bt1 << " ";
5282 str << energy_bt2 << " ";
5283 str << energy_bt3a << " " << energy_bt3b << " ";
5284 str << energy_bt4a << " " << energy_bt4b << " ";
5285 str << "NB: ";
5286 str << energy_nbt1a << " " << energy_nbt1b << " " << energy_nbt1c << " ";
5287 str << energy_nbt2a << " " << energy_nbt2b << " " << energy_nbt2c << " ";
5288 str << energy_nbt3a << " " << energy_nbt3b << " ";
5289 str << endl;
5290
5291 str << "Energy = " << energy << " kJ/mol + " << constraints << " kJ/mol = ";
5292 str << (energy + constraints) << " kJ/mol" << ends;
5293
5294 cout << str.str().c_str() << endl;*/
5295 /////////////////////////////////////////////
5296 /////////////////////////////////////////////
5297
5298 // for consistency, it seems we have to include constraints in energy...
5299 // for consistency, it seems we have to include constraints in energy...
5300 // for consistency, it seems we have to include constraints in energy...
5301
5302 energy += constraints;
5303 }
5304
5305 // f = sum[(r/a)^-3] = sum[a^3 * r^-3] // now seems to be r^-12
5306 // df/dr = -3 * sum[a^3 * r^-4]
5307
GetVDWSurf(fGL * pp,fGL * dd)5308 fGL eng1_sf::GetVDWSurf(fGL * pp, fGL * dd)
5309 {
5310 fGL vdwsv = 0.0;
5311 if (dd != NULL) dd[0] = dd[1] = dd[2] = 0.0;
5312
5313 atom ** atmtab = GetSetup()->GetSFAtoms();
5314 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount();n1++)
5315 {
5316 fGL tmp1[3]; fGL r2 = 0.0;
5317 const f64 * cdata = & crd[l2g_sf[n1] * 3];
5318 for (i32s n2 = 0;n2 < 3;n2++)
5319 {
5320 tmp1[n2] = pp[n2] - cdata[n2];
5321 r2 += tmp1[n2] * tmp1[n2];
5322 }
5323
5324 if (r2 == 0.0) return +1.0e+35; // numeric_limits<fGL>::max()?!?!?!
5325 fGL r1 = sqrt(r2);
5326
5327 fGL tmp2 = r1 / (atmtab[n1]->vdwr + 0.15); // solvent radius??? 0.15
5328 fGL qqq = tmp2 * tmp2 * tmp2 * tmp2;
5329
5330 fGL tmp3 = 1.0 / (qqq * qqq * qqq);
5331 vdwsv += tmp3;
5332
5333 if (dd != NULL) // sign ??? constant ???
5334 {
5335 for (i32s n2 = 0;n2 < 3;n2++)
5336 {
5337 fGL tmp4 = tmp1[n2] / r1;
5338 fGL tmp5 = tmp4 * tmp3 / tmp2;
5339 dd[n2] += tmp5;
5340 }
5341 }
5342 }
5343
5344 return vdwsv;
5345 }
5346
5347 // f = sum[Q/r] = sum[Q * r^-1]
5348 // df/dr = -1 * sum[Q * r^-2]
5349
GetESP(fGL * pp,fGL * dd)5350 fGL eng1_sf::GetESP(fGL * pp, fGL * dd)
5351 {
5352 fGL espv = 0.0;
5353 if (dd != NULL) dd[0] = dd[1] = dd[2] = 0.0;
5354
5355 atom ** atmtab = GetSetup()->GetSFAtoms();
5356 for (i32s n1 = 0;n1 < GetSetup()->GetSFAtomCount() - num_solvent;n1++)
5357 {
5358 fGL tmp1[3]; fGL r2 = 0.0;
5359 const f64 * cdata = & crd[l2g_sf[n1] * 3];
5360 for (i32s n2 = 0;n2 < 3;n2++)
5361 {
5362 tmp1[n2] = pp[n2] - cdata[n2];
5363 r2 += tmp1[n2] * tmp1[n2];
5364 }
5365
5366 if (r2 == 0.0) return +1.0e+35; // numeric_limits<fGL>::max()?!?!?!
5367 fGL r1 = sqrt(r2);
5368
5369 // e = 2 + 76 * (( r / A ) ^ n) / (1 + ( r / A ) ^ n) + Z/r^9 = 2 + 76 * f / g
5370 // de/dr = 76 * (( g * Df - f * Dg ) / g ^ 2 )
5371
5372 f64 eps_n = myprm->epsilon1 + r2 * myprm->epsilon2; // assume constant!!!
5373 f64 eps_A = 1.25; // assume constant!!!
5374
5375 f64 t7a = r1 / eps_A;
5376 f64 t7b = pow(t7a, eps_n);
5377 f64 t7c = 1.0 + t7b;
5378 f64 t7d = 2.0 + 76.0 * (t7b / t7c);
5379
5380 f64 t7x = eps_n * pow(t7a, eps_n - 1.0) / eps_A;
5381 f64 t7y = 76.0 * ((t7c * t7x - t7b * t7x) / (t7c * t7c));
5382
5383 // do we have a correct constant here??? I think so, if we define
5384 // electrostatic potential as potential energy of a unit positive charge.
5385
5386 // fGL tmp2 = 4.1868 * 33.20716 * atmtab[n1]->charge / r1;
5387 fGL tmp2 = 4.1868 * 33.20716 * atmtab[n1]->charge / (t7d * r1); // screened!
5388 espv += tmp2;
5389
5390 if (dd != NULL) // sign ??? constant ???
5391 {
5392 for (i32s n2 = 0;n2 < 3;n2++)
5393 {
5394 fGL tmp3 = tmp1[n2] / r1;
5395 // fGL tmp4 = tmp3 * tmp2 / r1;
5396 fGL tmp4 = tmp3 * (-tmp2 * (1.0 / (t7d * r2) + t7y / (t7d * t7d * r1))); // screened!
5397 dd[n2] += tmp4;
5398 }
5399 }
5400 }
5401
5402 return espv;
5403 }
5404
5405 /*################################################################################################*/
5406
5407 // eof
5408