1 //
2 // Copyright (c) 2016, Guillaume GODIN
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without
6 // modification, are permitted provided that the following conditions are
7 // met:
8 //
9 // * Redistributions of source code must retain the above copyright
10 // notice, this list of conditions and the following disclaimer.
11 // * Redistributions in binary form must reproduce the above
12 // copyright notice, this list of conditions and the following
13 // disclaimer in the documentation and/or other materials provided
14 // with the distribution.
15 // * Neither the name of Institue of Cancer Research.
16 // nor the names of its contributors may be used to endorse or promote
17 // products derived from this software without specific prior written
18 // permission.
19 //
20 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 //
32 // Adding GETAWAY descriptors by Guillaume Godin
33 // for build & set RDBASE! => export RDBASE=/Users/mbp/Github/rdkit_mine/
34
35 #include <GraphMol/RDKitBase.h>
36 #include <GraphMol/MolTransforms/MolTransforms.h>
37
38 #include "GETAWAY.h"
39 #include "PBF.h"
40 #include "MolData3Ddescriptors.h"
41
42 #include "GraphMol/PartialCharges/GasteigerCharges.h"
43 #include "GraphMol/PartialCharges/GasteigerParams.h"
44 #include <Numerics/EigenSolvers/PowerEigenSolver.h>
45
46 #include <Numerics/Matrix.h>
47 #include <Numerics/SquareMatrix.h>
48 #include <Numerics/SymmMatrix.h>
49 #include <cmath>
50 #include <Eigen/Dense>
51 #include <Eigen/SVD>
52 #include <iostream>
53 #include <deque>
54 #include <Eigen/Core>
55 #include <Eigen/QR>
56
57 using namespace Eigen;
58
59 namespace RDKit {
60
61 namespace Descriptors {
62
63 namespace {
64
65 MolData3Ddescriptors moldata3D;
66
roundn(double in,int factor)67 double roundn(double in, int factor) {
68 return std::round(in * pow(10., factor)) / pow(10., factor);
69 }
70
retreiveMat(MatrixXd matrix)71 double* retreiveMat(MatrixXd matrix) {
72 double* arrayd = matrix.data();
73 return arrayd;
74 }
75
retreiveVect(VectorXd matrix)76 double* retreiveVect(VectorXd matrix) {
77 double* arrayd = matrix.data();
78 return arrayd;
79 }
80
getEigenVect(std::vector<double> v)81 VectorXd getEigenVect(std::vector<double> v) {
82 double* varray_ptr = &v[0];
83 Map<VectorXd> V(varray_ptr, v.size());
84 return V;
85 }
86
round_to_n_digits_(double x,int n)87 double round_to_n_digits_(double x, int n) {
88 char buff[32];
89 sprintf(buff, "%.*g", n, x); // use g to have significant digits!
90 return atof(buff);
91 }
92
round_to_n_digits(double x,int n)93 double round_to_n_digits(double x, int n) {
94 double scale = pow(10.0, ceil(log10(fabs(x))) + n);
95 return std::round(x * scale) / scale;
96 }
97
IsClose(double a,double b,int n)98 bool IsClose(double a, double b, int n) {
99 return (fabs(a - b) <= pow(0.1, n) * 1.1);
100 }
101
countZeros(std::string ta)102 int countZeros(std::string ta) {
103 int nbzero = 0;
104 for (char i : ta) {
105 if (i != '0') {
106 break;
107 }
108 nbzero = nbzero + 1;
109 }
110
111 return nbzero;
112 }
113
IsClose2(double a,double b,unsigned int precision)114 bool IsClose2(double a, double b, unsigned int precision) {
115 RDUNUSED_PARAM(precision);
116 bool isclose = false;
117
118 std::string sa, sb;
119 std::string ta, tb;
120 std::stringstream outa, outb;
121 outa << a;
122 sa = outa.str();
123 outb << b;
124 sb = outb.str();
125 ta = sa.substr(sa.find('.') + 1);
126 tb = sb.substr(sb.find('.') + 1);
127
128 // same number of digits!
129 if (countZeros(ta) == countZeros(tb)) {
130 if (ta.length() > tb.length()) {
131 tb += '0';
132 }
133 if (ta.length() < tb.length()) {
134 ta += '0';
135 }
136 }
137 // std::cout << tb << ":"<< atoi(tb.c_str()) << ":" << countZeros(ta) <<
138 // "\n";
139 // std::cout << ta << ":"<< atoi(ta.c_str()) << ":" << countZeros(ta) <<
140 // "\n";
141
142 if (ta.length() == tb.length()) {
143 // last digit +/-1 deviation only!)
144 if (abs(atoi(ta.c_str()) - atoi(tb.c_str())) < 2) {
145 // std::cout << a << "," << b << " ta:" << ta << "," << "tb:" << tb<<
146 // "\n";
147 isclose = true;
148 }
149 }
150 return isclose;
151 }
152
IsClose3(double a,double b,unsigned int precision)153 bool IsClose3(double a, double b, unsigned int precision) {
154 RDUNUSED_PARAM(precision);
155 bool isclose = false;
156 std::string sa, sb;
157 std::string ta, tb;
158 std::stringstream outa, outb;
159 outa << a;
160 sa = outa.str();
161 outb << b;
162 sb = outb.str();
163 ta = sa.substr(sa.find('.') + 1);
164 tb = sb.substr(sb.find('.') + 1);
165 // same number of digits!
166 if (ta.length() == tb.length()) {
167 // last digit +/-1 deviation only!)
168 if (abs(atoi(ta.c_str()) - atoi(tb.c_str())) < 2) {
169 // std::cout << a << "," << b << " ta:" << ta << "," << "tb:" << tb<<
170 // "\n";
171 isclose = true;
172 }
173 }
174 return isclose;
175 }
176
177 // need to clean that code to have always the same output which is not the case
178 // the classes used the last 2 digit +/- 1 to cluster in the same group (with 3
179 // digit precision)
180 // 0.012 and 0.013 are in the same group while 0.014 are not!
clusterArray(std::vector<double> data,double precision)181 std::vector<double> clusterArray(std::vector<double> data, double precision) {
182 std::vector<double> Store;
183
184 // sort the input data
185 std::sort(data.begin(), data.end());
186
187 // find the difference between each number and its predecessor
188 std::vector<double> diffs;
189
190 std::adjacent_difference(data.begin(), data.end(), std::back_inserter(diffs));
191
192 // convert differences into percentage changes
193 // std::transform(diffs.begin(), diffs.end(), data.begin(), diffs.begin(),
194 // std::divides<double>());
195
196 int j = 0;
197 int count = 0;
198 for (unsigned int i = 0; i < data.size(); i++) {
199 // std::cout << diffs[i] << ",";
200 count++;
201 // if a difference exceeds 0.01 <=> 1%, start a new group: if transform is
202 // used!
203 // use diff not ratio (with 0.003 precision look like it's what it's used in
204 // Dragon 6!)
205 if (diffs[i] > pow(0.1, precision)) {
206 Store.push_back(count);
207 count = 0;
208 j++;
209 }
210 }
211
212 return Store;
213 }
214
215 // need to clean that code to have always the same output which is not the case
216 // the classes used the last 2 digit +/- 1 to cluster in the same group (with 3
217 // digit precision)
218 // 0.012 and 0.013 are in the same group while 0.014 are not!
219 // alternative clustering method not working as in Dragon
clusterArray2(std::vector<double> data,unsigned int precision)220 std::vector<double> clusterArray2(std::vector<double> data,
221 unsigned int precision) {
222 std::vector<double> Store;
223
224 // sort the input data descend order!
225 std::sort(data.begin(), data.end(), std::greater<double>());
226 std::deque<double> B(data.begin(), data.end());
227
228 int count = 0;
229 while (!B.empty()) {
230 double dk = B.front();
231 B.pop_front();
232 count++;
233 // case where the last value is not grouped it goes to the single last group
234 if (B.empty()) {
235 Store.push_back(count);
236 }
237
238 for (unsigned int i = 0; i < B.size(); i++) {
239 if (IsClose2(dk, B[i], precision)) {
240 count++;
241 } else {
242 Store.push_back(count);
243 for (unsigned int j = 0; j < i; j++) {
244 B.pop_front();
245 }
246 count = 0;
247 break;
248 }
249
250 if (i == B.size() - 1) {
251 Store.push_back(count);
252 B.pop_front();
253 count = 0;
254 }
255 }
256 }
257 return Store;
258 }
259
GetGeodesicMatrix(const double * dist,int lag,int numAtoms)260 std::vector<double> GetGeodesicMatrix(const double* dist, int lag,
261 int numAtoms) {
262 PRECONDITION(dist != nullptr, "bad array");
263
264 int sizeArray = numAtoms * numAtoms;
265 std::vector<double> Geodesic;
266 Geodesic.reserve(sizeArray);
267 std::transform(dist, dist + sizeArray, Geodesic.begin(),
268 [&lag](const double& dist) { return int(dist == lag); });
269
270 return Geodesic;
271 }
272
getSVD(const MatrixXd & A)273 JacobiSVD<MatrixXd> getSVD(const MatrixXd& A) {
274 JacobiSVD<MatrixXd> mysvd(A, ComputeThinU | ComputeThinV);
275 return mysvd;
276 }
277
GetPinv(const MatrixXd & A)278 MatrixXd GetPinv(const MatrixXd& A) {
279 JacobiSVD<MatrixXd> svd = getSVD(A);
280 double pinvtoler = 1.e-3; // choose your tolerance wisely!
281 VectorXd vs = svd.singularValues();
282 VectorXd vsinv = svd.singularValues();
283
284 for (unsigned int i = 0; i < A.cols(); ++i) {
285 if (vs(i) > pinvtoler) {
286 vsinv(i) = 1.0 / vs(i);
287 } else {
288 vsinv(i) = 0.0;
289 }
290 }
291
292 MatrixXd S = vsinv.asDiagonal();
293 MatrixXd Ap = svd.matrixV() * S * svd.matrixU().transpose();
294 return Ap;
295 }
296
GetCenterMatrix(MatrixXd Mat)297 MatrixXd GetCenterMatrix(MatrixXd Mat) {
298 VectorXd v = Mat.colwise().mean();
299 MatrixXd X = Mat.rowwise() - v.transpose();
300 return X;
301 }
302
GetHmatrix(MatrixXd X)303 MatrixXd GetHmatrix(MatrixXd X) {
304 MatrixXd Weighted = X.transpose() * X;
305 return X * GetPinv(Weighted) * X.transpose();
306 }
307
GetRmatrix(MatrixXd H,MatrixXd DM,int numAtoms)308 MatrixXd GetRmatrix(MatrixXd H, MatrixXd DM, int numAtoms) {
309 MatrixXd R = MatrixXd::Zero(numAtoms, numAtoms);
310 for (int i = 0; i < numAtoms - 1; i++) {
311 for (int j = i + 1; j < numAtoms; j++) {
312 R(i, j) = sqrt(H(i, i) * H(j, j)) / DM(i, j);
313 R(j, i) = R(i, j);
314 }
315 }
316
317 return R;
318 }
319
GetHeavyList(const ROMol & mol)320 std::vector<int> GetHeavyList(const ROMol& mol) {
321 int numAtoms = mol.getNumAtoms();
322 std::vector<int> HeavyList;
323 HeavyList.reserve(numAtoms);
324 for (int i = 0; i < numAtoms; ++i) {
325 const RDKit::Atom* atom = mol.getAtomWithIdx(i);
326 HeavyList.push_back(int(atom->getAtomicNum() > 1));
327 }
328 return HeavyList;
329 }
330
AppendDouble(double * w,const double * Append,int length,int pos)331 double* AppendDouble(double* w, const double* Append, int length, int pos) {
332 PRECONDITION(w != nullptr, "bad array");
333 PRECONDITION(Append != nullptr, "bad array");
334
335 for (int i = pos; i < pos + length; i++) {
336 w[i] = Append[i - pos];
337 }
338
339 return w;
340 }
341
getRCON(MatrixXd R,MatrixXd Adj,int numAtoms)342 double getRCON(MatrixXd R, MatrixXd Adj, int numAtoms) {
343 // similar implementation of J. Chem. Inf. Comput. Sci. 2004, 44, 200-209
344 // equation 1 or 2 page 201
345 // we use instead of atomic absolute values the atomic relative ones as in
346 // Dragon
347 double RCON = 0.0;
348 VectorXd VSR = R.rowwise().sum();
349 for (int i = 0; i < numAtoms - 1; ++i) {
350 for (int j = i + 1; j < numAtoms; ++j) {
351 if (Adj(i, j) > 0) {
352 RCON +=
353 sqrt(VSR(i) * VSR(j)); // the sqrt is in the sum not over the sum!
354 }
355 }
356 }
357 return RCON;
358 }
359
getHATS(double W1,double W2,double H1,double H2)360 double getHATS(double W1, double W2, double H1, double H2) {
361 return W1 * H1 * W2 * H2;
362 }
363
getH(double W1,double W2,double H)364 double getH(double W1, double W2, double H) { return W1 * H * W2; }
365
getMax(const double * Rk)366 double getMax(const double* Rk) {
367 PRECONDITION(Rk != nullptr, "bad rK");
368 double RTp = 0;
369 for (int j = 0; j < 8; j++) {
370 if (Rk[j] > RTp) {
371 RTp = Rk[j];
372 }
373 }
374 return RTp;
375 }
376
getGETAWAYDescCustom(MatrixXd H,MatrixXd R,MatrixXd Adj,int numAtoms,std::vector<int> Heavylist,const ROMol & mol,std::vector<double> & res,unsigned int precision,const std::string & customAtomPropName)377 void getGETAWAYDescCustom(MatrixXd H, MatrixXd R, MatrixXd Adj, int numAtoms,
378 std::vector<int> Heavylist, const ROMol& mol,
379 std::vector<double>& res, unsigned int precision,
380 const std::string& customAtomPropName) {
381 // prepare data for Getaway parameter computation
382 // compute parameters
383
384 VectorXd Lev = H.diagonal();
385 std::vector<double> heavyLev;
386
387 for (int i = 0; i < numAtoms; i++) {
388 if (Heavylist[i] == 1) {
389 heavyLev.push_back(round_to_n_digits_(Lev(i), precision));
390 }
391 }
392
393 std::vector<double> Clus = clusterArray2(heavyLev, precision);
394
395 double numHeavy = heavyLev.size();
396 double ITH0 = numHeavy * log(numHeavy) / log(2.0);
397 double ITH = ITH0;
398 for (double Clu : Clus) {
399 ITH -= Clu * log(Clu) / log(2.0);
400 }
401 res[0] = roundn(ITH, 3); // issue sometime with this due to cluster
402 double ISH = ITH / ITH0;
403 res[1] = roundn(ISH, 3); // issue sometime with this due to cluster
404
405 // use the PBF to determine 2D vs 3D (with Threshold)
406 // determine if this is a plane molecule
407 double pbf = RDKit::Descriptors::PBF(mol);
408 double D;
409 if (pbf < 1.e-5) {
410 D = 2.0;
411 } else {
412 D = 3.0;
413 }
414
415 // what about linear molecule ie (D=1) ?
416 double HIC = 0.0;
417 for (int i = 0; i < numAtoms; i++) {
418 HIC -= H(i, i) / D * log(H(i, i) / D) / log(2);
419 }
420 res[2] = roundn(HIC, 3);
421
422 double HGM = 1.0;
423 for (int i = 0; i < numAtoms; i++) {
424 HGM = HGM * H(i, i);
425 }
426 HGM = 100.0 * pow(HGM, 1.0 / numAtoms);
427 res[3] = roundn(HGM, 3);
428
429 double RARS = R.rowwise().sum().sum() / numAtoms;
430
431 JacobiSVD<MatrixXd> mysvd = getSVD(R);
432
433 VectorXd EIG = mysvd.singularValues();
434
435 double rcon = getRCON(R, std::move(Adj), numAtoms);
436
437 std::vector<double> customAtomArray =
438 moldata3D.GetCustomAtomProp(mol, customAtomPropName);
439 VectorXd Wc = getEigenVect(customAtomArray);
440
441 MatrixXd Bi;
442 MatrixXd RBw;
443 double HATSc, HATSct, H0ct, R0ct, H0c, R0c, Rkmaxc, tmpc;
444 double HATSk[9];
445 double Hk[9];
446 double Rk[8];
447 double Rp[8];
448
449 double* dist =
450 MolOps::getDistanceMat(mol, false); // need to be be set to false to have
451 // topological distance not weighted!
452
453 Map<MatrixXd> D2(dist, numAtoms, numAtoms);
454
455 double Dmax = D2.colwise().maxCoeff().maxCoeff();
456
457 HATSct = 0.0;
458 H0ct = 0.0;
459 R0ct = 0.0;
460
461 // need to loop other all the D values so <= not <!
462 for (int i = 0; i <= Dmax; i++) {
463 if (i == 0) {
464 Bi = H.diagonal().asDiagonal();
465 }
466
467 HATSc = 0.0;
468 H0c = 0.0;
469
470 if (i == 0) { // use Bi
471 for (int j = 0; j < numAtoms; j++) {
472 for (int k = j; k < numAtoms; k++) {
473 if (Bi(j, k) > 0) {
474 HATSc += getHATS((double)Wc(j), (double)Wc(j), (double)H(j, j),
475 (double)H(j, j));
476
477 if (H(j, k) > 0) {
478 H0c += getH((double)Wc(j), (double)Wc(k), (double)H(j, k));
479 }
480 }
481 }
482 }
483 }
484 auto Bimat = GetGeodesicMatrix(dist, i, numAtoms);
485 Map<MatrixXd> Bj(Bimat.data(), numAtoms, numAtoms);
486 if (i > 0 && i < 9) { // use Bj
487 for (int j = 0; j < numAtoms - 1; j++) {
488 for (int k = j + 1; k < numAtoms; k++) {
489 if (Bj(j, k) == 1) {
490 HATSc += getHATS((double)Wc(j), (double)Wc(k), (double)H(j, j),
491 (double)H(k, k));
492
493 if (H(j, k) > 0) {
494 H0c += getH((double)Wc(j), (double)Wc(k), (double)H(j, k));
495 }
496 }
497 }
498 }
499 }
500
501 if (i >= 9) { // Totals missing part
502 for (int j = 0; j < numAtoms - 1; j++) {
503 for (int k = j + 1; k < numAtoms; k++) {
504 if (Bj(j, k) == 1) {
505 HATSct += 2 * getHATS((double)Wc(j), (double)Wc(k), (double)H(j, j),
506 (double)H(k, k));
507
508 if (H(j, k) > 0) {
509 H0ct += 2 * getH((double)Wc(j), (double)Wc(k), (double)H(j, k));
510 }
511 }
512 }
513 }
514 }
515
516 if (i < 9) {
517 HATSk[i] = HATSc;
518 Hk[i] = H0c;
519 }
520
521 R0c = 0.0;
522 Rkmaxc = 0.0;
523
524 // from 1 to 8
525 if (i > 0 && i < 9) {
526 for (int j = 0; j < numAtoms - 1; j++) {
527 for (int k = j + 1; k < numAtoms; k++) {
528 if (Bj(j, k) == 1) {
529 tmpc = getH((double)Wc(j), (double)Wc(k),
530 (double)R(j, k)); // Use same function but on all R not
531 // "H>0" like in the previous loop &
532 // i>0!
533
534 R0c += tmpc;
535
536 if (tmpc > Rkmaxc) {
537 Rkmaxc = tmpc;
538 }
539 }
540 }
541 }
542
543 Rk[i - 1] = R0c;
544 Rp[i - 1] = Rkmaxc;
545 }
546
547 if (i >= 9) {
548 for (int j = 0; j < numAtoms - 1; j++) {
549 for (int k = j + 1; k < numAtoms; k++) {
550 if (Bj(j, k) == 1) {
551 R0ct += 2 * getH((double)Wc(j), (double)Wc(k),
552 (double)R(j, k)); // Use same function but on all
553 // R not "H>0" like in the
554 // previous loop & i>0!
555 }
556 }
557 }
558 }
559 }
560
561 // can be column vs row selected that can explain the issue!
562 double HATSTc = HATSk[0] + HATSct;
563 for (int ii = 1; ii < 9; ii++) {
564 HATSTc += 2 * HATSk[ii];
565 }
566
567 double HTc = Hk[0] + H0ct;
568 for (int ii = 1; ii < 9; ii++) {
569 HTc += 2.0 * Hk[ii];
570 }
571
572 double RTc = 0.0;
573 for (double rk_ii : Rk) {
574 RTc += 2.0 * rk_ii;
575 }
576
577 double RTMc = getMax(Rp);
578
579 // create the output vector...
580 for (int i = 0; i < 9; i++) {
581 res[i + 4] = roundn(Hk[i], 3);
582 res[i + 14] = roundn(HATSk[i], 3);
583 }
584
585 res[13] = roundn(HTc, 3);
586 res[23] = roundn(HATSTc, 3);
587
588 res[24] = roundn(rcon, 3);
589 res[25] = roundn(RARS, 3);
590 res[26] = roundn(EIG(0), 3);
591
592 for (int i = 0; i < 8; i++) {
593 res[i + 27] = roundn(Rk[i], 3);
594 res[i + 36] = roundn(Rp[i], 3);
595 }
596
597 res[35] = roundn(RTc + R0ct, 3);
598 res[44] = roundn(RTMc, 3);
599 }
600
getGETAWAYDesc(MatrixXd H,MatrixXd R,MatrixXd Adj,int numAtoms,std::vector<int> Heavylist,const ROMol & mol,std::vector<double> & res,unsigned int precision)601 void getGETAWAYDesc(MatrixXd H, MatrixXd R, MatrixXd Adj, int numAtoms,
602 std::vector<int> Heavylist, const ROMol& mol,
603 std::vector<double>& res, unsigned int precision) {
604 // prepare data for Getaway parameter computation
605 // compute parameters
606
607 VectorXd Lev = H.diagonal();
608 std::vector<double> heavyLev;
609
610 for (int i = 0; i < numAtoms; i++) {
611 if (Heavylist[i] == 1) {
612 heavyLev.push_back(round_to_n_digits_(Lev(i), precision));
613 }
614 }
615
616 std::vector<double> Clus = clusterArray2(heavyLev, precision);
617
618 double numHeavy = heavyLev.size();
619 double ITH0 = numHeavy * log(numHeavy) / log(2.0);
620 double ITH = ITH0;
621 for (double Clu : Clus) {
622 ITH -= Clu * log(Clu) / log(2.0);
623 }
624 res[0] = roundn(ITH, 3); // issue sometime with this due to cluster
625 double ISH = ITH / ITH0;
626 res[1] = roundn(ISH, 3); // issue sometime with this due to cluster
627
628 // use the PBF to determine 2D vs 3D (with Threshold)
629 // determine if this is a plane molecule
630 double pbf = RDKit::Descriptors::PBF(mol);
631 double D;
632 if (pbf < 1.e-5) {
633 D = 2.0;
634 } else {
635 D = 3.0;
636 }
637
638 // what about linear molecule ie (D=1) ?
639 double HIC = 0.0;
640 for (int i = 0; i < numAtoms; i++) {
641 HIC -= H(i, i) / D * log(H(i, i) / D) / log(2.);
642 }
643 res[2] = roundn(HIC, 3);
644
645 double HGM = 1.0;
646 for (int i = 0; i < numAtoms; i++) {
647 HGM = HGM * H(i, i);
648 }
649 HGM = 100.0 * pow(HGM, 1.0 / numAtoms);
650 res[3] = roundn(HGM, 3);
651
652 double RARS = R.rowwise().sum().sum() / numAtoms;
653
654 JacobiSVD<MatrixXd> mysvd = getSVD(R);
655
656 VectorXd EIG = mysvd.singularValues();
657
658 double rcon = getRCON(R, std::move(Adj), numAtoms);
659
660 std::vector<double> wp = moldata3D.GetRelativePol(mol);
661
662 VectorXd Wp = getEigenVect(wp);
663
664 std::vector<double> wm = moldata3D.GetRelativeMW(mol);
665
666 VectorXd Wm = getEigenVect(wm);
667
668 std::vector<double> wi = moldata3D.GetRelativeIonPol(mol);
669
670 VectorXd Wi = getEigenVect(wi);
671
672 std::vector<double> wv = moldata3D.GetRelativeVdW(mol);
673
674 VectorXd Wv = getEigenVect(wv);
675
676 std::vector<double> we = moldata3D.GetRelativeENeg(mol);
677
678 VectorXd We = getEigenVect(we);
679
680 std::vector<double> wu = moldata3D.GetUn(numAtoms);
681
682 VectorXd Wu = getEigenVect(wu);
683
684 std::vector<double> ws = moldata3D.GetIState(mol);
685
686 VectorXd Ws = getEigenVect(ws);
687
688 MatrixXd Bi;
689 MatrixXd RBw;
690 double HATSu, HATSm, HATSv, HATSe, HATSp, HATSi, HATSs;
691 double HATSut, HATSmt, HATSvt, HATSet, HATSpt, HATSit, HATSst;
692 double H0ut, H0mt, H0vt, H0et, H0pt, H0it, H0st;
693 double R0ut, R0mt, R0vt, R0et, R0pt, R0it, R0st;
694 double H0u, H0m, H0v, H0e, H0p, H0i, H0s;
695 double R0u, R0m, R0v, R0e, R0p, R0i, R0s;
696 double Rkmaxu, Rkmaxm, Rkmaxv, Rkmaxe, Rkmaxp, Rkmaxi, Rkmaxs;
697 double tmpu, tmpm, tmpv, tmpe, tmpp, tmpi, tmps;
698 double HATSk[7][9];
699 double Hk[7][9];
700 double Rk[7][8];
701 double Rp[7][8];
702
703 double* dist =
704 MolOps::getDistanceMat(mol, false); // need to be be set to false to have
705 // topological distance not weighted!
706
707 Map<MatrixXd> D2(dist, numAtoms, numAtoms);
708
709 double Dmax = D2.colwise().maxCoeff().maxCoeff();
710
711 HATSut = 0.0;
712 HATSmt = 0.0;
713 HATSvt = 0.0;
714 HATSet = 0.0;
715 HATSpt = 0.0;
716 HATSit = 0.0;
717 HATSst = 0.0;
718
719 H0ut = 0.0;
720 H0mt = 0.0;
721 H0vt = 0.0;
722 H0et = 0.0;
723 H0pt = 0.0;
724 H0it = 0.0;
725 H0st = 0.0;
726
727 R0ut = 0.0;
728 R0mt = 0.0;
729 R0vt = 0.0;
730 R0et = 0.0;
731 R0pt = 0.0;
732 R0it = 0.0;
733 R0st = 0.0;
734
735 // need to loop other all the D values so <= not <!
736 for (int i = 0; i <= Dmax; i++) {
737 if (i == 0) {
738 Bi = H.diagonal().asDiagonal();
739 }
740
741 HATSu = 0.0;
742 HATSm = 0.0;
743 HATSv = 0.0;
744 HATSe = 0.0;
745 HATSp = 0.0;
746 HATSi = 0.0;
747 HATSs = 0.0;
748
749 H0u = 0.0;
750 H0m = 0.0;
751 H0v = 0.0;
752 H0e = 0.0;
753 H0p = 0.0;
754 H0i = 0.0;
755 H0s = 0.0;
756
757 if (i == 0) { // use Bi
758 for (int j = 0; j < numAtoms; j++) {
759 for (int k = j; k < numAtoms; k++) {
760 if (Bi(j, k) > 0) {
761 HATSu += getHATS((double)Wu(j), (double)Wu(j), (double)H(j, j),
762 (double)H(j, j));
763 HATSm += getHATS((double)Wm(j), (double)Wm(j), (double)H(j, j),
764 (double)H(j, j));
765 HATSv += getHATS((double)Wv(j), (double)Wv(j), (double)H(j, j),
766 (double)H(j, j));
767 HATSe += getHATS((double)We(j), (double)We(j), (double)H(j, j),
768 (double)H(j, j));
769 HATSp += getHATS((double)Wp(j), (double)Wp(j), (double)H(j, j),
770 (double)H(j, j));
771 HATSi += getHATS((double)Wi(j), (double)Wi(j), (double)H(j, j),
772 (double)H(j, j));
773 HATSs += getHATS((double)Ws(j), (double)Ws(j), (double)H(j, j),
774 (double)H(j, j));
775
776 if (H(j, k) > 0) {
777 H0u += getH((double)Wu(j), (double)Wu(k), (double)H(j, k));
778 H0m += getH((double)Wm(j), (double)Wm(k), (double)H(j, k));
779 H0v += getH((double)Wv(j), (double)Wv(k), (double)H(j, k));
780 H0e += getH((double)We(j), (double)We(k), (double)H(j, k));
781 H0p += getH((double)Wp(j), (double)Wp(k), (double)H(j, k));
782 H0i += getH((double)Wi(j), (double)Wi(k), (double)H(j, k));
783 H0s += getH((double)Ws(j), (double)Ws(k), (double)H(j, k));
784 }
785 }
786 }
787 }
788 }
789 auto Bimat = GetGeodesicMatrix(dist, i, numAtoms);
790 Map<MatrixXd> Bj(Bimat.data(), numAtoms, numAtoms);
791 if (i > 0 && i < 9) { // use Bj
792 for (int j = 0; j < numAtoms - 1; j++) {
793 for (int k = j + 1; k < numAtoms; k++) {
794 if (Bj(j, k) == 1) {
795 HATSu += getHATS((double)Wu(j), (double)Wu(k), (double)H(j, j),
796 (double)H(k, k));
797 HATSm += getHATS((double)Wm(j), (double)Wm(k), (double)H(j, j),
798 (double)H(k, k));
799 HATSv += getHATS((double)Wv(j), (double)Wv(k), (double)H(j, j),
800 (double)H(k, k));
801 HATSe += getHATS((double)We(j), (double)We(k), (double)H(j, j),
802 (double)H(k, k));
803 HATSp += getHATS((double)Wp(j), (double)Wp(k), (double)H(j, j),
804 (double)H(k, k));
805 HATSi += getHATS((double)Wi(j), (double)Wi(k), (double)H(j, j),
806 (double)H(k, k));
807 HATSs += getHATS((double)Ws(j), (double)Ws(k), (double)H(j, j),
808 (double)H(k, k));
809
810 if (H(j, k) > 0) {
811 H0u += getH((double)Wu(j), (double)Wu(k), (double)H(j, k));
812 H0m += getH((double)Wm(j), (double)Wm(k), (double)H(j, k));
813 H0v += getH((double)Wv(j), (double)Wv(k), (double)H(j, k));
814 H0e += getH((double)We(j), (double)We(k), (double)H(j, k));
815 H0p += getH((double)Wp(j), (double)Wp(k), (double)H(j, k));
816 H0i += getH((double)Wi(j), (double)Wi(k), (double)H(j, k));
817 H0s += getH((double)Ws(j), (double)Ws(k), (double)H(j, k));
818 }
819 }
820 }
821 }
822 }
823
824 if (i >= 9) { // Totals missing part
825 for (int j = 0; j < numAtoms - 1; j++) {
826 for (int k = j + 1; k < numAtoms; k++) {
827 if (Bj(j, k) == 1) {
828 HATSut += 2 * getHATS((double)Wu(j), (double)Wu(k), (double)H(j, j),
829 (double)H(k, k));
830 HATSmt += 2 * getHATS((double)Wm(j), (double)Wm(k), (double)H(j, j),
831 (double)H(k, k));
832 HATSvt += 2 * getHATS((double)Wv(j), (double)Wv(k), (double)H(j, j),
833 (double)H(k, k));
834 HATSet += 2 * getHATS((double)We(j), (double)We(k), (double)H(j, j),
835 (double)H(k, k));
836 HATSpt += 2 * getHATS((double)Wp(j), (double)Wp(k), (double)H(j, j),
837 (double)H(k, k));
838 HATSit += 2 * getHATS((double)Wi(j), (double)Wi(k), (double)H(j, j),
839 (double)H(k, k));
840 HATSst += 2 * getHATS((double)Ws(j), (double)Ws(k), (double)H(j, j),
841 (double)H(k, k));
842
843 if (H(j, k) > 0) {
844 H0ut += 2 * getH((double)Wu(j), (double)Wu(k), (double)H(j, k));
845 H0mt += 2 * getH((double)Wm(j), (double)Wm(k), (double)H(j, k));
846 H0vt += 2 * getH((double)Wv(j), (double)Wv(k), (double)H(j, k));
847 H0et += 2 * getH((double)We(j), (double)We(k), (double)H(j, k));
848 H0pt += 2 * getH((double)Wp(j), (double)Wp(k), (double)H(j, k));
849 H0it += 2 * getH((double)Wi(j), (double)Wi(k), (double)H(j, k));
850 H0st += 2 * getH((double)Ws(j), (double)Ws(k), (double)H(j, k));
851 }
852 }
853 }
854 }
855 }
856
857 if (i < 9) {
858 HATSk[0][i] = HATSu;
859 HATSk[1][i] = HATSm;
860 HATSk[2][i] = HATSv;
861 HATSk[3][i] = HATSe;
862 HATSk[4][i] = HATSp;
863 HATSk[5][i] = HATSi;
864 HATSk[6][i] = HATSs;
865
866 Hk[0][i] = H0u;
867 Hk[1][i] = H0m;
868 Hk[2][i] = H0v;
869 Hk[3][i] = H0e;
870 Hk[4][i] = H0p;
871 Hk[5][i] = H0i;
872 Hk[6][i] = H0s;
873 }
874
875 R0u = 0.0;
876 R0m = 0.0;
877 R0v = 0.0;
878 R0e = 0.0;
879 R0p = 0.0;
880 R0i = 0.0;
881 R0s = 0.0;
882
883 Rkmaxu = 0;
884 Rkmaxm = 0;
885 Rkmaxv = 0;
886 Rkmaxe = 0;
887 Rkmaxp = 0;
888 Rkmaxi = 0;
889 Rkmaxs = 0;
890
891 // from 1 to 8
892 if (i > 0 && i < 9) {
893 for (int j = 0; j < numAtoms - 1; j++) {
894 for (int k = j + 1; k < numAtoms; k++) {
895 if (Bj(j, k) == 1) {
896 tmpu = getH((double)Wu(j), (double)Wu(k),
897 (double)R(j, k)); // Use same function but on all R not
898 // "H>0" like in the previous loop &
899 // i>0!
900 tmpm = getH((double)Wm(j), (double)Wm(k),
901 (double)R(j, k)); // Use same function but on all R not
902 // "H>0" like in the previous loop &
903 // i>0!
904 tmpv = getH((double)Wv(j), (double)Wv(k),
905 (double)R(j, k)); // Use same function but on all R not
906 // "H>0" like in the previous loop &
907 // i>0!
908 tmpe = getH((double)We(j), (double)We(k),
909 (double)R(j, k)); // Use same function but on all R not
910 // "H>0" like in the previous loop &
911 // i>0!
912 tmpp = getH((double)Wp(j), (double)Wp(k),
913 (double)R(j, k)); // Use same function but on all R not
914 // "H>0" like in the previous loop &
915 // i>0!
916 tmpi = getH((double)Wi(j), (double)Wi(k),
917 (double)R(j, k)); // Use same function but on all R not
918 // "H>0" like in the previous loop &
919 // i>0!
920 tmps = getH((double)Ws(j), (double)Ws(k),
921 (double)R(j, k)); // Use same function but on all R not
922 // "H>0" like in the previous loop &
923 // i>0!
924 R0u += tmpu;
925 R0m += tmpm;
926 R0v += tmpv;
927 R0e += tmpe;
928 R0p += tmpp;
929 R0i += tmpi;
930 R0s += tmps;
931 if (tmpu > Rkmaxu) {
932 Rkmaxu = tmpu;
933 }
934 if (tmpm > Rkmaxm) {
935 Rkmaxm = tmpm;
936 }
937 if (tmpv > Rkmaxv) {
938 Rkmaxv = tmpv;
939 }
940 if (tmpe > Rkmaxe) {
941 Rkmaxe = tmpe;
942 }
943 if (tmpp > Rkmaxp) {
944 Rkmaxp = tmpp;
945 }
946 if (tmpi > Rkmaxi) {
947 Rkmaxi = tmpi;
948 }
949 if (tmps > Rkmaxs) {
950 Rkmaxs = tmps;
951 }
952 }
953 }
954 }
955
956 Rk[0][i - 1] = R0u;
957 Rk[1][i - 1] = R0m;
958 Rk[2][i - 1] = R0v;
959 Rk[3][i - 1] = R0e;
960 Rk[4][i - 1] = R0p;
961 Rk[5][i - 1] = R0i;
962 Rk[6][i - 1] = R0s;
963
964 Rp[0][i - 1] = Rkmaxu;
965 Rp[1][i - 1] = Rkmaxm;
966 Rp[2][i - 1] = Rkmaxv;
967 Rp[3][i - 1] = Rkmaxe;
968 Rp[4][i - 1] = Rkmaxp;
969 Rp[5][i - 1] = Rkmaxi;
970 Rp[6][i - 1] = Rkmaxs;
971 }
972
973 if (i >= 9) {
974 for (int j = 0; j < numAtoms - 1; j++) {
975 for (int k = j + 1; k < numAtoms; k++) {
976 if (Bj(j, k) == 1) {
977 R0ut += 2 * getH((double)Wu(j), (double)Wu(k),
978 (double)R(j, k)); // Use same function but on all
979 // R not "H>0" like in the
980 // previous loop & i>0!
981 R0mt += 2 * getH((double)Wm(j), (double)Wm(k),
982 (double)R(j, k)); // Use same function but on all
983 // R not "H>0" like in the
984 // previous loop & i>0!
985 R0vt += 2 * getH((double)Wv(j), (double)Wv(k),
986 (double)R(j, k)); // Use same function but on all
987 // R not "H>0" like in the
988 // previous loop & i>0!
989 R0et += 2 * getH((double)We(j), (double)We(k),
990 (double)R(j, k)); // Use same function but on all
991 // R not "H>0" like in the
992 // previous loop & i>0!
993 R0pt += 2 * getH((double)Wp(j), (double)Wp(k),
994 (double)R(j, k)); // Use same function but on all
995 // R not "H>0" like in the
996 // previous loop & i>0!
997 R0it += 2 * getH((double)Wi(j), (double)Wi(k),
998 (double)R(j, k)); // Use same function but on all
999 // R not "H>0" like in the
1000 // previous loop & i>0!
1001 R0st += 2 * getH((double)Ws(j), (double)Ws(k),
1002 (double)R(j, k)); // Use same function but on all
1003 // R not "H>0" like in the
1004 // previous loop & i>0!
1005 }
1006 }
1007 }
1008 }
1009 }
1010
1011 // can be column vs row selected that can explain the issue!
1012 double HATSTu = HATSk[0][0] + HATSut;
1013 double HATSTm = HATSk[1][0] + HATSmt;
1014 double HATSTv = HATSk[2][0] + HATSvt;
1015 double HATSTe = HATSk[3][0] + HATSet;
1016 double HATSTp = HATSk[4][0] + HATSpt;
1017 double HATSTi = HATSk[5][0] + HATSit;
1018 double HATSTs = HATSk[6][0] + HATSst;
1019
1020 for (int ii = 1; ii < 9; ii++) {
1021 HATSTu += 2 * HATSk[0][ii];
1022 HATSTm += 2 * HATSk[1][ii];
1023 HATSTv += 2 * HATSk[2][ii];
1024 HATSTe += 2 * HATSk[3][ii];
1025 HATSTp += 2 * HATSk[4][ii];
1026 HATSTi += 2 * HATSk[5][ii];
1027 HATSTs += 2 * HATSk[6][ii];
1028 }
1029
1030 double HTu = Hk[0][0] + H0ut;
1031 double HTm = Hk[1][0] + H0mt;
1032 double HTv = Hk[2][0] + H0vt;
1033 double HTe = Hk[3][0] + H0et;
1034 double HTp = Hk[4][0] + H0pt;
1035 double HTi = Hk[5][0] + H0it;
1036 double HTs = Hk[6][0] + H0st;
1037
1038 for (int ii = 1; ii < 9; ii++) {
1039 HTu += 2.0 * Hk[0][ii];
1040 HTm += 2.0 * Hk[1][ii];
1041 HTv += 2.0 * Hk[2][ii];
1042 HTe += 2.0 * Hk[3][ii];
1043 HTp += 2.0 * Hk[4][ii];
1044 HTi += 2.0 * Hk[5][ii];
1045 HTs += 2.0 * Hk[6][ii];
1046 }
1047
1048 // 2*(Rk[1]+Rk[2]+Rk[3]+Rk[4]+Rk[5]+Rk[6]+Rk[7]+Rk[8]) // this is not true
1049 // this is on all the element
1050
1051 double RTu = 0.0;
1052 double RTm = 0.0;
1053 double RTv = 0.0;
1054 double RTe = 0.0;
1055 double RTp = 0.0;
1056 double RTi = 0.0;
1057 double RTs = 0.0;
1058
1059 for (int ii = 0; ii < 8; ii++) {
1060 RTu += 2.0 * Rk[0][ii];
1061 RTm += 2.0 * Rk[1][ii];
1062 RTv += 2.0 * Rk[2][ii];
1063 RTe += 2.0 * Rk[3][ii];
1064 RTp += 2.0 * Rk[4][ii];
1065 RTi += 2.0 * Rk[5][ii];
1066 RTs += 2.0 * Rk[6][ii];
1067 }
1068
1069 double RTMu = getMax(Rp[0]);
1070 double RTMm = getMax(Rp[1]);
1071 double RTMv = getMax(Rp[2]);
1072 double RTMe = getMax(Rp[3]);
1073 double RTMp = getMax(Rp[4]);
1074 double RTMi = getMax(Rp[5]);
1075 double RTMs = getMax(Rp[6]);
1076
1077 // create the output vector...
1078 for (int i = 0; i < 9; i++) {
1079 res[i + 4] = roundn(Hk[0][i], 3);
1080 res[i + 14] = roundn(HATSk[0][i], 3);
1081 res[i + 24] = roundn(Hk[1][i], 3);
1082 res[i + 34] = roundn(HATSk[1][i], 3);
1083 res[i + 44] = roundn(Hk[2][i], 3);
1084 res[i + 54] = roundn(HATSk[2][i], 3);
1085 res[i + 64] = roundn(Hk[3][i], 3);
1086 res[i + 74] = roundn(HATSk[3][i], 3);
1087 res[i + 84] = roundn(Hk[4][i], 3);
1088 res[i + 94] = roundn(HATSk[4][i], 3);
1089 res[i + 104] = roundn(Hk[5][i], 3);
1090 res[i + 114] = roundn(HATSk[5][i], 3);
1091 res[i + 124] = roundn(Hk[6][i], 3);
1092 res[i + 134] = roundn(HATSk[6][i], 3);
1093 }
1094
1095 res[13] = roundn(HTu, 3);
1096 res[23] = roundn(HATSTu, 3);
1097 res[33] = roundn(HTm, 3);
1098 res[43] = roundn(HATSTm, 3);
1099 res[53] = roundn(HTv, 3);
1100 res[63] = roundn(HATSTv, 3);
1101 res[73] = roundn(HTe, 3);
1102 res[83] = roundn(HATSTe, 3);
1103 res[93] = roundn(HTp, 3);
1104 res[103] = roundn(HATSTp, 3);
1105 res[113] = roundn(HTi, 3);
1106 res[123] = roundn(HATSTi, 3);
1107 res[133] = roundn(HTs, 3);
1108 res[143] = roundn(HATSTs, 3);
1109
1110 res[144] = roundn(rcon, 3);
1111 res[145] = roundn(RARS, 3);
1112 res[146] = roundn(EIG(0), 3);
1113
1114 for (int i = 0; i < 8; i++) {
1115 res[i + 147] = roundn(Rk[0][i], 3);
1116 res[i + 156] = roundn(Rp[0][i], 3);
1117 res[i + 165] = roundn(Rk[1][i], 3);
1118 res[i + 174] = roundn(Rp[1][i], 3);
1119 res[i + 183] = roundn(Rk[2][i], 3);
1120 res[i + 192] = roundn(Rp[2][i], 3);
1121 res[i + 201] = roundn(Rk[3][i], 3);
1122 res[i + 210] = roundn(Rp[3][i], 3);
1123 res[i + 219] = roundn(Rk[4][i], 3);
1124 res[i + 228] = roundn(Rp[4][i], 3);
1125 res[i + 237] = roundn(Rk[5][i], 3);
1126 res[i + 246] = roundn(Rp[5][i], 3);
1127 res[i + 255] = roundn(Rk[6][i], 3);
1128 res[i + 264] = roundn(Rp[6][i], 3);
1129 }
1130
1131 res[155] = roundn(RTu + R0ut, 3);
1132 res[164] = roundn(RTMu, 3);
1133 res[173] = roundn(RTm + R0mt, 3);
1134 res[182] = roundn(RTMm, 3);
1135 res[191] = roundn(RTv + R0vt, 3);
1136 res[200] = roundn(RTMv, 3);
1137 res[209] = roundn(RTe + R0et, 3);
1138 res[218] = roundn(RTMe, 3);
1139 res[227] = roundn(RTp + R0pt, 3);
1140 res[236] = roundn(RTMp, 3);
1141 res[245] = roundn(RTi + R0it, 3);
1142 res[254] = roundn(RTMi, 3);
1143 res[263] = roundn(RTs + R0st, 3);
1144 res[272] = roundn(RTMs, 3);
1145 }
1146
1147 /*
1148 std::vector<std::string>
1149 GETAWAYNAMES={"ITH","ISH","HIC","HGM","H0u","H1u","H2u","H3u","H4u","H5u","H6u","H7u","H8u","HTu",
1150 "HATS0u","HATS1u","HATS2u","HATS3u","HATS4u","HATS5u","HATS6u","HATS7u","HATS8u","HATSu","H0m","H1m","H2m","H3m","H4m","H5m",
1151 "H6m","H7m","H8m","HTm","HATS0m","HATS1m","HATS2m","HATS3m","HATS4m","HATS5m","HATS6m","HATS7m","HATS8m","HATSm","H0v","H1v",
1152 "H2v","H3v","H4v","H5v","H6v","H7v","H8v","HTv","HATS0v","HATS1v","HATS2v","HATS3v","HATS4v","HATS5v","HATS6v","HATS7v","HATS8v",
1153 "HATSv","H0e","H1e","H2e","H3e","H4e","H5e","H6e","H7e","H8e","HTe","HATS0e","HATS1e","HATS2e","HATS3e","HATS4e","HATS5e","HATS6e",
1154 "HATS7e","HATS8e","HATSe","H0p","H1p","H2p","H3p","H4p","H5p","H6p","H7p","H8p","HTp","HATS0p","HATS1p","HATS2p","HATS3p","HATS4p",
1155 "HATS5p","HATS6p","HATS7p","HATS8p","HATSp","H0i","H1i","H2i","H3i","H4i","H5i","H6i","H7i","H8i","HTi","HATS0i","HATS1i","HATS2i",
1156 "HATS3i","HATS4i","HATS5i","HATS6i","HATS7i","HATS8i","HATSi","H0s","H1s","H2s","H3s","H4s","H5s","H6s","H7s","H8s","HTs","HATS0s",
1157 "HATS1s","HATS2s","HATS3s","HATS4s","HATS5s","HATS6s","HATS7s","HATS8s","HATSs","RCON","RARS","REIG","R1u","R2u","R3u","R4u","R5u",
1158 "R6u","R7u","R8u","RTu","R1u+","R2u+","R3u+","R4u+","R5u+","R6u+","R7u+","R8u+","RTu+","R1m","R2m","R3m","R4m","R5m","R6m","R7m",
1159 "R8m","RTm","R1m+","R2m+","R3m+","R4m+","R5m+","R6m+","R7m+","R8m+","RTm+","R1v","R2v","R3v","R4v","R5v","R6v","R7v","R8v","RTv",
1160 "R1v+","R2v+","R3v+","R4v+","R5v+","R6v+","R7v+","R8v+","RTv+","R1e","R2e","R3e","R4e","R5e","R6e","R7e","R8e","RTe","R1e+","R2e+",
1161 "R3e+","R4e+","R5e+","R6e+","R7e+","R8e+","RTe+","R1p","R2p","R3p","R4p","R5p","R6p","R7p","R8p","RTp","R1p+","R2p+","R3p+","R4p+",
1162 "R5p+","R6p+","R7p+","R8p+","RTp+","R1i","R2i","R3i","R4i","R5i","R6i","R7i","R8i","RTi","R1i+","R2i+","R3i+","R4i+","R5i+","R6i+",
1163 "R7i+","R8i+","RTi+","R1s","R2s","R3s","R4s","R5s","R6s","R7s","R8s","RTs","R1s+","R2s+","R3s+","R4s+","R5s+","R6s+","R7s+","R8s+","RTs+"};
1164 */
1165
GetGETAWAYone(double * dist3D,double * AdjMat,std::vector<double> Vpoints,const ROMol & mol,const Conformer & conf,std::vector<int> Heavylist,std::vector<double> & res,unsigned int precision,const std::string & customAtomPropName)1166 void GetGETAWAYone(double* dist3D, double* AdjMat, std::vector<double> Vpoints,
1167 const ROMol& mol, const Conformer& conf,
1168 std::vector<int> Heavylist, std::vector<double>& res,
1169 unsigned int precision,
1170 const std::string& customAtomPropName) {
1171 PRECONDITION(dist3D != nullptr, "no distance matrix");
1172 PRECONDITION(AdjMat != nullptr, "no adjacency matrix");
1173 int numAtoms = conf.getNumAtoms();
1174
1175 Map<MatrixXd> ADJ(AdjMat, numAtoms, numAtoms);
1176
1177 Map<MatrixXd> DM(dist3D, numAtoms, numAtoms);
1178
1179 double* vpoints = &Vpoints[0];
1180
1181 Map<MatrixXd> matorigin(vpoints, 3, numAtoms);
1182
1183 MatrixXd MatOrigin = matorigin.transpose();
1184
1185 MatrixXd Xmean = GetCenterMatrix(MatOrigin);
1186
1187 MatrixXd H = GetHmatrix(Xmean);
1188
1189 MatrixXd R = GetRmatrix(H, DM, numAtoms);
1190
1191 getGETAWAYDescCustom(H, R, ADJ, numAtoms, std::move(Heavylist), mol, res,
1192 precision, customAtomPropName);
1193 }
1194
GetGETAWAY(double * dist3D,double * AdjMat,std::vector<double> Vpoints,const ROMol & mol,const Conformer & conf,std::vector<int> Heavylist,std::vector<double> & res,unsigned int precision)1195 void GetGETAWAY(double* dist3D, double* AdjMat, std::vector<double> Vpoints,
1196 const ROMol& mol, const Conformer& conf,
1197 std::vector<int> Heavylist, std::vector<double>& res,
1198 unsigned int precision) {
1199 PRECONDITION(dist3D != nullptr, "no distance matrix");
1200 PRECONDITION(AdjMat != nullptr, "no adjacency matrix");
1201
1202 int numAtoms = conf.getNumAtoms();
1203
1204 Map<MatrixXd> ADJ(AdjMat, numAtoms, numAtoms);
1205
1206 Map<MatrixXd> DM(dist3D, numAtoms, numAtoms);
1207
1208 double* vpoints = &Vpoints[0];
1209
1210 Map<MatrixXd> matorigin(vpoints, 3, numAtoms);
1211
1212 MatrixXd MatOrigin = matorigin.transpose();
1213
1214 MatrixXd Xmean = GetCenterMatrix(MatOrigin);
1215
1216 MatrixXd H = GetHmatrix(Xmean);
1217
1218 MatrixXd R = GetRmatrix(H, DM, numAtoms);
1219
1220 getGETAWAYDesc(H, R, ADJ, numAtoms, std::move(Heavylist), mol, res,
1221 precision);
1222 }
1223
1224 } // end of anonymous namespace
1225
GETAWAY(const ROMol & mol,std::vector<double> & res,int confId,unsigned int precision,const std::string & customAtomPropName)1226 void GETAWAY(const ROMol& mol, std::vector<double>& res, int confId,
1227 unsigned int precision, const std::string& customAtomPropName) {
1228 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers")
1229
1230 int numAtoms = mol.getNumAtoms();
1231
1232 const Conformer& conf = mol.getConformer(confId);
1233
1234 std::vector<double> Vpoints(3 * numAtoms);
1235
1236 for (int i = 0; i < numAtoms; ++i) {
1237 Vpoints[3 * i] = conf.getAtomPos(i).x;
1238 Vpoints[3 * i + 1] = conf.getAtomPos(i).y;
1239 Vpoints[3 * i + 2] = conf.getAtomPos(i).z;
1240 }
1241
1242 std::vector<int> Heavylist = GetHeavyList(mol);
1243 // should be the same as the List size upper!
1244 // int nHeavyAt= mol.getNumHeavyAtoms();
1245
1246 double* dist3D = MolOps::get3DDistanceMat(mol, confId);
1247
1248 double* AdjMat = MolOps::getAdjacencyMatrix(
1249 mol, false, 0, false,
1250 nullptr); // false to have only the 1,0 matrix unweighted
1251
1252 if (!customAtomPropName.empty()) {
1253 res.clear();
1254 res.resize(45);
1255 GetGETAWAYone(dist3D, AdjMat, Vpoints, mol, conf, Heavylist, res, precision,
1256 customAtomPropName);
1257 } else {
1258 res.clear();
1259 res.resize(273);
1260
1261 GetGETAWAY(dist3D, AdjMat, Vpoints, mol, conf, Heavylist, res, precision);
1262 }
1263 }
1264
1265 } // namespace Descriptors
1266 } // namespace RDKit
1267