1 // SASAEVAL.CPP
2
3 // Copyright (C) 2000 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 "sasaeval.h"
22
23 #include "engine.h"
24
25 #include "local_i18n.h"
26 #include "notice.h"
27
28 #include <vector>
29 #include <algorithm>
30 #include <sstream>
31 using namespace std;
32
33 /*################################################################################################*/
34
35 // the surface area code apparently contains some bugs, since it sometimes
36 // crashes. another possibility is that the surface area math contains some
37 // bad cases (like arcs/segments with zero length/area ???) which should be
38 // avoided somehow. either way, the numerical and analytical gradients of
39 // surface area seem to match.
40
41 // LOWLIMIT is a trick to prevent zero-divisions in surface-area calculations.
42 // if a zero-division seems to happen, the values are just changed to LOWLIMIT
43 // in as early stage as possible, thus making minimum effects on results...
44
45 #define LOWLIMIT 0.0000001 // 0.0000001 seems to work quite well...
46
47 /*################################################################################################*/
48
sasaeval(engine * tmpe)49 sasaeval::sasaeval(engine * tmpe)
50 {
51 eng = tmpe;
52
53 natm_GLOB = eng->GetAtomCount();
54 natm_loc = NOT_DEFINED;
55
56 radius_GLOB = new f64[natm_GLOB];
57 index_GLOB_2_LOC = new i32u[natm_GLOB];
58
59 for (i32u i = 0;i < natm_GLOB;i++)
60 {
61 radius_GLOB[i] = -1.0;
62 index_GLOB_2_LOC[i] = NOT_DEFINED;
63 }
64
65 index_l2g = NULL;
66
67 radius1 = NULL;
68 radius2 = NULL;
69
70 dist1 = NULL;
71 dist2 = NULL;
72
73 nl = NULL;
74
75 sasa = NULL;
76 d_sasa = NULL;
77 }
78
~sasaeval(void)79 sasaeval::~sasaeval(void)
80 {
81 delete[] radius_GLOB;
82 radius_GLOB = NULL;
83
84 delete[] index_GLOB_2_LOC;
85 index_GLOB_2_LOC = NULL;
86
87 if (index_l2g != NULL)
88 {
89 delete[] index_l2g;
90 index_l2g = NULL;
91 }
92
93 if (radius1 != NULL)
94 {
95 delete[] radius1;
96 radius1 = NULL;
97 }
98
99 if (radius2 != NULL)
100 {
101 delete[] radius2;
102 radius2 = NULL;
103 }
104
105 if (dist1 != NULL)
106 {
107 delete[] dist1;
108 dist1 = NULL;
109 }
110
111 if (dist2 != NULL)
112 {
113 delete[] dist2;
114 dist2 = NULL;
115 }
116
117 if (nl != NULL)
118 {
119 delete[] nl;
120 nl = NULL;
121 }
122
123 if (sasa != NULL)
124 {
125 delete[] sasa;
126 sasa = NULL;
127 }
128
129 if (d_sasa != NULL)
130 {
131 delete[] d_sasa;
132 d_sasa = NULL;
133 }
134 }
135
RegisterAtom(i32u atmi_GLOB,double r)136 bool sasaeval::RegisterAtom(i32u atmi_GLOB, double r)
137 {
138 if (atmi_GLOB >= (i32u) natm_GLOB)
139 {
140 assertion_failed(__FILE__, __LINE__, "atmi_GLOB overflow.");
141 }
142
143 if (r < 0.001)
144 {
145 assertion_failed(__FILE__, __LINE__, "bad radius.");
146 }
147
148 if (radius_GLOB[atmi_GLOB] >= 0.0)
149 {
150 cout << _("WARNING : sasaeval::RegisterAtom() : atom ") << atmi_GLOB << _(" is already registered!") << endl;
151 return false;
152 }
153
154 radius_GLOB[atmi_GLOB] = r;
155 return true;
156 }
157
RegisterAtomsFinished(void)158 void sasaeval::RegisterAtomsFinished(void)
159 {
160 natm_loc = 0;
161
162 for (i32s i = 0;i < natm_GLOB;i++)
163 {
164 if (radius_GLOB[i] < 0.0)
165 {
166 index_GLOB_2_LOC[i] = NOT_DEFINED; // should be already...
167 }
168 else
169 {
170 index_GLOB_2_LOC[i] = natm_loc;
171 natm_loc++;
172 }
173 }
174
175 index_l2g = new i32u[natm_loc];
176
177 radius1 = new f64[natm_loc];
178 radius2 = new f64[natm_loc];
179
180 int localcounter = 0;
181 for (i32s i = 0;i < natm_GLOB;i++)
182 {
183 if (radius_GLOB[i] < 0.0) continue;
184
185 index_l2g[localcounter] = i;
186
187 const f64 r = radius_GLOB[i];
188
189 radius1[localcounter] = r;
190 radius2[localcounter] = r * r;
191
192 localcounter++;
193 }
194
195 dist1 = new i32s[natm_loc];
196 dist2 = new f64[natm_loc * (natm_loc - 1) / 2];
197
198 i32s n1 = 0; i32s n2 = 0;
199 while (n2 < natm_loc)
200 {
201 dist1[n2++] = n1;
202 n1 += natm_loc - n2;
203 }
204
205 nl = new cg_nbt3_nl[natm_loc];
206
207 for (i32s i = 0;i < natm_loc;i++)
208 {
209 nl[i].index = new i32s[SIZE_NLI];
210 }
211
212 sasa = new f64[natm_loc];
213 d_sasa = new f64[natm_loc * 3];
214 }
215
HandleNL(i32u atmiA_GLOB,i32u atmiB_GLOB,f64 dist)216 void sasaeval::HandleNL(i32u atmiA_GLOB, i32u atmiB_GLOB, f64 dist)
217 {
218 if (natm_loc < 0)
219 {
220 assertion_failed(__FILE__, __LINE__, "atom registration is not finished!");
221 }
222
223 bool bad_atoms = false;
224 if (atmiA_GLOB >= (i32u) natm_GLOB) bad_atoms = true;
225 if (atmiB_GLOB >= (i32u) natm_GLOB) bad_atoms = true;
226 if (atmiA_GLOB == atmiB_GLOB) bad_atoms = true;
227
228 if (bad_atoms)
229 {
230 ostringstream msg;
231 msg << "bad atoms " << atmiA_GLOB << " " << atmiB_GLOB << ends;
232 assertion_failed(__FILE__, __LINE__, msg.str().c_str());
233 }
234
235 const i32u atmi[2] = { index_GLOB_2_LOC[atmiA_GLOB], index_GLOB_2_LOC[atmiB_GLOB] };
236
237 const bool first = (atmi[0] > atmi[1]);
238 dist2[dist1[atmi[first]] + (atmi[!first] - atmi[first]) - 1] = dist;
239
240 if (dist < (radius1[atmi[0]] + radius1[atmi[1]]))
241 {
242 nl[atmi[0]].index[nl[atmi[0]].index_count++] = atmi[1];
243 if (nl[atmi[0]].index_count >= SIZE_NLI)
244 {
245 assertion_failed(__FILE__, __LINE__, "SASA NL index table overflow!");
246 }
247
248 nl[atmi[1]].index[nl[atmi[1]].index_count++] = atmi[0];
249 if (nl[atmi[1]].index_count >= SIZE_NLI)
250 {
251 assertion_failed(__FILE__, __LINE__, "SASA NL index table overflow!");
252 }
253 }
254 }
255
Evaluate(i32s p1)256 void sasaeval::Evaluate(i32s p1)
257 {
258 if (natm_loc < 0)
259 {
260 assertion_failed(__FILE__, __LINE__, "atom registration is not finished!");
261 }
262
263 for (i32s n1 = 0;n1 < natm_loc;n1++)
264 {
265 nl[n1].index_count = 0;
266
267 sasa[n1] = 0.0;
268
269 if (p1 > 0)
270 {
271 d_sasa[n1 * 3 + 0] = 0.0;
272 d_sasa[n1 * 3 + 1] = 0.0;
273 d_sasa[n1 * 3 + 2] = 0.0;
274 }
275 }
276
277 cg_nbt3_nl * nlist = nl;
278
279 for (i32s n1 = 0;n1 < natm_loc;n1++)
280 {
281 cg_nbt3_nd ndt[SIZE_NLI];
282 for (i32s n2 = 0;n2 < nlist[n1].index_count;n2++)
283 {
284 ndt[n2].index = nlist[n1].index[n2];
285 i32s atmi[2] = { n1, ndt[n2].index }; bool first = (atmi[0] > atmi[1]);
286 ndt[n2].distance = dist2[dist1[atmi[first]] + (atmi[!first] - atmi[first]) - 1];
287 }
288
289 sort(ndt, ndt + nlist[n1].index_count);
290 i32s n_count = 0; i32s nt[SIZE_NT];
291
292 // neighbor-list reduction... THIS WON'T WORK IF ANY OF THE BT1/NBT1-TERMS ARE LEFT OUT!!!
293 // neighbor-list reduction... THIS WON'T WORK IF ANY OF THE BT1/NBT1-TERMS ARE LEFT OUT!!!
294 // neighbor-list reduction... THIS WON'T WORK IF ANY OF THE BT1/NBT1-TERMS ARE LEFT OUT!!!
295
296 // test this against a slow-but-simple implementation?!?!?!
297 // seems to be OK because different layers give similar results...
298
299 for (i32s n2 = 0;n2 < nlist[n1].index_count;n2++)
300 {
301 i32s ind1 = ndt[n2].index;
302 f64 dij = ndt[n2].distance;
303
304 bool flag = true;
305
306 for (i32s n3 = n2 + 1;n3 < nlist[n1].index_count;n3++)
307 {
308 i32s ind2 = ndt[n3].index;
309
310 i32s atmi[2] = { ind1, ind2 }; bool first = (atmi[0] > atmi[1]);
311 f64 djk = dist2[dist1[atmi[first]] + (atmi[!first] - atmi[first]) - 1];
312
313 if (djk > dij) continue;
314
315 f64 dij2 = dij * dij; f64 djk2 = djk * djk;
316 f64 dik = ndt[n3].distance; f64 dik2 = dik * dik;
317
318 // here dij and dik both represent distances which should never be
319 // very close to zero (if LJ-terms work as they should) -> no checking
320
321 f64 ca = (radius2[n1] + dij2 - radius2[ind1]) / (2.0 * radius1[n1] * dij);
322 f64 cb = (radius2[n1] + dik2 - radius2[ind2]) / (2.0 * radius1[n1] * dik);
323 f64 cg = (dij2 + dik2 - djk2) / (2.0 * dij * dik);
324
325 f64 sa2 = 1.0 - ca * ca;
326 f64 sg2 = 1.0 - cg * cg;
327
328 f64 dc = sa2 * sg2;
329 if (dc < 0.0) dc = 0.0; // domain check...
330
331 if (cb < ca * cg - sqrt(dc))
332 {
333 flag = false;
334 break;
335 }
336 }
337
338 if (flag)
339 {
340 nt[n_count++] = ind1;
341 if (n_count >= SIZE_NT)
342 {
343 assertion_failed(__FILE__, __LINE__, "NT overflow!");
344 }
345 }
346 }
347
348 i32s coi_count = 0; cg_nbt3_coi coit[SIZE_COI];
349
350 // next we will create the coi-table...
351 // next we will create the coi-table...
352 // next we will create the coi-table...
353
354 for (i32s n2 = 0;n2 < n_count;n2++)
355 {
356 coit[coi_count].index = nt[n2]; coit[coi_count].flag = false;
357
358 f64 t1a[3]; f64 t1b = 0.0;
359 for (i32s n3 = 0;n3 < 3;n3++)
360 {
361 i32s tmpi;
362
363 tmpi = index_l2g[n1];
364 f64 t9a = eng->crd[tmpi * 3 + n3];
365
366 tmpi = index_l2g[coit[coi_count].index];
367 f64 t9b = eng->crd[tmpi * 3 + n3];
368
369 t1a[n3] = t9b - t9a;
370 t1b += t1a[n3] * t1a[n3];
371 }
372
373 f64 t1c = sqrt(t1b);
374 coit[coi_count].dist = t1c;
375
376 // also t1c is a distance which should never be close to zero -> no checking
377
378 f64 t2a[3];
379 for (i32s n3 = 0;n3 < 3;n3++)
380 {
381 t2a[n3] = t1a[n3] / t1c;
382 coit[coi_count].dv[n3] = t2a[n3];
383 }
384
385 coit[coi_count].g = (t1b + radius2[n1] - radius2[coit[coi_count].index]) / (2.0 * t1c);
386 coit[coi_count].ct = coit[coi_count].g / radius1[n1];
387
388 if (p1 > 0)
389 {
390 for (i32s n3 = 0;n3 < 3;n3++)
391 {
392 for (i32s n4 = 0;n4 < 3;n4++)
393 {
394 f64 t9a = t2a[n3] * t2a[n4]; f64 t9b;
395 if (n3 != n4) t9b = -t9a; else t9b = 1.0 - t9a;
396 coit[coi_count].ddv[n3][n4] = t9b / t1c;
397 }
398 }
399
400 f64 t3a = (t1c - coit[coi_count].g) / t1c;
401 coit[coi_count].dg[0] = t3a * coit[coi_count].dv[0];
402 coit[coi_count].dg[1] = t3a * coit[coi_count].dv[1];
403 coit[coi_count].dg[2] = t3a * coit[coi_count].dv[2];
404
405 coit[coi_count].dct[0] = coit[coi_count].dg[0] / radius1[n1];
406 coit[coi_count].dct[1] = coit[coi_count].dg[1] / radius1[n1];
407 coit[coi_count].dct[2] = coit[coi_count].dg[2] / radius1[n1];
408 }
409
410 coit[coi_count++].ipd_count = 0;
411 if (coi_count >= SIZE_COI)
412 {
413 assertion_failed(__FILE__, __LINE__, "COI overflow!");
414 }
415 }
416
417 i32s ips_total_count = 0;
418 i32s ips_count = 0; cg_nbt3_ips ipst[SIZE_IPS];
419
420 // next we will create the ips-table...
421 // next we will create the ips-table...
422 // next we will create the ips-table...
423
424 for (i32s n2 = 0;n2 < coi_count - 1;n2++)
425 {
426 for (i32s n3 = n2 + 1;n3 < coi_count;n3++)
427 {
428 f64 t1a[3];
429 t1a[0] = coit[n2].dv[0] * coit[n3].dv[0];
430 t1a[1] = coit[n2].dv[1] * coit[n3].dv[1];
431 t1a[2] = coit[n2].dv[2] * coit[n3].dv[2];
432
433 f64 t1b = t1a[0] + t1a[1] + t1a[2]; // cos phi
434
435 if (t1b < -1.0) t1b = -1.0; // domain check...
436 if (t1b > +1.0) t1b = +1.0; // domain check...
437
438 f64 t1c = 1.0 - t1b * t1b; // sin^2 phi
439 if (t1c < LOWLIMIT) t1c = LOWLIMIT;
440
441 f64 t2a = (coit[n2].g - coit[n3].g * t1b) / t1c; // tau_kj
442 f64 t2b = (coit[n3].g - coit[n2].g * t1b) / t1c; // tau_jk
443
444 f64 t2c = radius2[n1] - coit[n2].g * t2a - coit[n3].g * t2b; // gamma^2
445 if (t2c < LOWLIMIT) continue; // these will not intercept...
446
447 ips_total_count++;
448 coit[n2].flag = true;
449 coit[n3].flag = true;
450
451 f64 t3a[3]; // eta
452 t3a[0] = coit[n2].dv[0] * t2a + coit[n3].dv[0] * t2b;
453 t3a[1] = coit[n2].dv[1] * t2a + coit[n3].dv[1] * t2b;
454 t3a[2] = coit[n2].dv[2] * t2a + coit[n3].dv[2] * t2b;
455
456 f64 t1d = sqrt(t1c); // sin phi
457
458 f64 t3b[3]; // omega
459 t3b[0] = (coit[n2].dv[1] * coit[n3].dv[2] - coit[n2].dv[2] * coit[n3].dv[1]) / t1d;
460 t3b[1] = (coit[n2].dv[2] * coit[n3].dv[0] - coit[n2].dv[0] * coit[n3].dv[2]) / t1d;
461 t3b[2] = (coit[n2].dv[0] * coit[n3].dv[1] - coit[n2].dv[1] * coit[n3].dv[0]) / t1d;
462
463 f64 t2d = sqrt(t2c); // gamma
464
465 for (i32s n4 = 0;n4 < 3;n4++)
466 {
467 f64 t9a = t3b[n4] * t2d;
468 ipst[ips_count].ipv[0][n4] = t3a[n4] - t9a;
469 ipst[ips_count].ipv[1][n4] = t3a[n4] + t9a;
470 }
471
472 // skip those intersection points that fall inside any other sphere...
473 // SKIP ALSO IF EQUAL DISTANCE??? i.e. compare using "<" or "<=" ???
474
475 bool skip_both = false;
476 bool skip[2] = { false, false };
477 for (i32s n4 = 0;n4 < n_count;n4++)
478 {
479 i32s n5 = nt[n4];
480 if (n5 == coit[n2].index || n5 == coit[n3].index) continue;
481
482 f64 t9a[3];
483 t9a[0] = (eng->crd[index_l2g[n1] * 3 + 0] + ipst[ips_count].ipv[0][0]) - eng->crd[index_l2g[n5] * 3 + 0];
484 t9a[1] = (eng->crd[index_l2g[n1] * 3 + 1] + ipst[ips_count].ipv[0][1]) - eng->crd[index_l2g[n5] * 3 + 1];
485 t9a[2] = (eng->crd[index_l2g[n1] * 3 + 2] + ipst[ips_count].ipv[0][2]) - eng->crd[index_l2g[n5] * 3 + 2];
486 f64 t9b = t9a[0] * t9a[0] + t9a[1] * t9a[1] + t9a[2] * t9a[2];
487 if (t9b < radius2[n5]) skip[0] = true;
488
489 f64 t9c[3];
490 t9c[0] = (eng->crd[index_l2g[n1] * 3 + 0] + ipst[ips_count].ipv[1][0]) - eng->crd[index_l2g[n5] * 3 + 0];
491 t9c[1] = (eng->crd[index_l2g[n1] * 3 + 1] + ipst[ips_count].ipv[1][1]) - eng->crd[index_l2g[n5] * 3 + 1];
492 t9c[2] = (eng->crd[index_l2g[n1] * 3 + 2] + ipst[ips_count].ipv[1][2]) - eng->crd[index_l2g[n5] * 3 + 2];
493 f64 t9d = t9c[0] * t9c[0] + t9c[1] * t9c[1] + t9c[2] * t9c[2];
494 if (t9d < radius2[n5]) skip[1] = true;
495
496 skip_both = (skip[0] && skip[1]);
497 if (skip_both) break;
498 }
499
500 if (skip_both) continue; // overwrite this one...
501
502 ipst[ips_count].coi[0] = n2;
503 ipst[ips_count].coi[1] = n3;
504
505 if (!skip[0])
506 {
507 coit[n2].AddIPD(ipst[ips_count].ipv[0], ips_count);
508 coit[n3].AddIPD(ipst[ips_count].ipv[0], ips_count | ORDER_FLAG);
509 }
510
511 if (!skip[1])
512 {
513 coit[n2].AddIPD(ipst[ips_count].ipv[1], ips_count | INDEX_FLAG | ORDER_FLAG);
514 coit[n3].AddIPD(ipst[ips_count].ipv[1], ips_count | INDEX_FLAG);
515 }
516
517 if (p1 > 0)
518 {
519 f64 t1f[3]; // d(cos phi) / dXk
520 t1f[0] = (coit[n3].dv[0] - t1b * coit[n2].dv[0]) / coit[n2].dist;
521 t1f[1] = (coit[n3].dv[1] - t1b * coit[n2].dv[1]) / coit[n2].dist;
522 t1f[2] = (coit[n3].dv[2] - t1b * coit[n2].dv[2]) / coit[n2].dist;
523
524 f64 t1g[3]; // d(cos phi) / dXj
525 t1g[0] = (coit[n2].dv[0] - t1b * coit[n3].dv[0]) / coit[n3].dist;
526 t1g[1] = (coit[n2].dv[1] - t1b * coit[n3].dv[1]) / coit[n3].dist;
527 t1g[2] = (coit[n2].dv[2] - t1b * coit[n3].dv[2]) / coit[n3].dist;
528
529 f64 t2e[3]; // d(tau_kj) / dXk
530 t2e[0] = (t1f[0] * (2.0 * t2a * t1b - coit[n3].g) + coit[n2].dg[0]) / t1c;
531 t2e[1] = (t1f[1] * (2.0 * t2a * t1b - coit[n3].g) + coit[n2].dg[1]) / t1c;
532 t2e[2] = (t1f[2] * (2.0 * t2a * t1b - coit[n3].g) + coit[n2].dg[2]) / t1c;
533
534 f64 t2f[3]; // d(tau_kj) / dXj
535 t2f[0] = (t1g[0] * (2.0 * t2a * t1b - coit[n3].g) - t1b * coit[n3].dg[0]) / t1c;
536 t2f[1] = (t1g[1] * (2.0 * t2a * t1b - coit[n3].g) - t1b * coit[n3].dg[1]) / t1c;
537 t2f[2] = (t1g[2] * (2.0 * t2a * t1b - coit[n3].g) - t1b * coit[n3].dg[2]) / t1c;
538
539 f64 t2g[3]; // d(tau_jk) / dXk
540 t2g[0] = (t1f[0] * (2.0 * t2b * t1b - coit[n2].g) - t1b * coit[n2].dg[0]) / t1c;
541 t2g[1] = (t1f[1] * (2.0 * t2b * t1b - coit[n2].g) - t1b * coit[n2].dg[1]) / t1c;
542 t2g[2] = (t1f[2] * (2.0 * t2b * t1b - coit[n2].g) - t1b * coit[n2].dg[2]) / t1c;
543
544 f64 t2h[3]; // d(tau_jk) / dXj
545 t2h[0] = (t1g[0] * (2.0 * t2b * t1b - coit[n2].g) + coit[n3].dg[0]) / t1c;
546 t2h[1] = (t1g[1] * (2.0 * t2b * t1b - coit[n2].g) + coit[n3].dg[1]) / t1c;
547 t2h[2] = (t1g[2] * (2.0 * t2b * t1b - coit[n2].g) + coit[n3].dg[2]) / t1c;
548
549 f64 t3c[3][3]; // d(eta) / dXk
550 f64 t3d[3][3]; // d(eta) / dXj
551
552 for (i32s n4 = 0;n4 < 3;n4++)
553 {
554 for (i32s n5 = 0;n5 < 3;n5++)
555 {
556 f64 t9a = coit[n2].dv[n5]; f64 t9b = coit[n3].dv[n5];
557 t3c[n4][n5] = t9a * t2e[n4] + t9b * t2g[n4] + t2a * coit[n2].ddv[n4][n5];
558 t3d[n4][n5] = t9a * t2f[n4] + t9b * t2h[n4] + t2b * coit[n3].ddv[n4][n5];
559 }
560 }
561
562 f64 t3e[3][3]; // d(omega) / dXk
563 f64 t3f[3][3]; // d(omega) / dXj
564
565 for (i32s n4 = 0;n4 < 3;n4++)
566 {
567 for (i32s n5 = 0;n5 < 3;n5++)
568 {
569 t3e[n4][n5] = t1b * t3b[n5] * t1f[n4] / t1c;
570 t3f[n4][n5] = t1b * t3b[n5] * t1g[n4] / t1c;
571 }
572
573 t3e[n4][0] += (coit[n2].ddv[n4][1] * coit[n3].dv[2] - coit[n2].ddv[n4][2] * coit[n3].dv[1]) / t1d;
574 t3e[n4][1] += (coit[n2].ddv[n4][2] * coit[n3].dv[0] - coit[n2].ddv[n4][0] * coit[n3].dv[2]) / t1d;
575 t3e[n4][2] += (coit[n2].ddv[n4][0] * coit[n3].dv[1] - coit[n2].ddv[n4][1] * coit[n3].dv[0]) / t1d;
576
577 t3f[n4][0] += (coit[n2].dv[1] * coit[n3].ddv[n4][2] - coit[n2].dv[2] * coit[n3].ddv[n4][1]) / t1d;
578 t3f[n4][1] += (coit[n2].dv[2] * coit[n3].ddv[n4][0] - coit[n2].dv[0] * coit[n3].ddv[n4][2]) / t1d;
579 t3f[n4][2] += (coit[n2].dv[0] * coit[n3].ddv[n4][1] - coit[n2].dv[1] * coit[n3].ddv[n4][0]) / t1d;
580
581 }
582
583 f64 t2i[3]; // d(gamma) / dXk
584 t2i[0] = -(coit[n2].g * t2e[0] + t2a * coit[n2].dg[0] + coit[n3].g * t2g[0]) / (2.0 * t2d);
585 t2i[1] = -(coit[n2].g * t2e[1] + t2a * coit[n2].dg[1] + coit[n3].g * t2g[1]) / (2.0 * t2d);
586 t2i[2] = -(coit[n2].g * t2e[2] + t2a * coit[n2].dg[2] + coit[n3].g * t2g[2]) / (2.0 * t2d);
587
588 f64 t2j[3]; // d(gamma) / dXj
589 t2j[0] = -(coit[n2].g * t2f[0] + coit[n3].g * t2h[0] + t2b * coit[n3].dg[0]) / (2.0 * t2d);
590 t2j[1] = -(coit[n2].g * t2f[1] + coit[n3].g * t2h[1] + t2b * coit[n3].dg[1]) / (2.0 * t2d);
591 t2j[2] = -(coit[n2].g * t2f[2] + coit[n3].g * t2h[2] + t2b * coit[n3].dg[2]) / (2.0 * t2d);
592
593 // the final result is derivatives for points dipv[2][2][3][3].
594 // indexes are as follows: [point][atom][variable][xyz].
595
596 for (i32s n4 = 0;n4 < 3;n4++)
597 {
598 for (i32s n5 = 0;n5 < 3;n5++)
599 {
600 ipst[ips_count].dipv[0][0][n4][n5] = t3c[n4][n5];
601 ipst[ips_count].dipv[0][1][n4][n5] = t3d[n4][n5];
602 ipst[ips_count].dipv[1][0][n4][n5] = t3c[n4][n5];
603 ipst[ips_count].dipv[1][1][n4][n5] = t3d[n4][n5];
604 }
605
606 for (i32s n5 = 0;n5 < 3;n5++)
607 {
608 f64 t9a = t3b[n5] * t2i[n4] + t2d * t3e[n4][n5];
609 f64 t9b = t3b[n5] * t2j[n4] + t2d * t3f[n4][n5];
610
611 ipst[ips_count].dipv[0][0][n4][n5] -= t9a;
612 ipst[ips_count].dipv[0][1][n4][n5] -= t9b;
613 ipst[ips_count].dipv[1][0][n4][n5] += t9a;
614 ipst[ips_count].dipv[1][1][n4][n5] += t9b;
615 }
616 }
617 }
618
619 ips_count++;
620 if (ips_count >= SIZE_IPS)
621 {
622 assertion_failed(__FILE__, __LINE__, "IPS overflow!");
623 }
624 }
625 }
626
627 i32s arc_count = 0; cg_nbt3_arc arct[SIZE_ARC];
628
629 // next we will create the arc-table...
630 // next we will create the arc-table...
631 // next we will create the arc-table...
632
633 for (i32s n2 = 0;n2 < coi_count;n2++)
634 {
635 f64 t1z = radius2[n1] - coit[n2].g * coit[n2].g;
636 if (t1z < 0.0) t1z = 0.0; // domain check...
637
638 f64 t1a = sqrt(t1z);
639 if (t1a < LOWLIMIT) t1a = LOWLIMIT;
640
641 sort(coit[n2].ipdt, coit[n2].ipdt + coit[n2].ipd_count);
642
643 for (i32s n3 = 0;n3 < coit[n2].ipd_count;n3++)
644 {
645 if (coit[n2].ipdt[n3].ipdata & ORDER_FLAG) continue;
646 i32s n4 = n3 + 1; if (n4 == coit[n2].ipd_count) n4 = 0;
647 if (!(coit[n2].ipdt[n4].ipdata & ORDER_FLAG)) continue;
648
649 arct[arc_count].coi = n2; arct[arc_count].flag = false;
650
651 arct[arc_count].ipdata[0] = (coit[n2].ipdt[n3].ipdata & ~ORDER_FLAG);
652 arct[arc_count].ipdata[1] = (coit[n2].ipdt[n4].ipdata & ~ORDER_FLAG);
653
654 i32s i1a = (arct[arc_count].ipdata[0] & FLAG_MASK);
655 bool i1b = (arct[arc_count].ipdata[0] & INDEX_FLAG ? 1 : 0);
656
657 i32s i2a = (arct[arc_count].ipdata[1] & FLAG_MASK);
658 bool i2b = (arct[arc_count].ipdata[1] & INDEX_FLAG ? 1 : 0);
659
660 arct[arc_count].index[0][0] = coit[ipst[i1a].coi[i1b]].index;
661 arct[arc_count].index[0][1] = coit[ipst[i1a].coi[!i1b]].index;
662
663 arct[arc_count].index[1][0] = coit[ipst[i2a].coi[!i2b]].index;
664 arct[arc_count].index[1][1] = coit[ipst[i2a].coi[i2b]].index;
665
666 // let's compute the tangent vectors...
667
668 f64 * ref1 = ipst[i1a].ipv[i1b];
669 arct[arc_count].tv[0][0] = (ref1[1] * coit[n2].dv[2] - ref1[2] * coit[n2].dv[1]) / t1a;
670 arct[arc_count].tv[0][1] = (ref1[2] * coit[n2].dv[0] - ref1[0] * coit[n2].dv[2]) / t1a;
671 arct[arc_count].tv[0][2] = (ref1[0] * coit[n2].dv[1] - ref1[1] * coit[n2].dv[0]) / t1a;
672
673 f64 * ref2 = ipst[i2a].ipv[i2b];
674 arct[arc_count].tv[1][0] = (ref2[1] * coit[n2].dv[2] - ref2[2] * coit[n2].dv[1]) / t1a;
675 arct[arc_count].tv[1][1] = (ref2[2] * coit[n2].dv[0] - ref2[0] * coit[n2].dv[2]) / t1a;
676 arct[arc_count].tv[1][2] = (ref2[0] * coit[n2].dv[1] - ref2[1] * coit[n2].dv[0]) / t1a;
677
678 if (p1 > 0)
679 {
680 for (i32s n4 = 0;n4 < 3;n4++)
681 {
682 f64 t9a = coit[n2].g * coit[n2].dg[n4] / t1a;
683 for (i32s n5 = 0;n5 < 3;n5++)
684 {
685 arct[arc_count].dtv[0][0][n4][n5] = t9a * arct[arc_count].tv[0][n5];
686 arct[arc_count].dtv[1][0][n4][n5] = t9a * arct[arc_count].tv[1][n5];
687 }
688
689 f64 * ref1a = ipst[i1a].dipv[i1b][i1b][n4]; // d(P1) / dXk
690 arct[arc_count].dtv[0][0][n4][0] += ref1a[1] * coit[n2].dv[2] - ref1a[2] * coit[n2].dv[1];
691 arct[arc_count].dtv[0][0][n4][1] += ref1a[2] * coit[n2].dv[0] - ref1a[0] * coit[n2].dv[2];
692 arct[arc_count].dtv[0][0][n4][2] += ref1a[0] * coit[n2].dv[1] - ref1a[1] * coit[n2].dv[0];
693
694 f64 * ref1b = ipst[i2a].dipv[i2b][!i2b][n4]; // d(P2) / dXk
695 arct[arc_count].dtv[1][0][n4][0] += ref1b[1] * coit[n2].dv[2] - ref1b[2] * coit[n2].dv[1];
696 arct[arc_count].dtv[1][0][n4][1] += ref1b[2] * coit[n2].dv[0] - ref1b[0] * coit[n2].dv[2];
697 arct[arc_count].dtv[1][0][n4][2] += ref1b[0] * coit[n2].dv[1] - ref1b[1] * coit[n2].dv[0];
698
699 f64 * ref2a = ipst[i1a].ipv[i1b];
700 arct[arc_count].dtv[0][0][n4][0] += ref2a[1] * coit[n2].ddv[n4][2] - ref2a[2] * coit[n2].ddv[n4][1];
701 arct[arc_count].dtv[0][0][n4][1] += ref2a[2] * coit[n2].ddv[n4][0] - ref2a[0] * coit[n2].ddv[n4][2];
702 arct[arc_count].dtv[0][0][n4][2] += ref2a[0] * coit[n2].ddv[n4][1] - ref2a[1] * coit[n2].ddv[n4][0];
703
704 f64 * ref2b = ipst[i2a].ipv[i2b];
705 arct[arc_count].dtv[1][0][n4][0] += ref2b[1] * coit[n2].ddv[n4][2] - ref2b[2] * coit[n2].ddv[n4][1];
706 arct[arc_count].dtv[1][0][n4][1] += ref2b[2] * coit[n2].ddv[n4][0] - ref2b[0] * coit[n2].ddv[n4][2];
707 arct[arc_count].dtv[1][0][n4][2] += ref2b[0] * coit[n2].ddv[n4][1] - ref2b[1] * coit[n2].ddv[n4][0];
708
709 for (i32s n5 = 0;n5 < 3;n5++)
710 {
711 arct[arc_count].dtv[0][0][n4][n5] /= t1a;
712 arct[arc_count].dtv[1][0][n4][n5] /= t1a;
713 }
714
715 f64 * ref3a = ipst[i1a].dipv[i1b][!i1b][n4]; // d(P1) / dXj
716 arct[arc_count].dtv[0][1][n4][0] = (ref3a[1] * coit[n2].dv[2] - ref3a[2] * coit[n2].dv[1]) / t1a;
717 arct[arc_count].dtv[0][1][n4][1] = (ref3a[2] * coit[n2].dv[0] - ref3a[0] * coit[n2].dv[2]) / t1a;
718 arct[arc_count].dtv[0][1][n4][2] = (ref3a[0] * coit[n2].dv[1] - ref3a[1] * coit[n2].dv[0]) / t1a;
719
720 f64 * ref3b = ipst[i2a].dipv[i2b][i2b][n4]; // d(P2) / dXj
721 arct[arc_count].dtv[1][1][n4][0] = (ref3b[1] * coit[n2].dv[2] - ref3b[2] * coit[n2].dv[1]) / t1a;
722 arct[arc_count].dtv[1][1][n4][1] = (ref3b[2] * coit[n2].dv[0] - ref3b[0] * coit[n2].dv[2]) / t1a;
723 arct[arc_count].dtv[1][1][n4][2] = (ref3b[0] * coit[n2].dv[1] - ref3b[1] * coit[n2].dv[0]) / t1a;
724 }
725 }
726
727 arc_count++;
728 if (arc_count >= SIZE_ARC)
729 {
730 assertion_failed(__FILE__, __LINE__, "ARC overflow!");
731 }
732 }
733 }
734
735 // all cases will pass through this point!!!
736 // all cases will pass through this point!!!
737 // all cases will pass through this point!!!
738
739 f64 area;
740 if (!arc_count)
741 {
742 if (ips_total_count)
743 {
744 // save the solv-exp value here if needed!!!
745 // save the solv-exp value here if needed!!!
746 // save the solv-exp value here if needed!!!
747 continue; // fully buried...
748 }
749 else area = 4.0 * M_PI;
750 }
751 else
752 {
753 area = 0.0;
754 i32s arc_counter = 0;
755
756 do
757 {
758 i32s prev; i32s curr = 0;
759 while (arct[curr].flag)
760 {
761 curr++;
762 if (curr == arc_count)
763 {
764 cout << "area_panic: can't find the first arc!!!" << endl;
765 goto area_panic;
766 }
767 }
768
769 i32s first = curr;
770
771 f64 sum1 = 0.0;
772 f64 sum2 = 0.0;
773
774 while (true)
775 {
776 i32s coi = arct[curr].coi;
777
778 f64 t1a[3];
779 t1a[0] = arct[curr].tv[1][1] * arct[curr].tv[0][2] - arct[curr].tv[1][2] * arct[curr].tv[0][1];
780 t1a[1] = arct[curr].tv[1][2] * arct[curr].tv[0][0] - arct[curr].tv[1][0] * arct[curr].tv[0][2];
781 t1a[2] = arct[curr].tv[1][0] * arct[curr].tv[0][1] - arct[curr].tv[1][1] * arct[curr].tv[0][0];
782
783 f64 t1b[3];
784 t1b[0] = coit[coi].dv[0] * t1a[0];
785 t1b[1] = coit[coi].dv[1] * t1a[1];
786 t1b[2] = coit[coi].dv[2] * t1a[2];
787
788 f64 t1c = (t1b[0] + t1b[1] + t1b[2] < 0.0 ? -1.0 : +1.0);
789
790 f64 t2a[3];
791 t2a[0] = arct[curr].tv[0][0] * arct[curr].tv[1][0];
792 t2a[1] = arct[curr].tv[0][1] * arct[curr].tv[1][1];
793 t2a[2] = arct[curr].tv[0][2] * arct[curr].tv[1][2];
794
795 f64 t2b = t2a[0] + t2a[1] + t2a[2];
796
797 if (t2b < -1.0) t2b = -1.0; // domain check...
798 if (t2b > +1.0) t2b = +1.0; // domain check...
799
800 f64 t2c = (1.0 - t1c) * M_PI + t1c * acos(t2b);
801 sum1 += t2c * coit[coi].ct;
802
803 if (p1 > 0)
804 {
805 f64 t2x = fabs(sin(t2c));
806 if (t2x < LOWLIMIT) t2x = LOWLIMIT;
807
808 f64 t2y = -coit[coi].ct * t1c / t2x;
809
810 // 1st are same points and 2nd are different ones...
811 // 1st are same points and 2nd are different ones...
812 // 1st are same points and 2nd are different ones...
813
814 for (i32s n2 = 0;n2 < 3;n2++)
815 {
816 f64 t3a[3];
817 t3a[0] = arct[curr].dtv[0][0][n2][0] * arct[curr].tv[1][0];
818 t3a[1] = arct[curr].dtv[0][0][n2][1] * arct[curr].tv[1][1];
819 t3a[2] = arct[curr].dtv[0][0][n2][2] * arct[curr].tv[1][2];
820 f64 t3b = t3a[0] + t3a[1] + t3a[2];
821
822 f64 t3c[3];
823 t3c[0] = arct[curr].tv[0][0] * arct[curr].dtv[1][0][n2][0];
824 t3c[1] = arct[curr].tv[0][1] * arct[curr].dtv[1][0][n2][1];
825 t3c[2] = arct[curr].tv[0][2] * arct[curr].dtv[1][0][n2][2];
826 f64 t3d = t3c[0] + t3c[1] + t3c[2];
827
828 f64 t4a[3];
829 t4a[0] = arct[curr].dtv[0][1][n2][0] * arct[curr].tv[1][0];
830 t4a[1] = arct[curr].dtv[0][1][n2][1] * arct[curr].tv[1][1];
831 t4a[2] = arct[curr].dtv[0][1][n2][2] * arct[curr].tv[1][2];
832 f64 t4b = t4a[0] + t4a[1] + t4a[2];
833
834 f64 t4c[3];
835 t4c[0] = arct[curr].tv[0][0] * arct[curr].dtv[1][1][n2][0];
836 t4c[1] = arct[curr].tv[0][1] * arct[curr].dtv[1][1][n2][1];
837 t4c[2] = arct[curr].tv[0][2] * arct[curr].dtv[1][1][n2][2];
838 f64 t4d = t4c[0] + t4c[1] + t4c[2];
839
840 f64 t3e = t2y * (t3b + t3d) + t2c * coit[coi].dct[n2];
841 f64 t5a = t2y * t4b; f64 t5b = t2y * t4d;
842
843 d_sasa[arct[curr].index[0][0] * 3 + n2] += t3e;
844 d_sasa[arct[curr].index[0][1] * 3 + n2] += t5a;
845 d_sasa[arct[curr].index[1][1] * 3 + n2] += t5b;
846 d_sasa[n1 * 3 + n2] -= t3e + t5a + t5b;
847 }
848 }
849
850 prev = curr; curr = 0;
851 i32u ipd = arct[prev].ipdata[1];
852 while (true)
853 {
854 if (arct[curr].ipdata[0] != ipd) curr++;
855 else break;
856
857 if (curr == arc_count)
858 {
859 cout << "area_panic: incomplete set of arcs!!!" << endl;
860 goto area_panic;
861 }
862 }
863
864 arc_counter++;
865 arct[curr].flag = true;
866
867 f64 t2d[3];
868 t2d[0] = arct[prev].tv[1][0] * arct[curr].tv[0][0];
869 t2d[1] = arct[prev].tv[1][1] * arct[curr].tv[0][1];
870 t2d[2] = arct[prev].tv[1][2] * arct[curr].tv[0][2];
871
872 f64 t2e = t2d[0] + t2d[1] + t2d[2];
873
874 if (t2e < -1.0) t2e = -1.0; // domain check...
875 if (t2e > +1.0) t2e = +1.0; // domain check...
876
877 f64 t2f = -acos(t2e); sum2 += t2f;
878
879 if (p1 > 0)
880 {
881 f64 t2x = fabs(sin(t2f));
882 if (t2x < LOWLIMIT) t2x = LOWLIMIT;
883
884 f64 t2y = 1.0 / t2x;
885
886 // prev_k = curr_j and prev_j = curr_k !!!
887 // prev_k = curr_j and prev_j = curr_k !!!
888 // prev_k = curr_j and prev_j = curr_k !!!
889
890 for (i32s n2 = 0;n2 < 3;n2++)
891 {
892 f64 t3a[3];
893 t3a[0] = arct[prev].dtv[1][0][n2][0] * arct[curr].tv[0][0];
894 t3a[1] = arct[prev].dtv[1][0][n2][1] * arct[curr].tv[0][1];
895 t3a[2] = arct[prev].dtv[1][0][n2][2] * arct[curr].tv[0][2];
896 f64 t3b = t3a[0] + t3a[1] + t3a[2];
897
898 f64 t3c[3];
899 t3c[0] = arct[prev].tv[1][0] * arct[curr].dtv[0][1][n2][0];
900 t3c[1] = arct[prev].tv[1][1] * arct[curr].dtv[0][1][n2][1];
901 t3c[2] = arct[prev].tv[1][2] * arct[curr].dtv[0][1][n2][2];
902 f64 t3d = t3c[0] + t3c[1] + t3c[2];
903
904 f64 t4a[3];
905 t4a[0] = arct[prev].dtv[1][1][n2][0] * arct[curr].tv[0][0];
906 t4a[1] = arct[prev].dtv[1][1][n2][1] * arct[curr].tv[0][1];
907 t4a[2] = arct[prev].dtv[1][1][n2][2] * arct[curr].tv[0][2];
908 f64 t4b = t4a[0] + t4a[1] + t4a[2];
909
910 f64 t4c[3];
911 t4c[0] = arct[prev].tv[1][0] * arct[curr].dtv[0][0][n2][0];
912 t4c[1] = arct[prev].tv[1][1] * arct[curr].dtv[0][0][n2][1];
913 t4c[2] = arct[prev].tv[1][2] * arct[curr].dtv[0][0][n2][2];
914 f64 t4d = t4c[0] + t4c[1] + t4c[2];
915
916 f64 t3e = t2y * (t3b + t3d);
917 f64 t4e = t2y * (t4b + t4d);
918
919 d_sasa[arct[prev].index[1][0] * 3 + n2] += t3e;
920 d_sasa[arct[prev].index[1][1] * 3 + n2] += t4e;
921 d_sasa[n1 * 3 + n2] -= t3e + t4e;
922 }
923 }
924
925 if (curr == first) break;
926 }
927
928 area += 2.0 * M_PI + sum1 + sum2;
929 } while (arc_counter < arc_count);
930
931 // when we have some problems somewhere above (for example, if we have
932 // an incomplete set of arcs or no arcs at all; these things are possible
933 // in rare special cases; for example we might have to reject some arcs
934 // if they contained some singular intermediate values) we will truncate
935 // the sum and jump right here.
936
937 // in this case we will calculate incorrect value for the area, but the
938 // good news is that the value and the gradient will still be consistent.
939
940 // since these cases are very rare, this probably won't make big problems
941 // in any applications...
942
943 area_panic: // we will jump here in all problematic cases...
944
945 while (area > 4.0 * M_PI) area -= 4.0 * M_PI;
946 }
947
948 // finally here we will handle the single separate patches...
949 // finally here we will handle the single separate patches...
950 // finally here we will handle the single separate patches...
951
952 for (i32s n2 = 0;n2 < coi_count;n2++)
953 {
954 if (coit[n2].flag) continue;
955
956 f64 t1a = 2.0 * M_PI / radius1[n1];
957 area -= t1a * (radius1[n1] - coit[n2].g);
958
959 if (p1 > 0)
960 {
961 for (i32s n3 = 0;n3 < 3;n3++)
962 {
963 f64 t1b = t1a * coit[n2].dg[n3];
964
965 d_sasa[coit[n2].index * 3 + n3] += t1b;
966 d_sasa[n1 * 3 + n3] -= t1b;
967 }
968 }
969 }
970
971 sasa[n1] = area;
972 }
973 }
974
975 /*################################################################################################*/
976
977 // eof
978
979