1 /********************************************************************** 2 obutil.cpp - Various utility methods. 3 4 Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc. 5 Some portions Copyright (C) 2001-2006 by Geoffrey R. Hutchison 6 7 This file is part of the Open Babel project. 8 For more information, see <http://openbabel.org/> 9 10 This program is free software; you can redistribute it and/or modify 11 it under the terms of the GNU General Public License as published by 12 the Free Software Foundation version 2 of the License. 13 14 This program is distributed in the hope that it will be useful, 15 but WITHOUT ANY WARRANTY; without even the implied warranty of 16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 17 GNU General Public License for more details. 18 ***********************************************************************/ 19 20 #include <openbabel/babelconfig.h> 21 #include <openbabel/math/matrix3x3.h> 22 #include <openbabel/mol.h> 23 #include <openbabel/atom.h> 24 #include <openbabel/obiter.h> 25 #include <openbabel/obutil.h> 26 #include <openbabel/internalcoord.h> 27 28 #include <cstring> 29 30 #ifdef HAVE_CONIO_H 31 #include <conio.h> 32 #endif 33 34 using namespace std; 35 namespace OpenBabel 36 { 37 38 /*! \class OBStopwatch obutil.h <openbabel/obutil.h> 39 \brief Stopwatch class used for timing length of execution 40 41 The OBStopwatch class makes timing the execution of blocks of 42 code to microsecond accuracy very simple. The class effectively 43 has two functions, Start() and Elapsed(). The usage of the 44 OBStopwatch class is demonstrated by the following code: 45 \code 46 OBStopwatch sw; 47 sw.Start(); 48 //insert code here 49 cout << "Elapsed time = " << sw.Elapsed() << endl; 50 \endcode 51 */ 52 53 //! Deprecated: use the OBMessageHandler class instead 54 //! \deprecated Throw an error through the OpenBabel::OBMessageHandler class ThrowError(char * str)55 void ThrowError(char *str) 56 { 57 obErrorLog.ThrowError("", str, obInfo); 58 } 59 60 //! Deprecated: use the OBMessageHandler class instead 61 //! \deprecated Throw an error through the OpenBabel::OBMessageHandler class ThrowError(std::string & str)62 void ThrowError(std::string &str) 63 { 64 obErrorLog.ThrowError("", str, obInfo); 65 } 66 67 // returns True if a < b, False otherwise. OBCompareInt(const int & a,const int & b)68 bool OBCompareInt(const int &a,const int &b) 69 { 70 return(a<b); 71 } 72 73 // Comparison function (for sorting unsigned ints) returns a < b OBCompareUnsigned(const unsigned int & a,const unsigned int & b)74 bool OBCompareUnsigned(const unsigned int &a,const unsigned int &b) 75 { 76 return(a<b); 77 } 78 79 //! Comparison for doubles: returns fabs(a - b) < epsilon IsNear(const double & a,const double & b,const double epsilon)80 bool IsNear(const double &a, const double &b, const double epsilon) 81 { 82 return (fabs(a - b) < epsilon); 83 } 84 85 //! Comparison for doubles: returns fabs(a) < epsilon IsNearZero(const double & a,const double epsilon)86 bool IsNearZero(const double &a, const double epsilon) 87 { 88 return (fabs(a) < epsilon); 89 } 90 91 //! Comparison for nan (not a number) IsNan(const double & a)92 bool IsNan(const double &a) 93 { 94 return ((a) != (a)); 95 } 96 97 //! Tests whether its argument can be squared without triggering an overflow or 98 //! underflow. CanBeSquared(const double & a)99 bool CanBeSquared(const double &a) 100 { 101 if( a == 0 ) return true; 102 const double max_squarable_double = 1e150; 103 const double min_squarable_double = 1e-150; 104 double abs_a = fabs(a); 105 return(abs_a < max_squarable_double && abs_a > min_squarable_double); 106 } 107 108 //! Utility function: replace the last extension in string &src with new extension char *ext. NewExtension(string & src,char * ext)109 string NewExtension(string &src,char *ext) 110 { 111 string::size_type pos = (unsigned int)src.find_last_of("."); 112 string dst; 113 114 if (pos != string::npos) 115 dst = src.substr(0,pos+1); 116 else 117 { 118 dst = src; 119 dst += "."; 120 } 121 122 dst += ext; 123 return(dst); 124 } 125 126 //! \return the geometric centroid to an array of coordinates in double* format 127 //! and center the coordinates to the origin. Operates on the first "size" 128 //! coordinates in the array. center_coords(double * c,unsigned int size)129 vector3 center_coords(double *c, unsigned int size) 130 { 131 if (size == 0) 132 { 133 return(VZero); 134 } 135 unsigned int i; 136 double x=0.0, y=0.0, z=0.0; 137 for (i = 0;i < size;++i) 138 { 139 x += c[i*3]; 140 y += c[i*3+1]; 141 z += c[i*3+2]; 142 } 143 x /= (double) size; 144 y /= (double) size; 145 z /= (double) size; 146 for (i = 0;i < size;++i) 147 { 148 c[i*3] -= x; 149 c[i*3+1] -= y; 150 c[i*3+2] -= z; 151 } 152 vector3 v(x,y,z); 153 return(v); 154 } 155 156 //! Rotates the coordinate set *c by the transformation matrix m[3][3] 157 //! Operates on the first "size" coordinates in the array. rotate_coords(double * c,double m[3][3],unsigned int size)158 void rotate_coords(double *c,double m[3][3],unsigned int size) 159 { 160 double x,y,z; 161 for (unsigned int i = 0;i < size;++i) 162 { 163 x = c[i*3]*m[0][0] + c[i*3+1]*m[0][1] + c[i*3+2]*m[0][2]; 164 y = c[i*3]*m[1][0] + c[i*3+1]*m[1][1] + c[i*3+2]*m[1][2]; 165 z = c[i*3]*m[2][0] + c[i*3+1]*m[2][1] + c[i*3+2]*m[2][2]; 166 c[i*3] = x; 167 c[i*3+1] = y; 168 c[i*3+2] = z; 169 } 170 } 171 172 //! Calculate the RMS deviation between the first N coordinates of *r and *f calc_rms(double * r,double * f,unsigned int N)173 double calc_rms(double *r,double *f, unsigned int N) 174 { 175 if (N == 0) 176 return 0.0; // no RMS deviation between two empty sets 177 178 double d2=0.0; 179 for (unsigned int i = 0;i < N;++i) 180 { 181 d2 += SQUARE(r[i*3] - f[i*3]) + 182 SQUARE(r[i*3+1] - f[i*3+1]) + 183 SQUARE(r[i*3+2] - f[i*3+2]); 184 } 185 186 d2 /= (double) N; 187 return(sqrt(d2)); 188 } 189 190 //! Rotate the coordinates of 'atoms' 191 //! such that tor == ang - atoms in 'tor' should be ordered such 192 //! that the 3rd atom is the pivot around which atoms rotate SetRotorToAngle(double * c,vector<int> & tor,double ang,vector<int> & atoms)193 void SetRotorToAngle(double *c,vector<int> &tor,double ang,vector<int> &atoms) 194 { 195 double v1x,v1y,v1z,v2x,v2y,v2z,v3x,v3y,v3z; 196 double c1x,c1y,c1z,c2x,c2y,c2z,c3x,c3y,c3z; 197 double c1mag,c2mag,radang,costheta,m[9]; 198 double x,y,z,mag,rotang,sn,cs,t,tx,ty,tz; 199 200 // 201 //calculate the torsion angle 202 // 203 v1x = c[tor[0]] - c[tor[1]]; 204 v2x = c[tor[1]] - c[tor[2]]; 205 v1y = c[tor[0]+1] - c[tor[1]+1]; 206 v2y = c[tor[1]+1] - c[tor[2]+1]; 207 v1z = c[tor[0]+2] - c[tor[1]+2]; 208 v2z = c[tor[1]+2] - c[tor[2]+2]; 209 v3x = c[tor[2]] - c[tor[3]]; 210 v3y = c[tor[2]+1] - c[tor[3]+1]; 211 v3z = c[tor[2]+2] - c[tor[3]+2]; 212 213 c1x = v1y*v2z - v1z*v2y; 214 c2x = v2y*v3z - v2z*v3y; 215 c1y = -v1x*v2z + v1z*v2x; 216 c2y = -v2x*v3z + v2z*v3x; 217 c1z = v1x*v2y - v1y*v2x; 218 c2z = v2x*v3y - v2y*v3x; 219 c3x = c1y*c2z - c1z*c2y; 220 c3y = -c1x*c2z + c1z*c2x; 221 c3z = c1x*c2y - c1y*c2x; 222 223 c1mag = SQUARE(c1x)+SQUARE(c1y)+SQUARE(c1z); 224 c2mag = SQUARE(c2x)+SQUARE(c2y)+SQUARE(c2z); 225 if (c1mag*c2mag < 0.01) 226 costheta = 1.0; //avoid div by zero error 227 else 228 costheta = (c1x*c2x + c1y*c2y + c1z*c2z)/(sqrt(c1mag*c2mag)); 229 230 if (costheta < -0.999999) 231 costheta = -0.999999; 232 if (costheta > 0.999999) 233 costheta = 0.999999; 234 235 if ((v2x*c3x + v2y*c3y + v2z*c3z) > 0.0) 236 radang = -acos(costheta); 237 else 238 radang = acos(costheta); 239 240 // 241 // now we have the torsion angle (radang) - set up the rot matrix 242 // 243 244 //find the difference between current and requested 245 rotang = ang - radang; 246 247 sn = sin(rotang); 248 cs = cos(rotang); 249 t = 1 - cs; 250 //normalize the rotation vector 251 mag = sqrt(SQUARE(v2x)+SQUARE(v2y)+SQUARE(v2z)); 252 x = v2x/mag; 253 y = v2y/mag; 254 z = v2z/mag; 255 256 //set up the rotation matrix 257 m[0]= t*x*x + cs; 258 m[1] = t*x*y + sn*z; 259 m[2] = t*x*z - sn*y; 260 m[3] = t*x*y - sn*z; 261 m[4] = t*y*y + cs; 262 m[5] = t*y*z + sn*x; 263 m[6] = t*x*z + sn*y; 264 m[7] = t*y*z - sn*x; 265 m[8] = t*z*z + cs; 266 267 // 268 //now the matrix is set - time to rotate the atoms 269 // 270 tx = c[tor[1]]; 271 ty = c[tor[1]+1]; 272 tz = c[tor[1]+2]; 273 vector<int>::iterator i; 274 int j; 275 for (i = atoms.begin();i != atoms.end();++i) 276 { 277 j = *i; 278 c[j] -= tx; 279 c[j+1] -= ty; 280 c[j+2]-= tz; 281 x = c[j]*m[0] + c[j+1]*m[1] + c[j+2]*m[2]; 282 y = c[j]*m[3] + c[j+1]*m[4] + c[j+2]*m[5]; 283 z = c[j]*m[6] + c[j+1]*m[7] + c[j+2]*m[8]; 284 c[j] = x; 285 c[j+1] = y; 286 c[j+2] = z; 287 c[j] += tx; 288 c[j+1] += ty; 289 c[j+2] += tz; 290 } 291 } 292 293 //! Safely open the supplied filename and return an ifstream, throwing an error 294 //! to the default OBMessageHandler error log if it fails. SafeOpen(std::ifstream & fs,const char * filename)295 bool SafeOpen(std::ifstream &fs, const char *filename) 296 { 297 #ifdef WIN32 298 string s(filename); 299 if (s.find(".bin") != string::npos) 300 fs.open(filename,ios::binary); 301 else 302 #endif 303 304 fs.open(filename); 305 306 if (!fs) 307 { 308 string error = "Unable to open file \'"; 309 error += filename; 310 error += "\' in read mode"; 311 obErrorLog.ThrowError(__FUNCTION__, error, obError); 312 return(false); 313 } 314 315 return(true); 316 } 317 318 319 //! Safely open the supplied filename and return an ofstream, throwing an error 320 //! to the default OBMessageHandler error log if it fails. SafeOpen(std::ofstream & fs,const char * filename)321 bool SafeOpen(std::ofstream &fs, const char *filename) 322 { 323 #ifdef WIN32 324 string s(filename); 325 if (s.find(".bin") != string::npos) 326 fs.open(filename,ios::binary); 327 else 328 #endif 329 330 fs.open(filename); 331 332 if (!fs) 333 { 334 string error = "Unable to open file \'"; 335 error += filename; 336 error += "\' in write mode"; 337 obErrorLog.ThrowError(__FUNCTION__, error, obError); 338 return(false); 339 } 340 341 return(true); 342 } 343 344 //! Safely open the supplied filename and return an ifstream, throwing an error 345 //! to the default OBMessageHandler error log if it fails. SafeOpen(std::ifstream & fs,const string & filename)346 bool SafeOpen(std::ifstream &fs, const string &filename) 347 { 348 return(SafeOpen(fs, filename.c_str())); 349 } 350 351 //! Safely open the supplied filename and return an ofstream, throwing an error 352 //! to the default OBMessageHandler error log if it fails. SafeOpen(std::ofstream & fs,const string & filename)353 bool SafeOpen(std::ofstream &fs, const string &filename) 354 { 355 return(SafeOpen(fs, filename.c_str())); 356 } 357 358 //! Shift the supplied string to uppercase ToUpper(std::string & s)359 void ToUpper(std::string &s) 360 { 361 if (s.empty()) 362 return; 363 unsigned int i; 364 for (i = 0;i < s.size();++i) 365 if (isalpha(s[i]) && !isdigit(s[i])) 366 s[i] = toupper(s[i]); 367 } 368 369 //! Shift the supplied char* to uppercase ToUpper(char * cptr)370 void ToUpper(char *cptr) 371 { 372 char *c; 373 for (c = cptr;*c != '\0';++c) 374 if (isalpha(*c) && !isdigit(*c)) 375 *c = toupper(*c); 376 } 377 378 //! Shift the supplied string to lowercase ToLower(std::string & s)379 void ToLower(std::string &s) 380 { 381 if (s.empty()) 382 return; 383 unsigned int i; 384 for (i = 0;i < s.size();++i) 385 if (isalpha(s[i]) && !isdigit(s[i])) 386 s[i] = tolower(s[i]); 387 } 388 389 //! Shift the supplied char* to lowercase ToLower(char * cptr)390 void ToLower(char *cptr) 391 { 392 char *c; 393 for (c = cptr;*c != '\0';++c) 394 if (isalpha(*c) && !isdigit(*c)) 395 *c = tolower(*c); 396 } 397 398 //! Shift the supplied string: lowercase to upper, and upper to lower 399 //! \param s - The string to switch case 400 //! \param start - The position to start inverting case InvertCase(std::string & s,unsigned int start)401 void InvertCase(std::string &s, unsigned int start) 402 { 403 if (start >= s.size()) 404 return; 405 unsigned int i; 406 for (i = start; i < s.size();++i) 407 if (isalpha(s[i]) && !isdigit(s[i])) { 408 if (isupper(s[i])) s[i] = tolower(s[i]); 409 else s[i] = toupper(s[i]); 410 } 411 } 412 413 //! Shift the supplied char*: lowercase to upper, and upper to lower InvertCase(char * cptr)414 void InvertCase(char *cptr) 415 { 416 char *c; 417 for (c = cptr;*c != '\0';++c) 418 if (isalpha(*c) && !isdigit(*c)) { 419 if (isupper(*c)) *c = tolower(*c); 420 else *c = toupper(*c); 421 } 422 } 423 424 //! "Clean" the supplied atom type, shifting the first character to uppercase, 425 //! the second character (if it's a letter) to lowercase, and terminating with a NULL 426 //! to strip off any trailing characters CleanAtomType(char * id)427 void CleanAtomType(char *id) 428 { 429 id[0] = toupper(id[0]); 430 if (isalpha(id[1]) == 0) 431 id[1] = '\0'; 432 else 433 { 434 id[1] = tolower(id[1]); 435 id[2] = '\0'; 436 } 437 } 438 439 //! Transform the supplied vector<OBInternalCoord*> into cartesian and update 440 //! the OBMol accordingly. The size of supplied internal coordinate vector 441 //! has to be the same as the number of atoms in molecule (+ NULL in the 442 //! beginning). 443 //! Implements <a href="http://qsar.sourceforge.net/dicts/blue-obelisk/index.xhtml#zmatrixCoordinatesIntoCartesianCoordinates">blue-obelisk:zmatrixCoordinatesIntoCartesianCoordinates</a> InternalToCartesian(std::vector<OBInternalCoord * > & vic,OBMol & mol)444 void InternalToCartesian(std::vector<OBInternalCoord*> &vic,OBMol &mol) 445 { 446 vector3 n,nn,v1,v2,v3,avec,bvec,cvec; 447 double dst = 0.0, ang = 0.0, tor = 0.0; 448 OBAtom *atom; 449 vector<OBAtom*>::iterator i; 450 unsigned int index; 451 452 if (vic.empty()) 453 return; 454 455 if (vic[0] != nullptr) { 456 std::vector<OBInternalCoord*>::iterator it = vic.begin(); 457 vic.insert(it, nullptr); 458 } 459 460 if (vic.size() != mol.NumAtoms() + 1) { 461 string error = "Number of internal coordinates is not the same as"; 462 error += " the number of atoms in molecule"; 463 obErrorLog.ThrowError(__FUNCTION__, error, obError); 464 return; 465 } 466 467 obErrorLog.ThrowError(__FUNCTION__, 468 "Ran OpenBabel::InternalToCartesian", obAuditMsg); 469 470 for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i)) 471 { 472 index = atom->GetIdx(); 473 474 // make sure we always have valid pointers 475 if (index >= vic.size() || !vic[index]) 476 return; 477 478 if (vic[index]->_a) // make sure we have a valid ptr 479 { 480 avec = vic[index]->_a->GetVector(); 481 dst = vic[index]->_dst; 482 } 483 else 484 { 485 // atom 1 486 atom->SetVector(0.0, 0.0, 0.0); 487 continue; 488 } 489 490 if (vic[index]->_b) 491 { 492 bvec = vic[index]->_b->GetVector(); 493 ang = vic[index]->_ang * DEG_TO_RAD; 494 } 495 else 496 { 497 // atom 2 498 atom->SetVector(dst, 0.0, 0.0); 499 continue; 500 } 501 502 if (vic[index]->_c) 503 { 504 cvec = vic[index]->_c->GetVector(); 505 tor = vic[index]->_tor * DEG_TO_RAD; 506 } 507 else 508 { 509 // atom 3 510 cvec = VY; 511 tor = 90. * DEG_TO_RAD; 512 } 513 514 v1 = avec - bvec; 515 v2 = avec - cvec; 516 n = cross(v1,v2); 517 nn = cross(v1,n); 518 n.normalize(); 519 nn.normalize(); 520 521 n *= -sin(tor); 522 nn *= cos(tor); 523 v3 = n + nn; 524 v3.normalize(); 525 v3 *= dst * sin(ang); 526 v1.normalize(); 527 v1 *= dst * cos(ang); 528 v2 = avec + v3 - v1; 529 530 atom->SetVector(v2); 531 } 532 533 // Delete dummy atoms 534 vector<OBAtom*> for_deletion; 535 FOR_ATOMS_OF_MOL(a, mol) 536 if (a->GetAtomicNum() == 0) 537 for_deletion.push_back(&(*a)); 538 for(vector<OBAtom*>::iterator a_it=for_deletion.begin(); a_it!=for_deletion.end(); ++a_it) 539 mol.DeleteAtom(*a_it); 540 541 } 542 543 //! Use the supplied OBMol and its Cartesian coordinates to generate 544 //! a set of internal (z-matrix) coordinates as supplied in the 545 //! vector<OBInternalCoord*> argument. 546 //! Implements <a href="http://qsar.sourceforge.net/dicts/blue-obelisk/index.xhtml#cartesianCoordinatesIntoZmatrixCoordinates">blue-obelisk:cartesianCoordinatesIntoZmatrixCoordinates</a>. 547 //! \todo Consider lengths, angles, and torsions for periodic systems CartesianToInternal(std::vector<OBInternalCoord * > & vic,OBMol & mol)548 void CartesianToInternal(std::vector<OBInternalCoord*> &vic,OBMol &mol) 549 { 550 double r,sum; 551 OBAtom *atom,*nbr,*ref; 552 vector<OBAtom*>::iterator i,j,m; 553 554 obErrorLog.ThrowError(__FUNCTION__, 555 "Ran OpenBabel::CartesianToInternal", obAuditMsg); 556 557 //set reference atoms 558 for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i)) 559 { 560 if (atom->GetIdx() == 1) 561 continue; 562 else if (atom->GetIdx() == 2) 563 { 564 vic[atom->GetIdx()]->_a = mol.GetAtom(1); 565 continue; 566 } 567 else if (atom->GetIdx() == 3) 568 { 569 if( (atom->GetVector()-mol.GetAtom(2)->GetVector()).length_2() 570 <(atom->GetVector()-mol.GetAtom(1)->GetVector()).length_2()) 571 { 572 vic[atom->GetIdx()]->_a = mol.GetAtom(2); 573 vic[atom->GetIdx()]->_b = mol.GetAtom(1); 574 } 575 else 576 { 577 vic[atom->GetIdx()]->_a = mol.GetAtom(1); 578 vic[atom->GetIdx()]->_b = mol.GetAtom(2); 579 } 580 continue; 581 } 582 sum=1.0E10; 583 ref = mol.GetAtom(1); 584 for(nbr = mol.BeginAtom(j);nbr && (i != j);nbr = mol.NextAtom(j)) 585 { 586 r = (atom->GetVector()-nbr->GetVector()).length_2(); 587 if((r < sum) && (vic[nbr->GetIdx()]->_a != nbr) && 588 (vic[nbr->GetIdx()]->_b != nbr)) 589 { 590 sum = r; 591 ref = nbr; 592 } 593 } 594 595 vic[atom->GetIdx()]->_a = ref; 596 if (ref->GetIdx() >= 3) 597 { 598 vic[atom->GetIdx()]->_b = vic[ref->GetIdx()]->_a; 599 vic[atom->GetIdx()]->_c = vic[ref->GetIdx()]->_b; 600 } 601 else 602 { 603 if(ref->GetIdx()== 1) 604 { 605 vic[atom->GetIdx()]->_b = mol.GetAtom(2); 606 vic[atom->GetIdx()]->_c = mol.GetAtom(3); 607 } 608 else 609 {//ref->GetIdx()== 2 610 vic[atom->GetIdx()]->_b = mol.GetAtom(1); 611 vic[atom->GetIdx()]->_c = mol.GetAtom(3); 612 } 613 } 614 } 615 616 //fill in geometries 617 unsigned int k; 618 vector3 v1,v2; 619 OBAtom *a,*b,*c; 620 for (k = 2;k <= mol.NumAtoms();++k) 621 { 622 atom = mol.GetAtom(k); 623 a = vic[k]->_a; 624 b = vic[k]->_b; 625 c = vic[k]->_c; 626 v1 = atom->GetVector() - a->GetVector(); 627 vic[k]->_dst = v1.length(); 628 if (k == 2) 629 continue; 630 631 v2 = b->GetVector() - a->GetVector(); 632 vic[k]->_ang = vectorAngle(v1,v2); 633 if (k == 3) 634 continue; 635 636 vic[k]->_tor = CalcTorsionAngle(atom->GetVector(), 637 a->GetVector(), 638 b->GetVector(), 639 c->GetVector()); 640 } 641 642 //check for linear geometries and try to correct if possible 643 bool done; 644 double ang; 645 for (k = 2;k <= mol.NumAtoms();++k) 646 { 647 ang = fabs(vic[k]->_ang); 648 if (ang > 5.0 && ang < 175.0) 649 continue; 650 atom = mol.GetAtom(k); 651 done = false; 652 for (a = mol.BeginAtom(i);a && a->GetIdx() < k && !done;a = mol.NextAtom(i)) 653 for (b=mol.BeginAtom(j);b && b->GetIdx()<a->GetIdx() && !done;b = mol.NextAtom(j)) 654 { 655 v1 = atom->GetVector() - a->GetVector(); 656 v2 = b->GetVector() - a->GetVector(); 657 ang = fabs(vectorAngle(v1,v2)); 658 if (ang < 5.0 || ang > 175.0) 659 continue; 660 661 // Also check length considerations -- don't bother if the length > 10.0 Angstroms 662 if (v1.length_2() > 99.999) 663 continue; 664 665 for (c = mol.BeginAtom(m);c && c->GetIdx() < atom->GetIdx();c = mol.NextAtom(m)) 666 if (c != atom && c != a && c != b) 667 break; 668 if (!c) 669 continue; 670 671 vic[k]->_a = a; 672 vic[k]->_b = b; 673 vic[k]->_c = c; 674 vic[k]->_dst = v1.length(); 675 vic[k]->_ang = vectorAngle(v1,v2); 676 vic[k]->_tor = CalcTorsionAngle(atom->GetVector(), 677 a->GetVector(), 678 b->GetVector(), 679 c->GetVector()); 680 if (!isfinite(vic[k]->_tor)) 681 vic[k]->_tor = 180.0; 682 done = true; 683 } 684 } 685 } 686 qtrfit(double * r,double * f,int size,double u[3][3])687 void qtrfit (double *r,double *f,int size, double u[3][3]) 688 { 689 int i; 690 double xxyx, xxyy, xxyz; 691 double xyyx, xyyy, xyyz; 692 double xzyx, xzyy, xzyz; 693 double d[4],q[4]; 694 double c[16],v[16]; 695 double rx,ry,rz,fx,fy,fz; 696 697 /* generate the upper triangle of the quadratic form matrix */ 698 699 xxyx = 0.0; 700 xxyy = 0.0; 701 xxyz = 0.0; 702 xyyx = 0.0; 703 xyyy = 0.0; 704 xyyz = 0.0; 705 xzyx = 0.0; 706 xzyy = 0.0; 707 xzyz = 0.0; 708 709 for (i = 0; i < size; ++i) 710 { 711 rx = r[i*3]; 712 ry = r[i*3+1]; 713 rz = r[i*3+2]; 714 fx = f[i*3]; 715 fy = f[i*3+1]; 716 fz = f[i*3+2]; 717 718 xxyx += fx * rx; 719 xxyy += fx * ry; 720 xxyz += fx * rz; 721 xyyx += fy * rx; 722 xyyy += fy * ry; 723 xyyz += fy * rz; 724 xzyx += fz * rx; 725 xzyy += fz * ry; 726 xzyz += fz * rz; 727 } 728 729 c[4*0+0] = xxyx + xyyy + xzyz; 730 731 c[4*0+1] = xzyy - xyyz; 732 c[4*1+1] = xxyx - xyyy - xzyz; 733 734 c[4*0+2] = xxyz - xzyx; 735 c[4*1+2] = xxyy + xyyx; 736 c[4*2+2] = xyyy - xzyz - xxyx; 737 738 c[4*0+3] = xyyx - xxyy; 739 c[4*1+3] = xzyx + xxyz; 740 c[4*2+3] = xyyz + xzyy; 741 c[4*3+3] = xzyz - xxyx - xyyy; 742 743 /* diagonalize c */ 744 745 matrix3x3::jacobi(4, c, d, v); 746 747 /* extract the desired quaternion */ 748 749 q[0] = v[4*0+3]; 750 q[1] = v[4*1+3]; 751 q[2] = v[4*2+3]; 752 q[3] = v[4*3+3]; 753 754 /* generate the rotation matrix */ 755 756 u[0][0] = q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3]; 757 u[1][0] = 2.0 * (q[1] * q[2] - q[0] * q[3]); 758 u[2][0] = 2.0 * (q[1] * q[3] + q[0] * q[2]); 759 760 u[0][1] = 2.0 * (q[2] * q[1] + q[0] * q[3]); 761 u[1][1] = q[0]*q[0] - q[1]*q[1] + q[2]*q[2] - q[3]*q[3]; 762 u[2][1] = 2.0 * (q[2] * q[3] - q[0] * q[1]); 763 764 u[0][2] = 2.0 * (q[3] * q[1] - q[0] * q[2]); 765 u[1][2] = 2.0 * (q[3] * q[2] + q[0] * q[1]); 766 u[2][2] = q[0]*q[0] - q[1]*q[1] - q[2]*q[2] + q[3]*q[3]; 767 } 768 769 770 771 static double Roots[4]; 772 773 #define ApproxZero 1E-7 774 #define IsZero(x) ((double)fabs(x)<ApproxZero) 775 #ifndef PI 776 #define PI 3.14159265358979323846226433 777 #endif 778 #define OneThird (1.0/3.0) 779 #define FourThirdsPI (4.0*PI/3.0) 780 #define TwoThirdsPI (2.0*PI/3.0) 781 782 /*FUNCTION */ 783 /* receives: the co-efficients for a general 784 * equation of degree one. 785 * Ax + B = 0 !! 786 */ SolveLinear(double A,double B)787 int SolveLinear(double A,double B) 788 { 789 if( IsZero(A) ) 790 return( 0 ); 791 Roots[0] = -B/A; 792 return( 1 ); 793 } 794 795 /*FUNCTION */ 796 /* receives: the co-efficients for a general 797 * linear equation of degree two. 798 * Ax^2 + Bx + C = 0 !! 799 */ SolveQuadratic(double A,double B,double C)800 int SolveQuadratic(double A,double B,double C) 801 { 802 double Descr, Temp, TwoA; 803 804 if( IsZero(A) ) 805 return( SolveLinear(B,C) ); 806 807 TwoA = A+A; 808 Temp = TwoA*C; 809 Descr = B*B - (Temp+Temp); 810 if( Descr<0.0 ) 811 return( 0 ); 812 813 if( Descr>0.0 ) 814 { 815 Descr = sqrt(Descr); 816 /* W. Press, B. Flannery, S. Teukolsky and W. Vetterling, 817 * "Quadratic and Cubic Equations", Numerical Recipes in C, 818 * Chapter 5, pp. 156-157, 1989. 819 */ 820 Temp = (B<0.0)? -0.5*(B-Descr) : -0.5*(B+Descr); 821 Roots[0] = Temp/A; 822 Roots[1] = C/Temp; 823 824 return( 2 ); 825 } 826 Roots[0] = -B/TwoA; 827 return( 1 ); 828 } 829 830 /*FUNCTION */ 831 /* task: to return the cube root of the 832 * given value taking into account 833 * that it may be negative. 834 */ CubeRoot(double X)835 double CubeRoot(double X) 836 { 837 if( X>=0.0 ) 838 { 839 return pow( X, OneThird ); 840 } 841 else 842 return -pow( -X, OneThird ); 843 } 844 SolveCubic(double A,double B,double C,double D)845 int SolveCubic(double A,double B,double C,double D) 846 { 847 double TwoA, ThreeA, BOver3A; 848 double Temp, POver3, QOver2; 849 double Desc, Rho, Psi; 850 851 852 if( IsZero(A) ) 853 { 854 return( SolveQuadratic(B,C,D) ); 855 } 856 857 TwoA = A+A; 858 ThreeA = TwoA+A; 859 BOver3A = B/ThreeA; 860 QOver2 = ((TwoA*BOver3A*BOver3A-C)*BOver3A+D)/TwoA; 861 POver3 = (C-B*BOver3A)/ThreeA; 862 863 864 Rho = POver3*POver3*POver3; 865 Desc = QOver2*QOver2 + Rho; 866 867 if( Desc<=0.0 ) 868 { 869 Rho = sqrt( -Rho ); 870 Psi = OneThird*acos(-QOver2/Rho); 871 Temp = CubeRoot( Rho ); 872 Temp = Temp+Temp; 873 874 Roots[0] = Temp*cos( Psi )-BOver3A; 875 Roots[1] = Temp*cos( Psi+TwoThirdsPI )-BOver3A; 876 Roots[2] = Temp*cos( Psi+FourThirdsPI )-BOver3A; 877 return( 3 ); 878 } 879 880 if( Desc> 0.0 ) 881 { 882 Temp = CubeRoot( -QOver2 ); 883 Roots[0] = Temp+Temp-BOver3A; 884 Roots[1] = -Temp-BOver3A; 885 return( 2 ); 886 } 887 888 Desc = sqrt( Desc ); 889 Roots[0] = CubeRoot(Desc-QOver2)-CubeRoot(Desc+QOver2) - BOver3A; 890 891 return( 1 ); 892 } 893 894 895 #define MAX_SWEEPS 50 896 ob_make_rmat(double a[3][3],double rmat[9])897 void ob_make_rmat(double a[3][3],double rmat[9]) 898 { 899 /* 900 onorm, dnorm - hold the sum of diagonals and off diagonals to check Jacobi completion 901 d[3] - holds the diagonals of the input vector, which transofrm to become the Eigenvalues 902 r1, r2 - hold 1st two Eigenvectors 903 v1,v2,v3 - hold orthogonal unit vectors derived from Eigenvectors 904 905 The junction performs a Jacobi Eigenvalue/vector determination 906 (https://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm) on the supplied 907 Inertial Tensor in a, and returns a unit transform matrix rmat as a row matrix. 908 To work, a must be diagonally symmetric (i.e a[i][j] = a[j][i]) 909 v starts out holding the unit matrix (i.e. no transform in co-ordinate frame), 910 and undergoes the same rotations as applied to the Inertial Tensor in the Jacobi 911 process to arrive at the new co-ordinate frame. 912 Finally, the eigenvalues are sorted in order that the largest principal moment aligns to the 913 new x-axis 914 */ 915 double onorm, dnorm; 916 double b, dma, q, t, c, s,d[3]; 917 double atemp, vtemp, dtemp,v[3][3]; 918 double r1[3],r2[3],v1[3],v2[3],v3[3]; 919 int i, j, k, l; 920 921 memset((char*)d,'\0',sizeof(double)*3); 922 923 for (j = 0; j < 3; ++j) 924 { 925 for (i = 0; i < 3; ++i) 926 v[i][j] = 0.0; 927 928 v[j][j] = 1.0; 929 d[j] = a[j][j]; 930 } 931 932 for (l = 1; l <= MAX_SWEEPS; ++l) 933 { 934 dnorm = 0.0; 935 onorm = 0.0; 936 for (j = 0; j < 3; ++j) 937 { 938 dnorm = dnorm + (double)fabs(d[j]); 939 for (i = 0; i <= j - 1; ++i) 940 { 941 onorm = onorm + (double)fabs(a[i][j]); 942 } 943 } 944 945 if((onorm/dnorm) <= 1.0e-12) 946 /* Completion achieved (i.e. off-diagonals are all 0.0 within error)*/ 947 goto Exit_now; 948 for (j = 1; j < 3; ++j) 949 { 950 for (i = 0; i <= j - 1; ++i) 951 { 952 b = a[i][j]; 953 if(fabs(b) > 0.0) 954 { 955 dma = d[j] - d[i]; 956 if((fabs(dma) + fabs(b)) <= fabs(dma)) 957 t = b / dma; 958 else 959 { 960 q = 0.5 * dma / b; 961 t = 1.0/((double)fabs(q) + (double)sqrt(1.0+q*q)); 962 if(q < 0.0) 963 t = -t; 964 } 965 c = 1.0/(double)sqrt(t * t + 1.0); 966 s = t * c; 967 a[i][j] = 0.0; 968 /* Perform a Jacobi rotation on the supplied matrix*/ 969 for (k = 0; k <= i-1; ++k) 970 { 971 atemp = c * a[k][i] - s * a[k][j]; 972 a[k][j] = s * a[k][i] + c * a[k][j]; 973 a[k][i] = atemp; 974 } 975 for (k = i+1; k <= j-1; ++k) 976 { 977 atemp = c * a[i][k] - s * a[k][j]; 978 a[k][j] = s * a[i][k] + c * a[k][j]; 979 a[i][k] = atemp; 980 } 981 for (k = j+1; k < 3; ++k) 982 { 983 atemp = c * a[i][k] - s * a[j][k]; 984 a[j][k] = s * a[i][k] + c * a[j][k]; 985 a[i][k] = atemp; 986 } 987 /* Rotate the reference frame */ 988 for (k = 0; k < 3; ++k) 989 { 990 vtemp = c * v[k][i] - s * v[k][j]; 991 v[k][j] = s * v[k][i] + c * v[k][j]; 992 v[k][i] = vtemp; 993 } 994 dtemp = c*c*d[i] + s*s*d[j] - 2.0*c*s*b; 995 d[j] = s*s*d[i] + c*c*d[j] + 2.0*c*s*b; 996 d[i] = dtemp; 997 } /* end if */ 998 } /* end for i */ 999 } /* end for j */ 1000 } /* end for l */ 1001 1002 Exit_now: 1003 1004 /* max_sweeps = l;*/ 1005 1006 /* Now sort the eigenvalues and eigenvectors*/ 1007 for (j = 0; j < 3-1; ++j) 1008 { 1009 k = j; 1010 dtemp = d[k]; 1011 for (i = j+1; i < 3; ++i) 1012 if(d[i] < dtemp) 1013 { 1014 k = i; 1015 dtemp = d[k]; 1016 } 1017 1018 if(k > j) 1019 { 1020 d[k] = d[j]; 1021 d[j] = dtemp; 1022 for (i = 0; i < 3 ; ++i) 1023 { 1024 dtemp = v[i][k]; 1025 v[i][k] = v[i][j]; 1026 v[i][j] = dtemp; 1027 } 1028 } 1029 } 1030 1031 /* Transfer the 1st two eigenvectors into r1 and r2*/ 1032 r1[0] = v[0][0]; 1033 r1[1] = v[1][0]; 1034 r1[2] = v[2][0]; 1035 r2[0] = v[0][1]; 1036 r2[1] = v[1][1]; 1037 r2[2] = v[2][1]; 1038 1039 /* Generate the 3rd unit vector for the new co-ordinate frame by cross product of r1 and r2*/ 1040 v3[0] = r1[1]*r2[2] - r1[2]*r2[1]; 1041 v3[1] = -r1[0]*r2[2] + r1[2]*r2[0]; 1042 v3[2] = r1[0]*r2[1] - r1[1]*r2[0]; 1043 /* Ensure it is normalised |v3|=1 */ 1044 s = (double)sqrt(v3[0]*v3[0] + v3[1]*v3[1] + v3[2]*v3[2]); 1045 v3[0] /= s; 1046 v3[1] /= s; 1047 v3[2] /= s; 1048 1049 /* Generate the 2nd unit vector for the new co-ordinate frame by cross product of v3 and r1*/ 1050 v2[0] = v3[1]*r1[2] - v3[2]*r1[1]; 1051 v2[1] = -v3[0]*r1[2] + v3[2]*r1[0]; 1052 v2[2] = v3[0]*r1[1] - v3[1]*r1[0]; 1053 /* Ensure it is normalised |v2|=1 */ 1054 s = (double)sqrt(v2[0]*v2[0] + v2[1]*v2[1] + v2[2]*v2[2]); 1055 v2[0] /= s; 1056 v2[1] /= s; 1057 v2[2] /= s; 1058 1059 /* Generate the 1st unit vector for the new co-ordinate frame by cross product of v2 and v3*/ 1060 v1[0] = v2[1]*v3[2] - v2[2]*v3[1]; 1061 v1[1] = -v2[0]*v3[2] + v2[2]*v3[0]; 1062 v1[2] = v2[0]*v3[1] - v2[1]*v3[0]; 1063 /* Ensure it is normalised |v1|=1 */ 1064 s = (double)sqrt(v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2]); 1065 v1[0] /= s; 1066 v1[1] /= s; 1067 v1[2] /= s; 1068 1069 /* Transfer to the row matrix form for the result*/ 1070 rmat[0] = v1[0]; 1071 rmat[1] = v1[1]; 1072 rmat[2] = v1[2]; 1073 rmat[3] = v2[0]; 1074 rmat[4] = v2[1]; 1075 rmat[5] = v2[2]; 1076 rmat[6] = v3[0]; 1077 rmat[7] = v3[1]; 1078 rmat[8] = v3[2]; 1079 } 1080 get_roots_3_3(double mat[3][3],double roots[3])1081 static int get_roots_3_3(double mat[3][3], double roots[3]) 1082 { 1083 double rmat[9]; 1084 1085 ob_make_rmat(mat,rmat); 1086 1087 mat[0][0]=rmat[0]; 1088 mat[0][1]=rmat[3]; 1089 mat[0][2]=rmat[6]; 1090 mat[1][0]=rmat[1]; 1091 mat[1][1]=rmat[4]; 1092 mat[1][2]=rmat[7]; 1093 mat[2][0]=rmat[2]; 1094 mat[2][1]=rmat[5]; 1095 mat[2][2]=rmat[8]; 1096 1097 roots[0]=(double)Roots[0]; 1098 roots[1]=(double)Roots[1]; 1099 roots[2]=(double)Roots[2]; 1100 1101 return 1; 1102 } 1103 superimpose(double * r,double * f,int size)1104 double superimpose(double *r,double *f,int size) 1105 { 1106 int i,j; 1107 double x,y,z,d2; 1108 double mat[3][3],rmat[3][3],mat2[3][3],roots[3]; 1109 1110 /* make inertial cross tensor */ 1111 for(i=0;i<3;++i) 1112 for(j=0;j<3;++j) 1113 mat[i][j]=0.0; 1114 1115 for(i=0;i < size;++i) 1116 { 1117 mat[0][0]+=r[3*i] *f[3*i]; 1118 mat[1][0]+=r[3*i+1]*f[3*i]; 1119 mat[2][0]+=r[3*i+2]*f[3*i]; 1120 mat[0][1]+=r[3*i] *f[3*i+1]; 1121 mat[1][1]+=r[3*i+1]*f[3*i+1]; 1122 mat[2][1]+=r[3*i+2]*f[3*i+1]; 1123 mat[0][2]+=r[3*i] *f[3*i+2]; 1124 mat[1][2]+=r[3*i+1]*f[3*i+2]; 1125 mat[2][2]+=r[3*i+2]*f[3*i+2]; 1126 } 1127 1128 d2=mat[0][0]*(mat[1][1]*mat[2][2]-mat[1][2]*mat[2][1]) 1129 -mat[0][1]*(mat[1][0]*mat[2][2]-mat[1][2]*mat[2][0]) 1130 +mat[0][2]*(mat[1][0]*mat[2][1]-mat[1][1]*mat[2][0]); 1131 1132 1133 /* square matrix= ((mat transpose) * mat) */ 1134 for(i=0;i<3;++i) 1135 for(j=0;j<3;++j) 1136 { 1137 x=mat[0][i]*mat[0][j]+mat[1][i]*mat[1][j]+mat[2][i]*mat[2][j]; 1138 mat2[i][j]=mat[i][j]; 1139 rmat[i][j]=x; 1140 } 1141 get_roots_3_3(rmat,roots); 1142 1143 roots[0]=(roots[0]<0.0001) ? 0.0: (roots[0]); 1144 roots[1]=(roots[1]<0.0001) ? 0.0: (roots[1]); 1145 roots[2]=(roots[2]<0.0001) ? 0.0: (roots[2]); 1146 1147 /* make sqrt of rmat, store in mat*/ 1148 1149 roots[0]=roots[0]<0.0001? 0.0: 1.0/(double)sqrt(roots[0]); 1150 roots[1]=roots[1]<0.0001? 0.0: 1.0/(double)sqrt(roots[1]); 1151 roots[2]=roots[2]<0.0001? 0.0: 1.0/(double)sqrt(roots[2]); 1152 1153 if(d2<0.0) 1154 { 1155 if( (roots[0]>=roots[1]) && (roots[0]>=roots[2]) ) 1156 roots[0]*=-1.0; 1157 if( (roots[1]>roots[0]) && (roots[1]>=roots[2]) ) 1158 roots[1]*=-1.0; 1159 if( (roots[2]>roots[1]) && (roots[2]>roots[0]) ) 1160 roots[2]*=-1.0; 1161 } 1162 1163 for(i=0;i<3;++i) 1164 for(j=0;j<3;++j) 1165 mat[i][j]=roots[0]*rmat[i][0]*rmat[j][0]+ 1166 roots[1]*rmat[i][1]*rmat[j][1]+ 1167 roots[2]*rmat[i][2]*rmat[j][2]; 1168 1169 /* and multiply into original inertial cross matrix, mat2 */ 1170 for(i=0;i<3;++i) 1171 for(j=0;j<3;++j) 1172 rmat[i][j]=mat[0][j]*mat2[i][0]+ 1173 mat[1][j]*mat2[i][1]+ 1174 mat[2][j]*mat2[i][2]; 1175 1176 /* rotate all coordinates */ 1177 d2 = 0.0; 1178 for(i=0;i<size;++i) 1179 { 1180 x=f[3*i]*rmat[0][0]+f[3*i+1]*rmat[0][1]+f[3*i+2]*rmat[0][2]; 1181 y=f[3*i]*rmat[1][0]+f[3*i+1]*rmat[1][1]+f[3*i+2]*rmat[1][2]; 1182 z=f[3*i]*rmat[2][0]+f[3*i+1]*rmat[2][1]+f[3*i+2]*rmat[2][2]; 1183 f[3*i ]=x; 1184 f[3*i+1]=y; 1185 f[3*i+2]=z; 1186 1187 x = r[i*3] - f[i*3]; 1188 y = r[i*3+1] - f[i*3+1]; 1189 z = r[i*3+2] - f[i*3+2]; 1190 d2 += x*x+y*y+z*z; 1191 } 1192 1193 d2 /= (double) size; 1194 1195 return((double)sqrt(d2)); 1196 } 1197 get_rmat(double * rvec,double * r,double * f,int size)1198 void get_rmat(double *rvec,double *r,double *f,int size) 1199 { 1200 int i,j; 1201 double x,d2; 1202 double mat[3][3],rmat[3][3],mat2[3][3],roots[3]; 1203 1204 /* make inertial cross tensor */ 1205 for(i=0;i<3;++i) 1206 for(j=0;j<3;++j) 1207 mat[i][j]=0.0; 1208 1209 for(i=0;i < size;++i) 1210 { 1211 mat[0][0]+=r[3*i] *f[3*i]; 1212 mat[1][0]+=r[3*i+1]*f[3*i]; 1213 mat[2][0]+=r[3*i+2]*f[3*i]; 1214 mat[0][1]+=r[3*i] *f[3*i+1]; 1215 mat[1][1]+=r[3*i+1]*f[3*i+1]; 1216 mat[2][1]+=r[3*i+2]*f[3*i+1]; 1217 mat[0][2]+=r[3*i] *f[3*i+2]; 1218 mat[1][2]+=r[3*i+1]*f[3*i+2]; 1219 mat[2][2]+=r[3*i+2]*f[3*i+2]; 1220 } 1221 1222 d2=mat[0][0]*(mat[1][1]*mat[2][2]-mat[1][2]*mat[2][1]) 1223 -mat[0][1]*(mat[1][0]*mat[2][2]-mat[1][2]*mat[2][0]) 1224 +mat[0][2]*(mat[1][0]*mat[2][1]-mat[1][1]*mat[2][0]); 1225 1226 /* square matrix= ((mat transpose) * mat) */ 1227 for(i=0;i<3;++i) 1228 for(j=0;j<3;++j) 1229 { 1230 x=mat[0][i]*mat[0][j]+mat[1][i]*mat[1][j]+mat[2][i]*mat[2][j]; 1231 mat2[i][j]=mat[i][j]; 1232 rmat[i][j]=x; 1233 } 1234 get_roots_3_3(rmat,roots); 1235 1236 roots[0]=(roots[0]<0.0001) ? 0.0: (roots[0]); 1237 roots[1]=(roots[1]<0.0001) ? 0.0: (roots[1]); 1238 roots[2]=(roots[2]<0.0001) ? 0.0: (roots[2]); 1239 1240 /* make sqrt of rmat, store in mat*/ 1241 1242 roots[0]=(roots[0]<0.0001) ? 0.0: 1.0/(double)sqrt(roots[0]); 1243 roots[1]=(roots[1]<0.0001) ? 0.0: 1.0/(double)sqrt(roots[1]); 1244 roots[2]=(roots[2]<0.0001) ? 0.0: 1.0/(double)sqrt(roots[2]); 1245 1246 if(d2<0.0) 1247 { 1248 if( (roots[0]>=roots[1]) && (roots[0]>=roots[2]) ) 1249 roots[0]*=-1.0; 1250 if( (roots[1]>roots[0]) && (roots[1]>=roots[2]) ) 1251 roots[1]*=-1.0; 1252 if( (roots[2]>roots[1]) && (roots[2]>roots[0]) ) 1253 roots[2]*=-1.0; 1254 } 1255 1256 for(i=0;i<3;++i) 1257 for(j=0;j<3;++j) 1258 mat[i][j]=roots[0]*rmat[i][0]*rmat[j][0]+ 1259 roots[1]*rmat[i][1]*rmat[j][1]+ 1260 roots[2]*rmat[i][2]*rmat[j][2]; 1261 1262 /* and multiply into original inertial cross matrix, mat2 */ 1263 for(i=0;i<3;++i) 1264 for(j=0;j<3;++j) 1265 rmat[i][j]=mat[0][j]*mat2[i][0]+ 1266 mat[1][j]*mat2[i][1]+ 1267 mat[2][j]*mat2[i][2]; 1268 1269 rvec[0] = rmat[0][0]; 1270 rvec[1] = rmat[0][1]; 1271 rvec[2] = rmat[0][2]; 1272 rvec[3] = rmat[1][0]; 1273 rvec[4] = rmat[1][1]; 1274 rvec[5] = rmat[1][2]; 1275 rvec[6] = rmat[2][0]; 1276 rvec[7] = rmat[2][1]; 1277 rvec[8] = rmat[2][2]; 1278 } 1279 1280 } // end namespace OpenBabel 1281 1282 //! \file obutil.cpp 1283 //! \brief Various utility methods. 1284