Coverage Report

Created: 2025-08-26 06:55

/src/openbabel/src/rotamer.cpp
Line
Count
Source (jump to first uncovered line)
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
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
104
  int Swab(int i)
105
0
  {
106
0
    unsigned char tmp[4],c;
107
0
    memcpy(tmp,(char*)&i,sizeof(int));
108
0
    c = tmp[0];
109
0
    tmp[0] = tmp[3];
110
0
    tmp[3] = c;
111
0
    c = tmp[1];
112
0
    tmp[1] = tmp[2];
113
0
    tmp[2] = c;
114
0
    memcpy((char*)&i,tmp,sizeof(int));
115
0
    return(i);
116
0
  }
117
118
  OBGenericData* OBRotamerList::Clone(OBBase* newparent) const
119
0
  {
120
    //Since the class contains OBAtom pointers, the new copy use
121
    //these from the new molecule, newparent
122
0
    OBMol* newmol = static_cast<OBMol*>(newparent);
123
124
0
    OBRotamerList *new_rml = new OBRotamerList;
125
0
    new_rml->_attr = _attr;
126
0
    new_rml->_type = _type;
127
128
    //Set base coordinates
129
0
    unsigned int k,l;
130
0
    vector<double*> bc;
131
0
    double *c = nullptr;
132
0
    double *cc = nullptr;
133
0
    for (k=0 ; k<NumBaseCoordinateSets() ; ++k)
134
0
      {
135
0
        c = new double [3*NumAtoms()];
136
0
        cc = GetBaseCoordinateSet(k);
137
0
        for (l=0 ; l<3*NumAtoms() ; ++l)
138
0
          c[l] = cc[l];
139
0
        bc.push_back(c);
140
0
      }
141
0
    if (NumBaseCoordinateSets())
142
0
      new_rml->SetBaseCoordinateSets(bc,NumAtoms());
143
144
    //Set reference array
145
0
    unsigned char *ref = new unsigned char [NumRotors()*4];
146
0
    if (ref)
147
0
      {
148
0
        GetReferenceArray(ref);
149
0
        new_rml->Setup(*newmol,ref,NumRotors());
150
0
        delete [] ref;
151
0
      }
152
153
    //Set Rotamers
154
0
    unsigned char *rotamers = new unsigned char [(NumRotors()+1)*NumRotamers()];
155
0
    if (rotamers)
156
0
      {
157
0
        vector<unsigned char*>::const_iterator kk;
158
0
        unsigned int idx=0;
159
0
        for (kk = _vrotamer.begin();kk != _vrotamer.end();++kk)
160
0
          {
161
0
            memcpy(&rotamers[idx],(const unsigned char*)*kk,sizeof(unsigned char)*(NumRotors()+1));
162
0
            idx += sizeof(unsigned char)*(NumRotors()+1);
163
0
          }
164
0
        new_rml->AddRotamers(rotamers,NumRotamers());
165
0
        delete [] rotamers;
166
0
      }
167
0
    return new_rml;
168
0
  }
169
170
  OBRotamerList::~OBRotamerList()
171
0
  {
172
0
    vector<unsigned char*>::iterator i;
173
0
    for (i = _vrotamer.begin();i != _vrotamer.end();++i)
174
0
      delete [] *i;
175
176
0
    vector<pair<OBAtom**,vector<int> > >::iterator j;
177
0
    for (j = _vrotor.begin();j != _vrotor.end();++j)
178
0
      delete [] j->first;
179
180
    //Delete the interal base coordinate list
181
0
    unsigned int k;
182
0
    for (k=0 ; k<_c.size() ; ++k)
183
0
      delete [] _c[k];
184
0
  }
185
186
  void OBRotamerList::GetReferenceArray(unsigned char *ref)const
187
0
  {
188
0
    int j;
189
0
    vector<pair<OBAtom**,vector<int> > >::const_iterator i;
190
0
    for (j=0,i = _vrotor.begin();i != _vrotor.end();++i)
191
0
      {
192
0
        ref[j++] = (unsigned char)(i->first[0])->GetIdx();
193
0
        ref[j++] = (unsigned char)(i->first[1])->GetIdx();
194
0
        ref[j++] = (unsigned char)(i->first[2])->GetIdx();
195
0
        ref[j++] = (unsigned char)(i->first[3])->GetIdx();
196
0
      }
197
0
  }
198
199
  void OBRotamerList::Setup(OBMol &mol,OBRotorList &rl)
200
0
  {
201
    //clear the old stuff out if necessary
202
0
    _vres.clear();
203
0
    vector<unsigned char*>::iterator j;
204
0
    for (j = _vrotamer.begin();j != _vrotamer.end();++j)
205
0
      delete [] *j;
206
0
    _vrotamer.clear();
207
208
0
    vector<pair<OBAtom**,vector<int> > >::iterator k;
209
0
    for (k = _vrotor.begin();k != _vrotor.end();++k)
210
0
      delete [] k->first;
211
0
    _vrotor.clear();
212
0
    _vrings.clear();
213
0
    _vringTors.clear();
214
215
    //create the new list
216
0
    OBRotor *rotor;
217
0
    vector<OBRotor*>::iterator i;
218
0
    vector<int> children;
219
220
0
    int ref[4];
221
0
    OBAtom **atomlist;
222
0
    for (rotor = rl.BeginRotor(i);rotor;rotor = rl.NextRotor(i))
223
0
      {
224
0
        atomlist = new OBAtom* [4];
225
0
        rotor->GetDihedralAtoms(ref);
226
0
        atomlist[0] = mol.GetAtom(ref[0]);
227
0
        atomlist[1] = mol.GetAtom(ref[1]);
228
0
        atomlist[2] = mol.GetAtom(ref[2]);
229
0
        atomlist[3] = mol.GetAtom(ref[3]);
230
0
        mol.FindChildren(children,ref[1],ref[2]);
231
0
        _vrotor.push_back(pair<OBAtom**,vector<int> > (atomlist,children));
232
0
        _vres.push_back(rotor->GetResolution());
233
0
      }
234
235
    // if the rotor list has ring bonds, build up an index
236
0
    if (rl.HasRingRotors()){
237
      // go through rings
238
      // for each step of the path, see if there's a matching rotor
239
0
      vector<int> path;
240
0
      int pSize;
241
0
      vector<double> ringTorsions;
242
0
      vector<int> ringRotors;
243
0
      FOR_RINGS_OF_MOL(r, mol)
244
0
      {
245
0
        ringTorsions.clear();
246
0
        ringRotors.clear();
247
248
0
        pSize = r->Size();
249
0
        if (pSize < 4)
250
0
          continue; // not rotatable
251
252
0
        path = r->_path;
253
0
        for (int j = 0; j < pSize; ++j) {
254
0
          double torsion = mol.GetTorsion(path[(j + pSize - 1) % pSize],
255
0
                                       path[(j + pSize) % pSize],
256
0
                                       path[(j + pSize + 1) % pSize],
257
0
                                       path[(j + pSize + 2) % pSize]);
258
0
          ringTorsions.push_back(torsion);
259
260
          // now check to see if any of these things are rotors
261
0
          int rotorIndex = -1; // not a rotor
262
0
          OBBond *bond = mol.GetBond(path[(j + pSize) % pSize], path[(j + pSize + 1) % pSize]);
263
0
          for (rotor = rl.BeginRotor(i);rotor;rotor = rl.NextRotor(i))
264
0
            {
265
0
              if (bond != rotor->GetBond())
266
0
                continue; // no match at all
267
268
              // Central bond matches, make sure 1..4 atoms are in the path
269
0
              rotor->GetDihedralAtoms(ref);
270
0
              if ( (ref[0] == path[(j + pSize - 1) % pSize] &&
271
0
                    ref[3] == path[(j + pSize + 2) % pSize])
272
0
                   ||
273
0
                   (ref[3] == path[(j + pSize - 1) % pSize] &&
274
0
                    ref[0] == path[(j + pSize + 2) % pSize]) ) {
275
0
                rotorIndex = rotor->GetIdx();
276
0
              }
277
0
            } // end checking all the rotors
278
0
          ringRotors.push_back(rotorIndex); // could be -1 if it's not rotatable
279
0
        }
280
281
0
        _vringTors.push_back(ringTorsions);
282
0
        _vrings.push_back(ringRotors);
283
0
      } // finished with the rings
284
0
    } // if (ring rotors)
285
286
0
    vector<double>::iterator n;
287
0
    vector<vector<double> >::iterator m;
288
0
    for (m = _vres.begin();m != _vres.end();++m)
289
0
      for (n = m->begin();n != m->end();++n)
290
0
        *n *= RAD_TO_DEG;
291
0
  }
292
293
  void OBRotamerList::Setup(OBMol &mol,unsigned char *ref,int nrotors)
294
0
  {
295
    //clear the old stuff out if necessary
296
0
    _vres.clear();
297
0
    vector<unsigned char*>::iterator j;
298
0
    for (j = _vrotamer.begin();j != _vrotamer.end();++j)
299
0
      delete [] *j;
300
0
    _vrotamer.clear();
301
302
0
    vector<pair<OBAtom**,vector<int> > >::iterator k;
303
0
    for (k = _vrotor.begin();k != _vrotor.end();++k)
304
0
      delete [] k->first;
305
0
    _vrotor.clear();
306
0
    _vrings.clear();
307
0
    _vringTors.clear();
308
309
    //create the new list
310
0
    int i;
311
0
    vector<int> children;
312
313
0
    int refatoms[4];
314
0
    OBAtom **atomlist;
315
0
    for (i = 0; i < nrotors; ++i)
316
0
      {
317
0
        atomlist = new OBAtom* [4];
318
0
        refatoms[0] = (int)ref[i*4  ];
319
0
        refatoms[1] = (int)ref[i*4+1];
320
0
        refatoms[2] = (int)ref[i*4+2];
321
0
        refatoms[3] = (int)ref[i*4+3];
322
0
        mol.FindChildren(children,refatoms[1],refatoms[2]);
323
0
        atomlist[0] = mol.GetAtom(refatoms[0]);
324
0
        atomlist[1] = mol.GetAtom(refatoms[1]);
325
0
        atomlist[2] = mol.GetAtom(refatoms[2]);
326
0
        atomlist[3] = mol.GetAtom(refatoms[3]);
327
0
        _vrotor.push_back(pair<OBAtom**,vector<int> > (atomlist,children));
328
0
      }
329
330
0
  }
331
332
  void OBRotamerList::AddRotamer(double *c)
333
0
  {
334
    // TODO: consider implementing periodic boundary conditions
335
0
    int idx,size;
336
0
    double angle,res=255.0/360.0;
337
0
    vector3 v1,v2,v3,v4;
338
339
0
    unsigned char *rot = new unsigned char [_vrotor.size()+1];
340
0
    rot[0] = (char) 0;
341
342
0
    vector<pair<OBAtom**,vector<int> > >::iterator i;
343
0
    for (size=1,i = _vrotor.begin();i != _vrotor.end();++i,++size)
344
0
      {
345
0
        idx = (i->first[0])->GetCoordinateIdx();
346
0
        v1.Set(c[idx],c[idx+1],c[idx+2]);
347
0
        idx = (i->first[1])->GetCoordinateIdx();
348
0
        v2.Set(c[idx],c[idx+1],c[idx+2]);
349
0
        idx = (i->first[2])->GetCoordinateIdx();
350
0
        v3.Set(c[idx],c[idx+1],c[idx+2]);
351
0
        idx = (i->first[3])->GetCoordinateIdx();
352
0
        v4.Set(c[idx],c[idx+1],c[idx+2]);
353
354
0
        angle = CalcTorsionAngle(v1,v2,v3,v4);
355
0
        while (angle < 0.0)
356
0
          angle += 360.0;
357
0
        while (angle > 360.0)
358
0
          angle -= 360.0;
359
0
        rot[size] = (unsigned char)rint(angle*res);
360
0
      }
361
362
0
    _vrotamer.push_back(rot);
363
0
  }
364
365
  void OBRotamerList::AddRotamer(int *arr)
366
0
  {
367
0
    unsigned int i;
368
0
    double angle,res=255.0/360.0;
369
370
0
    unsigned char *rot = new unsigned char [_vrotor.size()+1];
371
0
    rot[0] = (unsigned char)arr[0];
372
373
0
    for (i = 0;i < _vrotor.size();++i)
374
0
      {
375
0
        angle = _vres[i][arr[i+1]];
376
0
        while (angle < 0.0)
377
0
          angle += 360.0;
378
0
        while (angle > 360.0)
379
0
          angle -= 360.0;
380
0
        rot[i+1] = (unsigned char)rint(angle*res);
381
0
      }
382
0
    _vrotamer.push_back(rot);
383
0
  }
384
385
  void OBRotamerList::AddRotamer(std::vector<int> arr)
386
0
  {
387
0
    unsigned int i;
388
0
    double angle,res=255.0/360.0;
389
390
0
    if (arr.size() != (_vrotor.size() + 1))
391
0
      return; // wrong size key
392
393
    // gotta check for weird ring torsion combinations
394
0
    if (_vrings.size()) {
395
      // go through each ring and update the possible torsions
396
0
      for (unsigned int j = 0; j < _vrings.size(); ++j) {
397
0
        vector<int> path = _vrings[j];
398
0
        double torsionSum = 0.0;
399
400
        // go around the loop and add up the torsions
401
0
        for (unsigned int i = 0; i < path.size(); ++i) {
402
0
          if (path[i] == -1) {
403
            // not a rotor, use the fixed value
404
0
            torsionSum += _vringTors[j][i];
405
0
            continue;
406
0
          }
407
408
          // what angle are we trying to use with this key?
409
0
          angle = _vres[ path[i] ][arr[ path[i]+1 ]]*res;
410
0
          while (angle < 0.0)
411
0
            angle += 360.0;
412
0
          while (angle > 360.0)
413
0
            angle -= 360.0;
414
415
          // update the ring torsion for this setting
416
0
          _vringTors[j][i] = angle;
417
0
          torsionSum += angle;
418
0
        }
419
420
        // if the sum of the ring torsions is not ~0, bad move
421
0
        if (fabs(torsionSum) > 45.0) {
422
          //          cerr << " Bad move! " << fabs(torsionSum) << endl;
423
0
          return; // don't make the move
424
0
        }
425
        //        cerr << " Good move!" << endl;
426
0
      }
427
0
    }
428
429
0
    unsigned char *rot = new unsigned char [_vrotor.size()+1];
430
0
    rot[0] = (unsigned char)arr[0];
431
432
0
    for (i = 0;i < _vrotor.size();++i)
433
0
      {
434
0
        angle = _vres[i][arr[i+1]];
435
0
        while (angle < 0.0)
436
0
          angle += 360.0;
437
0
        while (angle > 360.0)
438
0
          angle -= 360.0;
439
0
        rot[i+1] = (unsigned char)rint(angle*res);
440
0
      }
441
0
    _vrotamer.push_back(rot);
442
0
  }
443
444
  void OBRotamerList::AddRotamer(unsigned char *arr)
445
0
  {
446
0
    unsigned int i;
447
0
    double angle,res=255.0/360.0;
448
449
0
    unsigned char *rot = new unsigned char [_vrotor.size()+1];
450
0
    rot[0] = (unsigned char)arr[0];
451
452
0
    for (i = 0;i < _vrotor.size();++i)
453
0
      {
454
0
        angle = _vres[i][(int)arr[i+1]];
455
0
        while (angle < 0.0)
456
0
          angle += 360.0;
457
0
        while (angle > 360.0)
458
0
          angle -= 360.0;
459
0
        rot[i+1] = (unsigned char)rint(angle*res);
460
0
      }
461
0
    _vrotamer.push_back(rot);
462
0
  }
463
464
  void OBRotamerList::AddRotamers(unsigned char *arr,int nrotamers)
465
0
  {
466
0
    unsigned int size;
467
0
    int i;
468
469
0
    size = (unsigned int)_vrotor.size()+1;
470
0
    for (i = 0;i < nrotamers;++i)
471
0
      {
472
0
        unsigned char *rot = new unsigned char [size];
473
0
        memcpy(rot,&arr[i*size],sizeof(char)*size);
474
0
        _vrotamer.push_back(rot);
475
0
      }
476
0
  }
477
478
  void OBRotamerList::ExpandConformerList(OBMol &mol,vector<double*> &clist)
479
0
  {
480
0
    vector<double*> tmpclist = CreateConformerList(mol);
481
482
    //transfer the conf list
483
0
    vector<double*>::iterator k;
484
0
    for (k = clist.begin();k != clist.end();++k)
485
0
      delete [] *k;
486
0
    clist = tmpclist;
487
0
  }
488
489
  //! Create a conformer list using the internal base set of coordinates
490
  vector<double*> OBRotamerList::CreateConformerList(OBMol& mol)
491
0
  {
492
0
    unsigned int j;
493
0
    double angle,invres=360.0/255.0;
494
0
    unsigned char *conf;
495
0
    vector<double*> tmpclist;
496
0
    vector<unsigned char*>::iterator i;
497
498
0
    for (i = _vrotamer.begin();i != _vrotamer.end();++i)
499
0
      {
500
0
        conf = *i;
501
0
        double *c = new double [mol.NumAtoms()*3];
502
0
        memcpy(c,_c[(int)conf[0]],sizeof(double)*mol.NumAtoms()*3);
503
504
0
        for (j = 0;j < _vrotor.size();++j)
505
0
          {
506
0
            angle = invres*((double)conf[j+1]);
507
0
            if (angle > 180.0)
508
0
              angle -= 360.0;
509
0
            SetRotorToAngle(c,_vrotor[j].first,angle,_vrotor[j].second);
510
0
          }
511
0
        tmpclist.push_back(c);
512
0
      }
513
514
0
    return tmpclist;
515
0
  }
516
517
  //! Change the current coordinate set
518
  //! \since version 2.2
519
  bool OBRotamerList::SetCurrentCoordinates(OBMol &mol, std::vector<int> arr)
520
0
  {
521
0
    double angle;
522
523
0
    if (arr.size() != (_vrotor.size() + 1))
524
0
      return false; // wrong size key
525
526
    // gotta check for weird ring torsion combinations
527
0
    if (_vrings.size()) {
528
      // go through each ring and update the possible torsions
529
0
      for (unsigned int j = 0; j < _vrings.size(); ++j) {
530
0
        vector<int> path = _vrings[j];
531
0
        double torsionSum = 0.0;
532
533
        // go around the loop and add up the torsions
534
0
        for (unsigned int i = 0; i < path.size(); ++i) {
535
0
          if (path[i] == -1) {
536
            // not a rotor, use the fixed value
537
0
            torsionSum += _vringTors[j][i];
538
0
            continue;
539
0
          }
540
541
          // what angle are we trying to use with this key?
542
0
          angle = _vres[ path[i] ][arr[ path[i]+1 ]];
543
0
          while (angle < 0.0)
544
0
            angle += 360.0;
545
0
          while (angle > 360.0)
546
0
            angle -= 360.0;
547
548
          // update the ring torsion for this setting
549
0
          _vringTors[j][i] = angle;
550
0
          torsionSum += angle;
551
0
        }
552
553
        // if the sum of the ring torsions is not ~0, bad move
554
0
        if (fabs(torsionSum) > 45.0) {
555
          //          cerr << " Bad move!" << endl;
556
0
          return false; // don't make the move
557
0
        }
558
0
      }
559
0
    }
560
561
0
    double *c = mol.GetCoordinates();
562
0
    for (unsigned int i = 0;i < _vrotor.size();++i)
563
0
      {
564
0
        if (arr[i+1] == -1) // skip this rotor
565
0
          continue;
566
0
        else {
567
0
          angle = _vres[i][arr[i+1]];
568
0
          while (angle < 0.0)
569
0
            angle += 360.0;
570
0
          while (angle > 360.0)
571
0
            angle -= 360.0;
572
0
          SetRotorToAngle(c,_vrotor[i].first,angle,_vrotor[i].second);
573
0
        } // set an angle
574
0
      } // for rotors
575
576
0
    return true;
577
0
  }
578
579
  void OBRotamerList::SetBaseCoordinateSets(OBMol& mol)
580
0
  {
581
0
    SetBaseCoordinateSets(mol.GetConformers(), mol.NumAtoms());
582
0
  }
583
584
  //Copies the coordinates in bc, NOT the pointers, into the object
585
  void OBRotamerList::SetBaseCoordinateSets(vector<double*> bc, unsigned int N)
586
0
  {
587
0
    unsigned int i,j;
588
589
    //Clear out old data
590
0
    for (i=0 ; i<_c.size() ; ++i)
591
0
      delete [] _c[i];
592
0
    _c.clear();
593
594
    //Copy new data
595
0
    double *c  = nullptr;
596
0
    double *cc = nullptr;
597
0
    for (i=0 ; i<bc.size() ; ++i)
598
0
      {
599
0
        c = new double [3*N];
600
0
        cc = bc[i];
601
0
        for (j=0 ; j<3*N ; ++j)
602
0
          c[j] = cc[j];
603
0
        _c.push_back(c);
604
0
      }
605
0
    _NBaseCoords = N;
606
0
  }
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
613
  void SetRotorToAngle(double *c, OBAtom **ref,double ang,vector<int> atoms)
614
0
  {
615
0
    double v1x,v1y,v1z,v2x,v2y,v2z,v3x,v3y,v3z;
616
0
    double c1x,c1y,c1z,c2x,c2y,c2z,c3x,c3y,c3z;
617
0
    double c1mag,c2mag,radang,costheta,m[9];
618
0
    double x,y,z,mag,rotang,sn,cs,t,tx,ty,tz;
619
620
0
    int tor[4];
621
0
    tor[0] = ref[0]->GetCoordinateIdx();
622
0
    tor[1] = ref[1]->GetCoordinateIdx();
623
0
    tor[2] = ref[2]->GetCoordinateIdx();
624
0
    tor[3] = ref[3]->GetCoordinateIdx();
625
626
    //
627
    //calculate the torsion angle
628
    //
629
0
    v1x = c[tor[0]]   - c[tor[1]];   v2x = c[tor[1]]   - c[tor[2]];
630
0
    v1y = c[tor[0]+1] - c[tor[1]+1]; v2y = c[tor[1]+1] - c[tor[2]+1];
631
0
    v1z = c[tor[0]+2] - c[tor[1]+2]; v2z = c[tor[1]+2] - c[tor[2]+2];
632
0
    v3x = c[tor[2]]   - c[tor[3]];
633
0
    v3y = c[tor[2]+1] - c[tor[3]+1];
634
0
    v3z = c[tor[2]+2] - c[tor[3]+2];
635
636
0
    c1x = v1y*v2z - v1z*v2y;   c2x = v2y*v3z - v2z*v3y;
637
0
    c1y = -v1x*v2z + v1z*v2x;  c2y = -v2x*v3z + v2z*v3x;
638
0
    c1z = v1x*v2y - v1y*v2x;   c2z = v2x*v3y - v2y*v3x;
639
0
    c3x = c1y*c2z - c1z*c2y;
640
0
    c3y = -c1x*c2z + c1z*c2x;
641
0
    c3z = c1x*c2y - c1y*c2x;
642
643
0
    c1mag = c1x*c1x + c1y*c1y + c1z*c1z;
644
0
    c2mag = c2x*c2x + c2y*c2y + c2z*c2z;
645
0
    if (c1mag*c2mag < 0.01) costheta = 1.0; //avoid div by zero error
646
0
    else costheta = (c1x*c2x + c1y*c2y + c1z*c2z)/(sqrt(c1mag*c2mag));
647
648
0
    if (costheta < -0.999999) costheta = -0.999999;
649
0
    if (costheta >  0.999999) costheta =  0.999999;
650
651
0
    if ((v2x*c3x + v2y*c3y + v2z*c3z) > 0.0) radang = -acos(costheta);
652
0
    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
0
    rotang = (DEG_TO_RAD*ang) - radang;
660
661
0
    sn = sin(rotang); cs = cos(rotang);t = 1 - cs;
662
    //normalize the rotation vector
663
0
    mag = sqrt(v2x*v2x + v2y*v2y + v2z*v2z);
664
0
    if (mag < 0.1) mag = 0.1; // avoid divide by zero error
665
0
    x = v2x/mag; y = v2y/mag; z = v2z/mag;
666
667
    //set up the rotation matrix
668
0
    m[0]= t*x*x + cs;     m[1] = t*x*y + sn*z;  m[2] = t*x*z - sn*y;
669
0
    m[3] = t*x*y - sn*z;  m[4] = t*y*y + cs;    m[5] = t*y*z + sn*x;
670
0
    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
0
    tx = c[tor[1]];ty = c[tor[1]+1];tz = c[tor[1]+2];
676
0
    vector<int>::iterator i;int j;
677
0
    for (i = atoms.begin();i != atoms.end();++i)
678
0
      {
679
0
        j = ((*i)-1)*3;
680
0
        c[j] -= tx;c[j+1] -= ty;c[j+2]-= tz;
681
0
        x = c[j]*m[0] + c[j+1]*m[1] + c[j+2]*m[2];
682
0
        y = c[j]*m[3] + c[j+1]*m[4] + c[j+2]*m[5];
683
0
        z = c[j]*m[6] + c[j+1]*m[7] + c[j+2]*m[8];
684
0
        c[j] = x; c[j+1] = y; c[j+2] = z;
685
0
        c[j] += tx;c[j+1] += ty;c[j+2] += tz;
686
0
      }
687
0
  }
688
689
  int PackCoordinate(double c[3],double max[3])
690
0
  {
691
0
    int tmp;
692
0
    double cf;
693
0
    cf = c[0];
694
0
    tmp  = ((int)(cf*max[0])) << 20;
695
0
    cf = c[1];
696
0
    tmp |= ((int)(cf*max[1])) << 10;
697
0
    cf = c[2];
698
0
    tmp |= ((int)(cf*max[2]));
699
0
    return(tmp);
700
0
  }
701
702
  void UnpackCoordinate(double c[3],double max[3],int tmp)
703
0
  {
704
0
    double cf;
705
0
    cf = (double)(tmp>>20);
706
0
    c[0] = cf;
707
0
    c[0] *= max[0];
708
0
    cf = (double)((tmp&0xffc00)>>10);
709
0
    c[1] = cf;
710
0
    c[1] *= max[1];
711
0
    cf = (double)(tmp&0x3ff);
712
0
    c[2] = cf;
713
0
    c[2] *= max[2];
714
0
  }
715
716
} //namespace OpenBabel
717
718
//! \file rotamer.cpp
719
//! \brief Handle rotamer list data.