1 /********************************************************************** 2 rotamer.cpp - Handle rotamer list data. 3 4 Copyright (C) 1998, 1999, 2000-2002 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 #include <openbabel/babelconfig.h> 20 21 #include <openbabel/mol.h> 22 #include <openbabel/atom.h> 23 #include <openbabel/ring.h> 24 #include <openbabel/obiter.h> 25 #include <openbabel/rotamer.h> 26 27 #define OB_TITLE_SIZE 254 28 #define OB_BINARY_SETWORD 32 29 30 using namespace std; 31 32 namespace OpenBabel 33 { 34 35 /** \class OBRotamerList rotamer.h <openbabel/rotamer.h> 36 37 A high-level class for rotamer / conformer generation, intended mainly 38 for use with the related class OBRotorList and the OBRotorRules database 39 40 Rotamers represent conformational isomers formed simply by rotation of 41 dihedral angles. On the other hand, conformers may include geometric 42 relaxation (i.e., slight modification of bond lengths, bond angles, etc.) 43 44 The following shows an example of generating 2 conformers using different 45 rotor states. Similar code could be used for systematic or Monte Carlo 46 conformer sampling when combined with energy evaluation (molecular 47 mechanics or otherwise). 48 49 \code 50 OBRotorList rl; // used to sample all rotatable bonds via the OBRotorRules data 51 // If you want to "fix" any particular atoms (i.e., freeze them in space) 52 // then set up an OBBitVec of the fixed atoms and call 53 // rl.SetFixAtoms(bitvec); 54 rl.Setup(mol); 55 56 // How many rotatable bonds are there? 57 cerr << " Number of rotors: " << rl.Size() << endl; 58 59 // indexed from 1, rotorKey[0] = 0 60 std::vector<int> rotorKey(rl.Size() + 1, 0); 61 62 // each entry represents the configuration of a rotor 63 // e.g. indexes into OBRotor::GetResolution() -- the different angles 64 // to sample for a rotamer search 65 for (unsigned int i = 0; i < rl.Size() + 1; ++i) 66 rotorKey[i] = 0; // could be anything from 0 .. OBRotor->GetResolution().size() 67 // -1 is for no rotation 68 69 // The OBRotamerList can generate conformations (i.e., coordinate sets) 70 OBRotamerList rotamers; 71 rotamers.SetBaseCoordinateSets(mol); 72 rotamers.Setup(mol, rl); 73 74 rotamers.AddRotamer(rotorKey); 75 rotorKey[1] = 2; // switch one rotor 76 rotamers.AddRotamer(rotorKey); 77 78 rotamers.ExpandConformerList(mol, mol.GetConformers()); 79 80 // change the molecule conformation 81 mol.SetConformer(0); // rotorKey 0, 0, ... 82 conv.Write(&mol); 83 84 mol.SetConformer(1); // rotorKey 0, 2, ... 85 86 \endcode 87 88 **/ 89 90 //test byte ordering 91 static int SINT = 0x00000001; 92 static unsigned char *STPTR = (unsigned char*)&SINT; 93 const bool SwabInt = (STPTR[0]!=0); 94 95 #if !HAVE_RINT rint(double x)96 inline double rint(double x) 97 { 98 return ( (x < 0.0) ? ceil(x-0.5) : floor(x+0.5)); 99 } 100 #endif 101 102 void SetRotorToAngle(double *c,OBAtom **ref,double ang,vector<int> atoms); 103 Swab(int i)104 int Swab(int i) 105 { 106 unsigned char tmp[4],c; 107 memcpy(tmp,(char*)&i,sizeof(int)); 108 c = tmp[0]; 109 tmp[0] = tmp[3]; 110 tmp[3] = c; 111 c = tmp[1]; 112 tmp[1] = tmp[2]; 113 tmp[2] = c; 114 memcpy((char*)&i,tmp,sizeof(int)); 115 return(i); 116 } 117 Clone(OBBase * newparent) const118 OBGenericData* OBRotamerList::Clone(OBBase* newparent) const 119 { 120 //Since the class contains OBAtom pointers, the new copy use 121 //these from the new molecule, newparent 122 OBMol* newmol = static_cast<OBMol*>(newparent); 123 124 OBRotamerList *new_rml = new OBRotamerList; 125 new_rml->_attr = _attr; 126 new_rml->_type = _type; 127 128 //Set base coordinates 129 unsigned int k,l; 130 vector<double*> bc; 131 double *c = nullptr; 132 double *cc = nullptr; 133 for (k=0 ; k<NumBaseCoordinateSets() ; ++k) 134 { 135 c = new double [3*NumAtoms()]; 136 cc = GetBaseCoordinateSet(k); 137 for (l=0 ; l<3*NumAtoms() ; ++l) 138 c[l] = cc[l]; 139 bc.push_back(c); 140 } 141 if (NumBaseCoordinateSets()) 142 new_rml->SetBaseCoordinateSets(bc,NumAtoms()); 143 144 //Set reference array 145 unsigned char *ref = new unsigned char [NumRotors()*4]; 146 if (ref) 147 { 148 GetReferenceArray(ref); 149 new_rml->Setup(*newmol,ref,NumRotors()); 150 delete [] ref; 151 } 152 153 //Set Rotamers 154 unsigned char *rotamers = new unsigned char [(NumRotors()+1)*NumRotamers()]; 155 if (rotamers) 156 { 157 vector<unsigned char*>::const_iterator kk; 158 unsigned int idx=0; 159 for (kk = _vrotamer.begin();kk != _vrotamer.end();++kk) 160 { 161 memcpy(&rotamers[idx],(const unsigned char*)*kk,sizeof(unsigned char)*(NumRotors()+1)); 162 idx += sizeof(unsigned char)*(NumRotors()+1); 163 } 164 new_rml->AddRotamers(rotamers,NumRotamers()); 165 delete [] rotamers; 166 } 167 return new_rml; 168 } 169 ~OBRotamerList()170 OBRotamerList::~OBRotamerList() 171 { 172 vector<unsigned char*>::iterator i; 173 for (i = _vrotamer.begin();i != _vrotamer.end();++i) 174 delete [] *i; 175 176 vector<pair<OBAtom**,vector<int> > >::iterator j; 177 for (j = _vrotor.begin();j != _vrotor.end();++j) 178 delete [] j->first; 179 180 //Delete the interal base coordinate list 181 unsigned int k; 182 for (k=0 ; k<_c.size() ; ++k) 183 delete [] _c[k]; 184 } 185 GetReferenceArray(unsigned char * ref) const186 void OBRotamerList::GetReferenceArray(unsigned char *ref)const 187 { 188 int j; 189 vector<pair<OBAtom**,vector<int> > >::const_iterator i; 190 for (j=0,i = _vrotor.begin();i != _vrotor.end();++i) 191 { 192 ref[j++] = (unsigned char)(i->first[0])->GetIdx(); 193 ref[j++] = (unsigned char)(i->first[1])->GetIdx(); 194 ref[j++] = (unsigned char)(i->first[2])->GetIdx(); 195 ref[j++] = (unsigned char)(i->first[3])->GetIdx(); 196 } 197 } 198 Setup(OBMol & mol,OBRotorList & rl)199 void OBRotamerList::Setup(OBMol &mol,OBRotorList &rl) 200 { 201 //clear the old stuff out if necessary 202 _vres.clear(); 203 vector<unsigned char*>::iterator j; 204 for (j = _vrotamer.begin();j != _vrotamer.end();++j) 205 delete [] *j; 206 _vrotamer.clear(); 207 208 vector<pair<OBAtom**,vector<int> > >::iterator k; 209 for (k = _vrotor.begin();k != _vrotor.end();++k) 210 delete [] k->first; 211 _vrotor.clear(); 212 _vrings.clear(); 213 _vringTors.clear(); 214 215 //create the new list 216 OBRotor *rotor; 217 vector<OBRotor*>::iterator i; 218 vector<int> children; 219 220 int ref[4]; 221 OBAtom **atomlist; 222 for (rotor = rl.BeginRotor(i);rotor;rotor = rl.NextRotor(i)) 223 { 224 atomlist = new OBAtom* [4]; 225 rotor->GetDihedralAtoms(ref); 226 atomlist[0] = mol.GetAtom(ref[0]); 227 atomlist[1] = mol.GetAtom(ref[1]); 228 atomlist[2] = mol.GetAtom(ref[2]); 229 atomlist[3] = mol.GetAtom(ref[3]); 230 mol.FindChildren(children,ref[1],ref[2]); 231 _vrotor.push_back(pair<OBAtom**,vector<int> > (atomlist,children)); 232 _vres.push_back(rotor->GetResolution()); 233 } 234 235 // if the rotor list has ring bonds, build up an index 236 if (rl.HasRingRotors()){ 237 // go through rings 238 // for each step of the path, see if there's a matching rotor 239 vector<int> path; 240 int pSize; 241 vector<double> ringTorsions; 242 vector<int> ringRotors; 243 FOR_RINGS_OF_MOL(r, mol) 244 { 245 ringTorsions.clear(); 246 ringRotors.clear(); 247 248 pSize = r->Size(); 249 if (pSize < 4) 250 continue; // not rotatable 251 252 path = r->_path; 253 for (int j = 0; j < pSize; ++j) { 254 double torsion = mol.GetTorsion(path[(j + pSize - 1) % pSize], 255 path[(j + pSize) % pSize], 256 path[(j + pSize + 1) % pSize], 257 path[(j + pSize + 2) % pSize]); 258 ringTorsions.push_back(torsion); 259 260 // now check to see if any of these things are rotors 261 int rotorIndex = -1; // not a rotor 262 OBBond *bond = mol.GetBond(path[(j + pSize) % pSize], path[(j + pSize + 1) % pSize]); 263 for (rotor = rl.BeginRotor(i);rotor;rotor = rl.NextRotor(i)) 264 { 265 if (bond != rotor->GetBond()) 266 continue; // no match at all 267 268 // Central bond matches, make sure 1..4 atoms are in the path 269 rotor->GetDihedralAtoms(ref); 270 if ( (ref[0] == path[(j + pSize - 1) % pSize] && 271 ref[3] == path[(j + pSize + 2) % pSize]) 272 || 273 (ref[3] == path[(j + pSize - 1) % pSize] && 274 ref[0] == path[(j + pSize + 2) % pSize]) ) { 275 rotorIndex = rotor->GetIdx(); 276 } 277 } // end checking all the rotors 278 ringRotors.push_back(rotorIndex); // could be -1 if it's not rotatable 279 } 280 281 _vringTors.push_back(ringTorsions); 282 _vrings.push_back(ringRotors); 283 } // finished with the rings 284 } // if (ring rotors) 285 286 vector<double>::iterator n; 287 vector<vector<double> >::iterator m; 288 for (m = _vres.begin();m != _vres.end();++m) 289 for (n = m->begin();n != m->end();++n) 290 *n *= RAD_TO_DEG; 291 } 292 Setup(OBMol & mol,unsigned char * ref,int nrotors)293 void OBRotamerList::Setup(OBMol &mol,unsigned char *ref,int nrotors) 294 { 295 //clear the old stuff out if necessary 296 _vres.clear(); 297 vector<unsigned char*>::iterator j; 298 for (j = _vrotamer.begin();j != _vrotamer.end();++j) 299 delete [] *j; 300 _vrotamer.clear(); 301 302 vector<pair<OBAtom**,vector<int> > >::iterator k; 303 for (k = _vrotor.begin();k != _vrotor.end();++k) 304 delete [] k->first; 305 _vrotor.clear(); 306 _vrings.clear(); 307 _vringTors.clear(); 308 309 //create the new list 310 int i; 311 vector<int> children; 312 313 int refatoms[4]; 314 OBAtom **atomlist; 315 for (i = 0; i < nrotors; ++i) 316 { 317 atomlist = new OBAtom* [4]; 318 refatoms[0] = (int)ref[i*4 ]; 319 refatoms[1] = (int)ref[i*4+1]; 320 refatoms[2] = (int)ref[i*4+2]; 321 refatoms[3] = (int)ref[i*4+3]; 322 mol.FindChildren(children,refatoms[1],refatoms[2]); 323 atomlist[0] = mol.GetAtom(refatoms[0]); 324 atomlist[1] = mol.GetAtom(refatoms[1]); 325 atomlist[2] = mol.GetAtom(refatoms[2]); 326 atomlist[3] = mol.GetAtom(refatoms[3]); 327 _vrotor.push_back(pair<OBAtom**,vector<int> > (atomlist,children)); 328 } 329 330 } 331 AddRotamer(double * c)332 void OBRotamerList::AddRotamer(double *c) 333 { 334 // TODO: consider implementing periodic boundary conditions 335 int idx,size; 336 double angle,res=255.0/360.0; 337 vector3 v1,v2,v3,v4; 338 339 unsigned char *rot = new unsigned char [_vrotor.size()+1]; 340 rot[0] = (char) 0; 341 342 vector<pair<OBAtom**,vector<int> > >::iterator i; 343 for (size=1,i = _vrotor.begin();i != _vrotor.end();++i,++size) 344 { 345 idx = (i->first[0])->GetCoordinateIdx(); 346 v1.Set(c[idx],c[idx+1],c[idx+2]); 347 idx = (i->first[1])->GetCoordinateIdx(); 348 v2.Set(c[idx],c[idx+1],c[idx+2]); 349 idx = (i->first[2])->GetCoordinateIdx(); 350 v3.Set(c[idx],c[idx+1],c[idx+2]); 351 idx = (i->first[3])->GetCoordinateIdx(); 352 v4.Set(c[idx],c[idx+1],c[idx+2]); 353 354 angle = CalcTorsionAngle(v1,v2,v3,v4); 355 while (angle < 0.0) 356 angle += 360.0; 357 while (angle > 360.0) 358 angle -= 360.0; 359 rot[size] = (unsigned char)rint(angle*res); 360 } 361 362 _vrotamer.push_back(rot); 363 } 364 AddRotamer(int * arr)365 void OBRotamerList::AddRotamer(int *arr) 366 { 367 unsigned int i; 368 double angle,res=255.0/360.0; 369 370 unsigned char *rot = new unsigned char [_vrotor.size()+1]; 371 rot[0] = (unsigned char)arr[0]; 372 373 for (i = 0;i < _vrotor.size();++i) 374 { 375 angle = _vres[i][arr[i+1]]; 376 while (angle < 0.0) 377 angle += 360.0; 378 while (angle > 360.0) 379 angle -= 360.0; 380 rot[i+1] = (unsigned char)rint(angle*res); 381 } 382 _vrotamer.push_back(rot); 383 } 384 AddRotamer(std::vector<int> arr)385 void OBRotamerList::AddRotamer(std::vector<int> arr) 386 { 387 unsigned int i; 388 double angle,res=255.0/360.0; 389 390 if (arr.size() != (_vrotor.size() + 1)) 391 return; // wrong size key 392 393 // gotta check for weird ring torsion combinations 394 if (_vrings.size()) { 395 // go through each ring and update the possible torsions 396 for (unsigned int j = 0; j < _vrings.size(); ++j) { 397 vector<int> path = _vrings[j]; 398 double torsionSum = 0.0; 399 400 // go around the loop and add up the torsions 401 for (unsigned int i = 0; i < path.size(); ++i) { 402 if (path[i] == -1) { 403 // not a rotor, use the fixed value 404 torsionSum += _vringTors[j][i]; 405 continue; 406 } 407 408 // what angle are we trying to use with this key? 409 angle = _vres[ path[i] ][arr[ path[i]+1 ]]*res; 410 while (angle < 0.0) 411 angle += 360.0; 412 while (angle > 360.0) 413 angle -= 360.0; 414 415 // update the ring torsion for this setting 416 _vringTors[j][i] = angle; 417 torsionSum += angle; 418 } 419 420 // if the sum of the ring torsions is not ~0, bad move 421 if (fabs(torsionSum) > 45.0) { 422 // cerr << " Bad move! " << fabs(torsionSum) << endl; 423 return; // don't make the move 424 } 425 // cerr << " Good move!" << endl; 426 } 427 } 428 429 unsigned char *rot = new unsigned char [_vrotor.size()+1]; 430 rot[0] = (unsigned char)arr[0]; 431 432 for (i = 0;i < _vrotor.size();++i) 433 { 434 angle = _vres[i][arr[i+1]]; 435 while (angle < 0.0) 436 angle += 360.0; 437 while (angle > 360.0) 438 angle -= 360.0; 439 rot[i+1] = (unsigned char)rint(angle*res); 440 } 441 _vrotamer.push_back(rot); 442 } 443 AddRotamer(unsigned char * arr)444 void OBRotamerList::AddRotamer(unsigned char *arr) 445 { 446 unsigned int i; 447 double angle,res=255.0/360.0; 448 449 unsigned char *rot = new unsigned char [_vrotor.size()+1]; 450 rot[0] = (unsigned char)arr[0]; 451 452 for (i = 0;i < _vrotor.size();++i) 453 { 454 angle = _vres[i][(int)arr[i+1]]; 455 while (angle < 0.0) 456 angle += 360.0; 457 while (angle > 360.0) 458 angle -= 360.0; 459 rot[i+1] = (unsigned char)rint(angle*res); 460 } 461 _vrotamer.push_back(rot); 462 } 463 AddRotamers(unsigned char * arr,int nrotamers)464 void OBRotamerList::AddRotamers(unsigned char *arr,int nrotamers) 465 { 466 unsigned int size; 467 int i; 468 469 size = (unsigned int)_vrotor.size()+1; 470 for (i = 0;i < nrotamers;++i) 471 { 472 unsigned char *rot = new unsigned char [size]; 473 memcpy(rot,&arr[i*size],sizeof(char)*size); 474 _vrotamer.push_back(rot); 475 } 476 } 477 ExpandConformerList(OBMol & mol,vector<double * > & clist)478 void OBRotamerList::ExpandConformerList(OBMol &mol,vector<double*> &clist) 479 { 480 vector<double*> tmpclist = CreateConformerList(mol); 481 482 //transfer the conf list 483 vector<double*>::iterator k; 484 for (k = clist.begin();k != clist.end();++k) 485 delete [] *k; 486 clist = tmpclist; 487 } 488 489 //! Create a conformer list using the internal base set of coordinates CreateConformerList(OBMol & mol)490 vector<double*> OBRotamerList::CreateConformerList(OBMol& mol) 491 { 492 unsigned int j; 493 double angle,invres=360.0/255.0; 494 unsigned char *conf; 495 vector<double*> tmpclist; 496 vector<unsigned char*>::iterator i; 497 498 for (i = _vrotamer.begin();i != _vrotamer.end();++i) 499 { 500 conf = *i; 501 double *c = new double [mol.NumAtoms()*3]; 502 memcpy(c,_c[(int)conf[0]],sizeof(double)*mol.NumAtoms()*3); 503 504 for (j = 0;j < _vrotor.size();++j) 505 { 506 angle = invres*((double)conf[j+1]); 507 if (angle > 180.0) 508 angle -= 360.0; 509 SetRotorToAngle(c,_vrotor[j].first,angle,_vrotor[j].second); 510 } 511 tmpclist.push_back(c); 512 } 513 514 return tmpclist; 515 } 516 517 //! Change the current coordinate set 518 //! \since version 2.2 SetCurrentCoordinates(OBMol & mol,std::vector<int> arr)519 bool OBRotamerList::SetCurrentCoordinates(OBMol &mol, std::vector<int> arr) 520 { 521 double angle; 522 523 if (arr.size() != (_vrotor.size() + 1)) 524 return false; // wrong size key 525 526 // gotta check for weird ring torsion combinations 527 if (_vrings.size()) { 528 // go through each ring and update the possible torsions 529 for (unsigned int j = 0; j < _vrings.size(); ++j) { 530 vector<int> path = _vrings[j]; 531 double torsionSum = 0.0; 532 533 // go around the loop and add up the torsions 534 for (unsigned int i = 0; i < path.size(); ++i) { 535 if (path[i] == -1) { 536 // not a rotor, use the fixed value 537 torsionSum += _vringTors[j][i]; 538 continue; 539 } 540 541 // what angle are we trying to use with this key? 542 angle = _vres[ path[i] ][arr[ path[i]+1 ]]; 543 while (angle < 0.0) 544 angle += 360.0; 545 while (angle > 360.0) 546 angle -= 360.0; 547 548 // update the ring torsion for this setting 549 _vringTors[j][i] = angle; 550 torsionSum += angle; 551 } 552 553 // if the sum of the ring torsions is not ~0, bad move 554 if (fabs(torsionSum) > 45.0) { 555 // cerr << " Bad move!" << endl; 556 return false; // don't make the move 557 } 558 } 559 } 560 561 double *c = mol.GetCoordinates(); 562 for (unsigned int i = 0;i < _vrotor.size();++i) 563 { 564 if (arr[i+1] == -1) // skip this rotor 565 continue; 566 else { 567 angle = _vres[i][arr[i+1]]; 568 while (angle < 0.0) 569 angle += 360.0; 570 while (angle > 360.0) 571 angle -= 360.0; 572 SetRotorToAngle(c,_vrotor[i].first,angle,_vrotor[i].second); 573 } // set an angle 574 } // for rotors 575 576 return true; 577 } 578 SetBaseCoordinateSets(OBMol & mol)579 void OBRotamerList::SetBaseCoordinateSets(OBMol& mol) 580 { 581 SetBaseCoordinateSets(mol.GetConformers(), mol.NumAtoms()); 582 } 583 584 //Copies the coordinates in bc, NOT the pointers, into the object SetBaseCoordinateSets(vector<double * > bc,unsigned int N)585 void OBRotamerList::SetBaseCoordinateSets(vector<double*> bc, unsigned int N) 586 { 587 unsigned int i,j; 588 589 //Clear out old data 590 for (i=0 ; i<_c.size() ; ++i) 591 delete [] _c[i]; 592 _c.clear(); 593 594 //Copy new data 595 double *c = nullptr; 596 double *cc = nullptr; 597 for (i=0 ; i<bc.size() ; ++i) 598 { 599 c = new double [3*N]; 600 cc = bc[i]; 601 for (j=0 ; j<3*N ; ++j) 602 c[j] = cc[j]; 603 _c.push_back(c); 604 } 605 _NBaseCoords = N; 606 } 607 608 //! Rotate the coordinates of 'atoms' 609 //! such that tor == ang. 610 //! Atoms in 'tor' should be ordered such that the 3rd atom is 611 //! the pivot around which atoms rotate (ang is in degrees) 612 //! \todo This code is identical to OBMol::SetTorsion() and should be combined SetRotorToAngle(double * c,OBAtom ** ref,double ang,vector<int> atoms)613 void SetRotorToAngle(double *c, OBAtom **ref,double ang,vector<int> atoms) 614 { 615 double v1x,v1y,v1z,v2x,v2y,v2z,v3x,v3y,v3z; 616 double c1x,c1y,c1z,c2x,c2y,c2z,c3x,c3y,c3z; 617 double c1mag,c2mag,radang,costheta,m[9]; 618 double x,y,z,mag,rotang,sn,cs,t,tx,ty,tz; 619 620 int tor[4]; 621 tor[0] = ref[0]->GetCoordinateIdx(); 622 tor[1] = ref[1]->GetCoordinateIdx(); 623 tor[2] = ref[2]->GetCoordinateIdx(); 624 tor[3] = ref[3]->GetCoordinateIdx(); 625 626 // 627 //calculate the torsion angle 628 // 629 v1x = c[tor[0]] - c[tor[1]]; v2x = c[tor[1]] - c[tor[2]]; 630 v1y = c[tor[0]+1] - c[tor[1]+1]; v2y = c[tor[1]+1] - c[tor[2]+1]; 631 v1z = c[tor[0]+2] - c[tor[1]+2]; v2z = c[tor[1]+2] - c[tor[2]+2]; 632 v3x = c[tor[2]] - c[tor[3]]; 633 v3y = c[tor[2]+1] - c[tor[3]+1]; 634 v3z = c[tor[2]+2] - c[tor[3]+2]; 635 636 c1x = v1y*v2z - v1z*v2y; c2x = v2y*v3z - v2z*v3y; 637 c1y = -v1x*v2z + v1z*v2x; c2y = -v2x*v3z + v2z*v3x; 638 c1z = v1x*v2y - v1y*v2x; c2z = v2x*v3y - v2y*v3x; 639 c3x = c1y*c2z - c1z*c2y; 640 c3y = -c1x*c2z + c1z*c2x; 641 c3z = c1x*c2y - c1y*c2x; 642 643 c1mag = c1x*c1x + c1y*c1y + c1z*c1z; 644 c2mag = c2x*c2x + c2y*c2y + c2z*c2z; 645 if (c1mag*c2mag < 0.01) costheta = 1.0; //avoid div by zero error 646 else costheta = (c1x*c2x + c1y*c2y + c1z*c2z)/(sqrt(c1mag*c2mag)); 647 648 if (costheta < -0.999999) costheta = -0.999999; 649 if (costheta > 0.999999) costheta = 0.999999; 650 651 if ((v2x*c3x + v2y*c3y + v2z*c3z) > 0.0) radang = -acos(costheta); 652 else radang = acos(costheta); 653 654 // 655 // now we have the torsion angle (radang) - set up the rot matrix 656 // 657 658 //find the difference between current and requested 659 rotang = (DEG_TO_RAD*ang) - radang; 660 661 sn = sin(rotang); cs = cos(rotang);t = 1 - cs; 662 //normalize the rotation vector 663 mag = sqrt(v2x*v2x + v2y*v2y + v2z*v2z); 664 if (mag < 0.1) mag = 0.1; // avoid divide by zero error 665 x = v2x/mag; y = v2y/mag; z = v2z/mag; 666 667 //set up the rotation matrix 668 m[0]= t*x*x + cs; m[1] = t*x*y + sn*z; m[2] = t*x*z - sn*y; 669 m[3] = t*x*y - sn*z; m[4] = t*y*y + cs; m[5] = t*y*z + sn*x; 670 m[6] = t*x*z + sn*y; m[7] = t*y*z - sn*x; m[8] = t*z*z + cs; 671 672 // 673 //now the matrix is set - time to rotate the atoms 674 // 675 tx = c[tor[1]];ty = c[tor[1]+1];tz = c[tor[1]+2]; 676 vector<int>::iterator i;int j; 677 for (i = atoms.begin();i != atoms.end();++i) 678 { 679 j = ((*i)-1)*3; 680 c[j] -= tx;c[j+1] -= ty;c[j+2]-= tz; 681 x = c[j]*m[0] + c[j+1]*m[1] + c[j+2]*m[2]; 682 y = c[j]*m[3] + c[j+1]*m[4] + c[j+2]*m[5]; 683 z = c[j]*m[6] + c[j+1]*m[7] + c[j+2]*m[8]; 684 c[j] = x; c[j+1] = y; c[j+2] = z; 685 c[j] += tx;c[j+1] += ty;c[j+2] += tz; 686 } 687 } 688 PackCoordinate(double c[3],double max[3])689 int PackCoordinate(double c[3],double max[3]) 690 { 691 int tmp; 692 double cf; 693 cf = c[0]; 694 tmp = ((int)(cf*max[0])) << 20; 695 cf = c[1]; 696 tmp |= ((int)(cf*max[1])) << 10; 697 cf = c[2]; 698 tmp |= ((int)(cf*max[2])); 699 return(tmp); 700 } 701 UnpackCoordinate(double c[3],double max[3],int tmp)702 void UnpackCoordinate(double c[3],double max[3],int tmp) 703 { 704 double cf; 705 cf = (double)(tmp>>20); 706 c[0] = cf; 707 c[0] *= max[0]; 708 cf = (double)((tmp&0xffc00)>>10); 709 c[1] = cf; 710 c[1] *= max[1]; 711 cf = (double)(tmp&0x3ff); 712 c[2] = cf; 713 c[2] *= max[2]; 714 } 715 716 } //namespace OpenBabel 717 718 //! \file rotamer.cpp 719 //! \brief Handle rotamer list data. 720