| 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-2005 by Geoffrey R. Hutchison |
| 6 |
|
| 7 |
This file is part of the Open Babel project. |
| 8 |
For more information, see <http://openbabel.sourceforge.net/> |
| 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 "rotamer.hpp" |
| 21 |
#include "mol.hpp" |
| 22 |
#include "obutil.hpp" |
| 23 |
#include "rotor.hpp" |
| 24 |
|
| 25 |
#define OB_TITLE_SIZE 254 |
| 26 |
#define OB_BINARY_SETWORD 32 |
| 27 |
|
| 28 |
using namespace std; |
| 29 |
|
| 30 |
namespace OpenBabel |
| 31 |
{ |
| 32 |
|
| 33 |
//test byte ordering |
| 34 |
static int SINT = 0x00000001; |
| 35 |
static unsigned char *STPTR = (unsigned char*)&SINT; |
| 36 |
const bool SwabInt = (STPTR[0]!=0); |
| 37 |
|
| 38 |
#if !HAVE_RINT |
| 39 |
inline double rint(double x) |
| 40 |
{ |
| 41 |
return ( (x < 0.0) ? ceil(x-0.5) : floor(x+0.5)); |
| 42 |
} |
| 43 |
#endif |
| 44 |
|
| 45 |
void SetRotorToAngle(double *c,OBAtom **ref,double ang,vector<int> atoms); |
| 46 |
|
| 47 |
int Swab(int i) |
| 48 |
{ |
| 49 |
unsigned char tmp[4],c; |
| 50 |
memcpy(tmp,(char*)&i,sizeof(int)); |
| 51 |
c = tmp[0]; |
| 52 |
tmp[0] = tmp[3]; |
| 53 |
tmp[3] = c; |
| 54 |
c = tmp[1]; |
| 55 |
tmp[1] = tmp[2]; |
| 56 |
tmp[2] = c; |
| 57 |
memcpy((char*)&i,tmp,sizeof(int)); |
| 58 |
return(i); |
| 59 |
} |
| 60 |
|
| 61 |
OBRotamerList::~OBRotamerList() |
| 62 |
{ |
| 63 |
vector<unsigned char*>::iterator i; |
| 64 |
for (i = _vrotamer.begin();i != _vrotamer.end();i++) |
| 65 |
delete [] *i; |
| 66 |
|
| 67 |
vector<pair<OBAtom**,vector<int> > >::iterator j; |
| 68 |
for (j = _vrotor.begin();j != _vrotor.end();j++) |
| 69 |
delete [] j->first; |
| 70 |
|
| 71 |
//Delete the interal base coordinate list |
| 72 |
unsigned int k; |
| 73 |
for (k=0 ; k<_c.size() ; k++) |
| 74 |
delete [] _c[k]; |
| 75 |
} |
| 76 |
|
| 77 |
void OBRotamerList::GetReferenceArray(unsigned char *ref) |
| 78 |
{ |
| 79 |
int j; |
| 80 |
vector<pair<OBAtom**,vector<int> > >::iterator i; |
| 81 |
for (j=0,i = _vrotor.begin();i != _vrotor.end();i++) |
| 82 |
{ |
| 83 |
ref[j++] = (unsigned char)(i->first[0])->GetIdx(); |
| 84 |
ref[j++] = (unsigned char)(i->first[1])->GetIdx(); |
| 85 |
ref[j++] = (unsigned char)(i->first[2])->GetIdx(); |
| 86 |
ref[j++] = (unsigned char)(i->first[3])->GetIdx(); |
| 87 |
} |
| 88 |
} |
| 89 |
|
| 90 |
void OBRotamerList::Setup(OBMol &mol,OBRotorList &rl) |
| 91 |
{ |
| 92 |
//clear the old stuff out if necessary |
| 93 |
_vres.clear(); |
| 94 |
vector<unsigned char*>::iterator j; |
| 95 |
for (j = _vrotamer.begin();j != _vrotamer.end();j++) |
| 96 |
delete [] *j; |
| 97 |
_vrotamer.clear(); |
| 98 |
|
| 99 |
vector<pair<OBAtom**,vector<int> > >::iterator k; |
| 100 |
for (k = _vrotor.begin();k != _vrotor.end();k++) |
| 101 |
delete [] k->first; |
| 102 |
_vrotor.clear(); |
| 103 |
|
| 104 |
//create the new list |
| 105 |
OBRotor *rotor; |
| 106 |
vector<OBRotor*>::iterator i; |
| 107 |
vector<int> children; |
| 108 |
|
| 109 |
int ref[4]; |
| 110 |
OBAtom **atomlist; |
| 111 |
for (rotor = rl.BeginRotor(i);rotor;rotor = rl.NextRotor(i)) |
| 112 |
{ |
| 113 |
atomlist = new OBAtom* [4]; |
| 114 |
rotor->GetDihedralAtoms(ref); |
| 115 |
atomlist[0] = mol.GetAtom(ref[0]); |
| 116 |
atomlist[1] = mol.GetAtom(ref[1]); |
| 117 |
atomlist[2] = mol.GetAtom(ref[2]); |
| 118 |
atomlist[3] = mol.GetAtom(ref[3]); |
| 119 |
mol.FindChildren(children,ref[1],ref[2]); |
| 120 |
_vrotor.push_back(pair<OBAtom**,vector<int> > (atomlist,children)); |
| 121 |
_vres.push_back(rotor->GetResolution()); |
| 122 |
} |
| 123 |
|
| 124 |
vector<double>::iterator n; |
| 125 |
vector<vector<double> >::iterator m; |
| 126 |
for (m = _vres.begin();m != _vres.end();m++) |
| 127 |
for (n = m->begin();n != m->end();n++) |
| 128 |
*n *= RAD_TO_DEG; |
| 129 |
} |
| 130 |
|
| 131 |
void OBRotamerList::Setup(OBMol &mol,unsigned char *ref,int nrotors) |
| 132 |
{ |
| 133 |
//clear the old stuff out if necessary |
| 134 |
_vres.clear(); |
| 135 |
vector<unsigned char*>::iterator j; |
| 136 |
for (j = _vrotamer.begin();j != _vrotamer.end();j++) |
| 137 |
delete [] *j; |
| 138 |
_vrotamer.clear(); |
| 139 |
|
| 140 |
vector<pair<OBAtom**,vector<int> > >::iterator k; |
| 141 |
for (k = _vrotor.begin();k != _vrotor.end();k++) |
| 142 |
delete [] k->first; |
| 143 |
_vrotor.clear(); |
| 144 |
|
| 145 |
//create the new list |
| 146 |
int i; |
| 147 |
vector<int> children; |
| 148 |
|
| 149 |
int refatoms[4]; |
| 150 |
OBAtom **atomlist; |
| 151 |
for (i = 0; i < nrotors; i++) |
| 152 |
{ |
| 153 |
atomlist = new OBAtom* [4]; |
| 154 |
refatoms[0] = (int)ref[i*4 ]; |
| 155 |
refatoms[1] = (int)ref[i*4+1]; |
| 156 |
refatoms[2] = (int)ref[i*4+2]; |
| 157 |
refatoms[3] = (int)ref[i*4+3]; |
| 158 |
mol.FindChildren(children,refatoms[1],refatoms[2]); |
| 159 |
atomlist[0] = mol.GetAtom(refatoms[0]); |
| 160 |
atomlist[1] = mol.GetAtom(refatoms[1]); |
| 161 |
atomlist[2] = mol.GetAtom(refatoms[2]); |
| 162 |
atomlist[3] = mol.GetAtom(refatoms[3]); |
| 163 |
_vrotor.push_back(pair<OBAtom**,vector<int> > (atomlist,children)); |
| 164 |
} |
| 165 |
|
| 166 |
} |
| 167 |
|
| 168 |
void OBRotamerList::AddRotamer(double *c) |
| 169 |
{ |
| 170 |
int idx,size; |
| 171 |
double angle,res=255.0f/360.0f; |
| 172 |
vector3 v1,v2,v3,v4; |
| 173 |
|
| 174 |
unsigned char *rot = new unsigned char [_vrotor.size()+1]; |
| 175 |
rot[0] = (char) 0; |
| 176 |
|
| 177 |
vector<pair<OBAtom**,vector<int> > >::iterator i; |
| 178 |
for (size=1,i = _vrotor.begin();i != _vrotor.end();i++,size++) |
| 179 |
{ |
| 180 |
idx = (i->first[0])->GetCIdx(); |
| 181 |
v1.Set(c[idx],c[idx+1],c[idx+2]); |
| 182 |
idx = (i->first[1])->GetCIdx(); |
| 183 |
v2.Set(c[idx],c[idx+1],c[idx+2]); |
| 184 |
idx = (i->first[2])->GetCIdx(); |
| 185 |
v3.Set(c[idx],c[idx+1],c[idx+2]); |
| 186 |
idx = (i->first[3])->GetCIdx(); |
| 187 |
v4.Set(c[idx],c[idx+1],c[idx+2]); |
| 188 |
|
| 189 |
angle = CalcTorsionAngle(v1,v2,v3,v4); |
| 190 |
while (angle < 0.0f) |
| 191 |
angle += 360.0f; |
| 192 |
while (angle > 360.0f) |
| 193 |
angle -= 360.0f; |
| 194 |
rot[size] = (unsigned char)rint(angle*res); |
| 195 |
} |
| 196 |
|
| 197 |
_vrotamer.push_back(rot); |
| 198 |
} |
| 199 |
|
| 200 |
void OBRotamerList::AddRotamer(int *arr) |
| 201 |
{ |
| 202 |
unsigned int i; |
| 203 |
double angle,res=255.0f/360.0f; |
| 204 |
|
| 205 |
unsigned char *rot = new unsigned char [_vrotor.size()+1]; |
| 206 |
rot[0] = (unsigned char)arr[0]; |
| 207 |
|
| 208 |
for (i = 0;i < _vrotor.size();i++) |
| 209 |
{ |
| 210 |
angle = _vres[i][arr[i+1]]; |
| 211 |
while (angle < 0.0f) |
| 212 |
angle += 360.0f; |
| 213 |
while (angle > 360.0f) |
| 214 |
angle -= 360.0f; |
| 215 |
rot[i+1] = (unsigned char)rint(angle*res); |
| 216 |
} |
| 217 |
_vrotamer.push_back(rot); |
| 218 |
} |
| 219 |
|
| 220 |
void OBRotamerList::AddRotamer(unsigned char *arr) |
| 221 |
{ |
| 222 |
unsigned int i; |
| 223 |
double angle,res=255.0f/360.0f; |
| 224 |
|
| 225 |
unsigned char *rot = new unsigned char [_vrotor.size()+1]; |
| 226 |
rot[0] = (unsigned char)arr[0]; |
| 227 |
|
| 228 |
for (i = 0;i < _vrotor.size();i++) |
| 229 |
{ |
| 230 |
angle = _vres[i][(int)arr[i+1]]; |
| 231 |
while (angle < 0.0f) |
| 232 |
angle += 360.0f; |
| 233 |
while (angle > 360.0f) |
| 234 |
angle -= 360.0f; |
| 235 |
rot[i+1] = (unsigned char)rint(angle*res); |
| 236 |
} |
| 237 |
_vrotamer.push_back(rot); |
| 238 |
} |
| 239 |
|
| 240 |
void OBRotamerList::AddRotamers(unsigned char *arr,int nrotamers) |
| 241 |
{ |
| 242 |
unsigned int size; |
| 243 |
int i; |
| 244 |
|
| 245 |
size = (unsigned int)_vrotor.size()+1; |
| 246 |
for (i = 0;i < nrotamers;i++) |
| 247 |
{ |
| 248 |
unsigned char *rot = new unsigned char [size]; |
| 249 |
memcpy(rot,&arr[i*size],sizeof(char)*size); |
| 250 |
_vrotamer.push_back(rot); |
| 251 |
} |
| 252 |
} |
| 253 |
|
| 254 |
void OBRotamerList::ExpandConformerList(OBMol &mol,vector<double*> &clist) |
| 255 |
{ |
| 256 |
unsigned int j; |
| 257 |
double angle,invres=360.0f/255.0f; |
| 258 |
unsigned char *conf; |
| 259 |
vector<double*> tmpclist; |
| 260 |
vector<unsigned char*>::iterator i; |
| 261 |
|
| 262 |
for (i = _vrotamer.begin();i != _vrotamer.end();i++) |
| 263 |
{ |
| 264 |
conf = *i; |
| 265 |
double *c = new double [mol.NumAtoms()*3]; |
| 266 |
memcpy(c,clist[(int)conf[0]],sizeof(double)*mol.NumAtoms()*3); |
| 267 |
|
| 268 |
for (j = 0;j < _vrotor.size();j++) |
| 269 |
{ |
| 270 |
angle = invres*((double)conf[j+1]); |
| 271 |
if (angle > 180.0) |
| 272 |
angle -= 360.0; |
| 273 |
SetRotorToAngle(c,_vrotor[j].first,angle,_vrotor[j].second); |
| 274 |
} |
| 275 |
tmpclist.push_back(c); |
| 276 |
} |
| 277 |
|
| 278 |
//transfer the conf list |
| 279 |
vector<double*>::iterator k; |
| 280 |
for (k = clist.begin();k != clist.end();k++) |
| 281 |
delete [] *k; |
| 282 |
clist = tmpclist; |
| 283 |
} |
| 284 |
|
| 285 |
//! Create a conformer list using the internal base set of coordinates |
| 286 |
vector<double*> OBRotamerList::CreateConformerList(OBMol& mol) |
| 287 |
{ |
| 288 |
unsigned int j; |
| 289 |
double angle,invres=360.0f/255.0f; |
| 290 |
unsigned char *conf; |
| 291 |
vector<double*> tmpclist; |
| 292 |
vector<unsigned char*>::iterator i; |
| 293 |
|
| 294 |
for (i = _vrotamer.begin();i != _vrotamer.end();i++) |
| 295 |
{ |
| 296 |
conf = *i; |
| 297 |
double *c = new double [mol.NumAtoms()*3]; |
| 298 |
memcpy(c,_c[(int)conf[0]],sizeof(double)*mol.NumAtoms()*3); |
| 299 |
|
| 300 |
for (j = 0;j < _vrotor.size();j++) |
| 301 |
{ |
| 302 |
angle = invres*((double)conf[j+1]); |
| 303 |
if (angle > 180.0) |
| 304 |
angle -= 360.0; |
| 305 |
SetRotorToAngle(c,_vrotor[j].first,angle,_vrotor[j].second); |
| 306 |
} |
| 307 |
tmpclist.push_back(c); |
| 308 |
} |
| 309 |
|
| 310 |
return tmpclist; |
| 311 |
} |
| 312 |
|
| 313 |
//Copies the coordinates in bc, NOT the pointers, into the object |
| 314 |
void OBRotamerList::SetBaseCoordinateSets(vector<double*> bc, unsigned int N) |
| 315 |
{ |
| 316 |
unsigned int i,j; |
| 317 |
|
| 318 |
//Clear out old data |
| 319 |
for (i=0 ; i<_c.size() ; i++) |
| 320 |
delete [] _c[i]; |
| 321 |
_c.clear(); |
| 322 |
|
| 323 |
//Copy new data |
| 324 |
double *c = NULL; |
| 325 |
double *cc= NULL; |
| 326 |
for (i=0 ; i<bc.size() ; i++) |
| 327 |
{ |
| 328 |
c = new double [3*N]; |
| 329 |
cc = bc[i]; |
| 330 |
for (j=0 ; j<3*N ; j++) |
| 331 |
c[j] = cc[j]; |
| 332 |
_c.push_back(c); |
| 333 |
} |
| 334 |
_NBaseCoords = N; |
| 335 |
} |
| 336 |
|
| 337 |
//! Rotate the coordinates of 'atoms' |
| 338 |
//! such that tor == ang - atoms in 'tor' should be ordered such |
| 339 |
//! that the 3rd atom is the pivot around which atoms rotate |
| 340 |
//! ang is in degrees |
| 341 |
void SetRotorToAngle(double *c, OBAtom **ref,double ang,vector<int> atoms) |
| 342 |
{ |
| 343 |
double v1x,v1y,v1z,v2x,v2y,v2z,v3x,v3y,v3z; |
| 344 |
double c1x,c1y,c1z,c2x,c2y,c2z,c3x,c3y,c3z; |
| 345 |
double c1mag,c2mag,radang,costheta,m[9]; |
| 346 |
double x,y,z,mag,rotang,sn,cs,t,tx,ty,tz; |
| 347 |
|
| 348 |
int tor[4]; |
| 349 |
tor[0] = ref[0]->GetCIdx(); |
| 350 |
tor[1] = ref[1]->GetCIdx(); |
| 351 |
tor[2] = ref[2]->GetCIdx(); |
| 352 |
tor[3] = ref[3]->GetCIdx(); |
| 353 |
|
| 354 |
// |
| 355 |
//calculate the torsion angle |
| 356 |
// |
| 357 |
v1x = c[tor[0]] - c[tor[1]]; v2x = c[tor[1]] - c[tor[2]]; |
| 358 |
v1y = c[tor[0]+1] - c[tor[1]+1]; v2y = c[tor[1]+1] - c[tor[2]+1]; |
| 359 |
v1z = c[tor[0]+2] - c[tor[1]+2]; v2z = c[tor[1]+2] - c[tor[2]+2]; |
| 360 |
v3x = c[tor[2]] - c[tor[3]]; |
| 361 |
v3y = c[tor[2]+1] - c[tor[3]+1]; |
| 362 |
v3z = c[tor[2]+2] - c[tor[3]+2]; |
| 363 |
|
| 364 |
c1x = v1y*v2z - v1z*v2y; c2x = v2y*v3z - v2z*v3y; |
| 365 |
c1y = -v1x*v2z + v1z*v2x; c2y = -v2x*v3z + v2z*v3x; |
| 366 |
c1z = v1x*v2y - v1y*v2x; c2z = v2x*v3y - v2y*v3x; |
| 367 |
c3x = c1y*c2z - c1z*c2y; |
| 368 |
c3y = -c1x*c2z + c1z*c2x; |
| 369 |
c3z = c1x*c2y - c1y*c2x; |
| 370 |
|
| 371 |
c1mag = c1x*c1x + c1y*c1y + c1z*c1z; |
| 372 |
c2mag = c2x*c2x + c2y*c2y + c2z*c2z; |
| 373 |
if (c1mag*c2mag < 0.01) costheta = 1.0; //avoid div by zero error |
| 374 |
else costheta = (c1x*c2x + c1y*c2y + c1z*c2z)/(sqrt(c1mag*c2mag)); |
| 375 |
|
| 376 |
if (costheta < -0.999999) costheta = -0.999999f; |
| 377 |
if (costheta > 0.999999) costheta = 0.999999f; |
| 378 |
|
| 379 |
if ((v2x*c3x + v2y*c3y + v2z*c3z) > 0.0) radang = -acos(costheta); |
| 380 |
else radang = acos(costheta); |
| 381 |
|
| 382 |
// |
| 383 |
// now we have the torsion angle (radang) - set up the rot matrix |
| 384 |
// |
| 385 |
|
| 386 |
//find the difference between current and requested |
| 387 |
rotang = (DEG_TO_RAD*ang) - radang; |
| 388 |
|
| 389 |
sn = sin(rotang); cs = cos(rotang);t = 1 - cs; |
| 390 |
//normalize the rotation vector |
| 391 |
mag = sqrt(v2x*v2x + v2y*v2y + v2z*v2z); |
| 392 |
x = v2x/mag; y = v2y/mag; z = v2z/mag; |
| 393 |
|
| 394 |
//set up the rotation matrix |
| 395 |
m[0]= t*x*x + cs; m[1] = t*x*y + sn*z; m[2] = t*x*z - sn*y; |
| 396 |
m[3] = t*x*y - sn*z; m[4] = t*y*y + cs; m[5] = t*y*z + sn*x; |
| 397 |
m[6] = t*x*z + sn*y; m[7] = t*y*z - sn*x; m[8] = t*z*z + cs; |
| 398 |
|
| 399 |
// |
| 400 |
//now the matrix is set - time to rotate the atoms |
| 401 |
// |
| 402 |
tx = c[tor[1]];ty = c[tor[1]+1];tz = c[tor[1]+2]; |
| 403 |
vector<int>::iterator i;int j; |
| 404 |
for (i = atoms.begin();i != atoms.end();i++) |
| 405 |
{ |
| 406 |
j = ((*i)-1)*3; |
| 407 |
c[j] -= tx;c[j+1] -= ty;c[j+2]-= tz; |
| 408 |
x = c[j]*m[0] + c[j+1]*m[1] + c[j+2]*m[2]; |
| 409 |
y = c[j]*m[3] + c[j+1]*m[4] + c[j+2]*m[5]; |
| 410 |
z = c[j]*m[6] + c[j+1]*m[7] + c[j+2]*m[8]; |
| 411 |
c[j] = x; c[j+1] = y; c[j+2] = z; |
| 412 |
c[j] += tx;c[j+1] += ty;c[j+2] += tz; |
| 413 |
} |
| 414 |
} |
| 415 |
|
| 416 |
int PackCoordinate(double c[3],double max[3]) |
| 417 |
{ |
| 418 |
int tmp; |
| 419 |
float cf; |
| 420 |
cf = c[0]; |
| 421 |
tmp = ((int)(cf*max[0])) << 20; |
| 422 |
cf = c[1]; |
| 423 |
tmp |= ((int)(cf*max[1])) << 10; |
| 424 |
cf = c[2]; |
| 425 |
tmp |= ((int)(cf*max[2])); |
| 426 |
return(tmp); |
| 427 |
} |
| 428 |
|
| 429 |
void UnpackCoordinate(double c[3],double max[3],int tmp) |
| 430 |
{ |
| 431 |
float cf; |
| 432 |
cf = (float)(tmp>>20); |
| 433 |
c[0] = cf; |
| 434 |
c[0] *= max[0]; |
| 435 |
cf = (float)((tmp&0xffc00)>>10); |
| 436 |
c[1] = cf; |
| 437 |
c[1] *= max[1]; |
| 438 |
cf = (float)(tmp&0x3ff); |
| 439 |
c[2] = cf; |
| 440 |
c[2] *= max[2]; |
| 441 |
} |
| 442 |
|
| 443 |
} //namespace OpenBabel |
| 444 |
|
| 445 |
//! \file rotamer.cpp |
| 446 |
//! \brief Handle rotamer list data. |
