← Back to branch summary

~ubuntu-branches/debian/stretch/openbabel/stretch

« back to all changes in this revision

Viewing changes to src/rotamer.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2008-07-22 23:54:58 UTC
  • mfrom: (3.1.10 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080722235458-3o606czluviz4akx
Tags: 2.2.0-2
* Upload to unstable.
* debian/control: Updated descriptions.
* debian/patches/gauss_cube_format.patch: New patch, makes the 
  gaussian cube format available again.
* debian/rules (DEB_DH_MAKESHLIBS_ARGS_libopenbabel3): Removed.
* debian/rules (DEB_CONFIGURE_EXTRA_FLAGS): Likewise.
* debian/libopenbabel3.install: Adjust formats directory to single 
  version hierarchy.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/rotamer.cpp
2
2
rotamer.cpp - Handle rotamer list data.
3
3
 
4
4
Copyright (C) 1998, 1999, 2000-2002 OpenEye Scientific Software, Inc.
5
 
Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
 
5
Some portions Copyright (C) 2001-2006 by Geoffrey R. Hutchison
6
6
 
7
7
This file is part of the Open Babel project.
8
8
For more information, see <http://openbabel.sourceforge.net/>
16
16
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17
17
GNU General Public License for more details.
18
18
***********************************************************************/
 
19
#include <openbabel/babelconfig.h>
19
20
 
20
 
#include "rotamer.h"
21
 
#include "mol.h"
22
 
#include "obutil.h"
23
 
#include "rotor.h"
 
21
#include <openbabel/rotamer.h>
24
22
 
25
23
#define OB_TITLE_SIZE     254
26
24
#define OB_BINARY_SETWORD 32
30
28
namespace OpenBabel
31
29
{
32
30
 
33
 
//test byte ordering
34
 
static int SINT = 0x00000001;
35
 
static unsigned char *STPTR = (unsigned char*)&SINT;
36
 
const bool SwabInt = (STPTR[0]!=0);
 
31
  /** \class OBRotamerList rotamer.h <openbabel/rotamer.h>
 
32
 
 
33
      A high-level class for rotamer / conformer generation, intended mainly
 
34
      for use with the related class OBRotorList and the OBRotorRules database
 
35
 
 
36
      Rotamers represent conformational isomers formed simply by rotation of
 
37
      dihedral angles. On the other hand, conformers may include geometric
 
38
      relaxation (i.e., slight modification of bond lengths, bond angles, etc.)
 
39
 
 
40
      The following shows an example of generating 2 conformers using different
 
41
      rotor states. Similar code could be used for systematic or Monte Carlo
 
42
      conformer sampling when combined with energy evaluation (molecular
 
43
      mechanics or otherwise).
 
44
 
 
45
      \code
 
46
      OBRotorList rl; // used to sample all rotatable bonds via the OBRotorRules data
 
47
      // If you want to "fix" any particular atoms (i.e., freeze them in space)
 
48
      // then set up an OBBitVec of the fixed atoms and call
 
49
      // rl.SetFixAtoms(bitvec);
 
50
      rl.Setup(mol);
 
51
 
 
52
      // How many rotatable bonds are there?
 
53
      cerr << " Number of rotors: " << rl.Size() << endl;
 
54
 
 
55
      // indexed from 1, rotorKey[0] = 0
 
56
      std::vector<int> rotorKey(rl.Size() + 1, 0);
 
57
 
 
58
      // each entry represents the configuration of a rotor
 
59
      // e.g. indexes into OBRotor::GetResolution() -- the different angles
 
60
      //   to sample for a rotamer search
 
61
      for (unsigned int i = 0; i < rl.Size() + 1; ++i)
 
62
      rotorKey[i] = 0; // could be anything from 0 .. OBRotor->GetResolution().size()
 
63
      // -1 is for no rotation
 
64
 
 
65
      // The OBRotamerList can generate conformations (i.e., coordinate sets)
 
66
      OBRotamerList rotamers;
 
67
      rotamers.SetBaseCoordinateSets(mol);
 
68
      rotamers.Setup(mol, rl);
 
69
 
 
70
      rotamers.AddRotamer(rotorKey);
 
71
      rotorKey[1] = 2; // switch one rotor
 
72
      rotamers.AddRotamer(rotorKey);
 
73
 
 
74
      rotamers.ExpandConformerList(mol, mol.GetConformers());
 
75
 
 
76
      // change the molecule conformation
 
77
      mol.SetConformer(0); // rotorKey 0, 0, ...
 
78
      conv.Write(&mol);
 
79
 
 
80
      mol.SetConformer(1); // rotorKey 0, 2, ...
 
81
 
 
82
      \endcode
 
83
 
 
84
  **/
 
85
 
 
86
  //test byte ordering
 
87
  static int SINT = 0x00000001;
 
88
  static unsigned char *STPTR = (unsigned char*)&SINT;
 
89
  const bool SwabInt = (STPTR[0]!=0);
37
90
 
38
91
#if !HAVE_RINT
39
 
inline double rint(double x)
40
 
{
 
92
  inline double rint(double x)
 
93
  {
41
94
    return ( (x < 0.0) ? ceil(x-0.5) : floor(x+0.5));
42
 
}
 
95
  }
43
96
#endif
44
97
 
45
 
void SetRotorToAngle(double *c,OBAtom **ref,double ang,vector<int> atoms);
 
98
  void SetRotorToAngle(double *c,OBAtom **ref,double ang,vector<int> atoms);
46
99
 
47
 
int Swab(int i)
48
 
{
 
100
  int Swab(int i)
 
101
  {
49
102
    unsigned char tmp[4],c;
50
103
    memcpy(tmp,(char*)&i,sizeof(int));
51
104
    c = tmp[0];
56
109
    tmp[2] = c;
57
110
    memcpy((char*)&i,tmp,sizeof(int));
58
111
    return(i);
59
 
}
60
 
 
61
 
OBRotamerList::~OBRotamerList()
62
 
{
 
112
  }
 
113
 
 
114
  OBGenericData* OBRotamerList::Clone(OBBase* newparent) const
 
115
  {
 
116
    //Since the class contains OBAtom pointers, the new copy use
 
117
    //these from the new molecule, newparent
 
118
    OBMol* newmol = static_cast<OBMol*>(newparent);
 
119
        
 
120
    OBRotamerList *new_rml = new OBRotamerList;
 
121
    new_rml->_attr = _attr;
 
122
    new_rml->_type = _type;
 
123
 
 
124
    //Set base coordinates
 
125
    unsigned int k,l;
 
126
    vector<double*> bc;
 
127
    double *c=NULL;
 
128
    double *cc=NULL;
 
129
    for (k=0 ; k<NumBaseCoordinateSets() ; ++k)
 
130
      {
 
131
        c = new double [3*NumAtoms()];
 
132
        cc = GetBaseCoordinateSet(k);
 
133
        for (l=0 ; l<3*NumAtoms() ; ++l)
 
134
                                        c[l] = cc[l];
 
135
        bc.push_back(c);
 
136
      }
 
137
    if (NumBaseCoordinateSets())
 
138
                        new_rml->SetBaseCoordinateSets(bc,NumAtoms());
 
139
 
 
140
    //Set reference array
 
141
    unsigned char *ref = new unsigned char [NumRotors()*4];
 
142
    if (ref)
 
143
      {
 
144
        GetReferenceArray(ref);
 
145
        new_rml->Setup(*newmol,ref,NumRotors()); 
 
146
        delete [] ref;
 
147
      }
 
148
 
 
149
    //Set Rotamers
 
150
    unsigned char *rotamers = new unsigned char [(NumRotors()+1)*NumRotamers()];
 
151
    if (rotamers)
 
152
      {
 
153
        vector<unsigned char*>::const_iterator kk;
 
154
        unsigned int idx=0;
 
155
        for (kk = _vrotamer.begin();kk != _vrotamer.end();++kk)
 
156
          {
 
157
            memcpy(&rotamers[idx],(const unsigned char*)*kk,sizeof(unsigned char)*(NumRotors()+1));
 
158
            idx += sizeof(unsigned char)*(NumRotors()+1);
 
159
          }
 
160
        new_rml->AddRotamers(rotamers,NumRotamers());
 
161
        delete [] rotamers;
 
162
      }
 
163
    return new_rml;
 
164
  }
 
165
 
 
166
  OBRotamerList::~OBRotamerList()
 
167
  {
63
168
    vector<unsigned char*>::iterator i;
64
 
    for (i = _vrotamer.begin();i != _vrotamer.end();i++)
65
 
        delete [] *i;
 
169
    for (i = _vrotamer.begin();i != _vrotamer.end();++i)
 
170
      delete [] *i;
66
171
 
67
172
    vector<pair<OBAtom**,vector<int> > >::iterator j;
68
 
    for (j = _vrotor.begin();j != _vrotor.end();j++)
69
 
        delete [] j->first;
 
173
    for (j = _vrotor.begin();j != _vrotor.end();++j)
 
174
      delete [] j->first;
70
175
 
71
176
    //Delete the interal base coordinate list
72
177
    unsigned int k;
73
 
    for (k=0 ; k<_c.size() ; k++)
74
 
        delete [] _c[k];
75
 
}
 
178
    for (k=0 ; k<_c.size() ; ++k)
 
179
      delete [] _c[k];
 
180
  }
76
181
 
77
 
void OBRotamerList::GetReferenceArray(unsigned char *ref)
78
 
{
 
182
  void OBRotamerList::GetReferenceArray(unsigned char *ref)const
 
183
  {
79
184
    int j;
80
 
    vector<pair<OBAtom**,vector<int> > >::iterator i;
81
 
    for (j=0,i = _vrotor.begin();i != _vrotor.end();i++)
82
 
    {
 
185
                vector<pair<OBAtom**,vector<int> > >::const_iterator i;
 
186
    for (j=0,i = _vrotor.begin();i != _vrotor.end();++i)
 
187
      {
83
188
        ref[j++] = (unsigned char)(i->first[0])->GetIdx();
84
189
        ref[j++] = (unsigned char)(i->first[1])->GetIdx();
85
190
        ref[j++] = (unsigned char)(i->first[2])->GetIdx();
86
191
        ref[j++] = (unsigned char)(i->first[3])->GetIdx();
87
 
    }
88
 
}
 
192
      }
 
193
  }
89
194
 
90
 
void OBRotamerList::Setup(OBMol &mol,OBRotorList &rl)
91
 
{
 
195
  void OBRotamerList::Setup(OBMol &mol,OBRotorList &rl)
 
196
  {
92
197
    //clear the old stuff out if necessary
93
198
    _vres.clear();
94
199
    vector<unsigned char*>::iterator j;
95
 
    for (j = _vrotamer.begin();j != _vrotamer.end();j++)
96
 
        delete [] *j;
 
200
    for (j = _vrotamer.begin();j != _vrotamer.end();++j)
 
201
      delete [] *j;
97
202
    _vrotamer.clear();
98
203
 
99
204
    vector<pair<OBAtom**,vector<int> > >::iterator k;
100
 
    for (k = _vrotor.begin();k != _vrotor.end();k++)
101
 
        delete [] k->first;
 
205
    for (k = _vrotor.begin();k != _vrotor.end();++k)
 
206
      delete [] k->first;
102
207
    _vrotor.clear();
103
208
 
104
209
    //create the new list
109
214
    int ref[4];
110
215
    OBAtom **atomlist;
111
216
    for (rotor = rl.BeginRotor(i);rotor;rotor = rl.NextRotor(i))
112
 
    {
 
217
      {
113
218
        atomlist = new OBAtom* [4];
114
219
        rotor->GetDihedralAtoms(ref);
115
220
        atomlist[0] = mol.GetAtom(ref[0]);
119
224
        mol.FindChildren(children,ref[1],ref[2]);
120
225
        _vrotor.push_back(pair<OBAtom**,vector<int> > (atomlist,children));
121
226
        _vres.push_back(rotor->GetResolution());
122
 
    }
 
227
      }
123
228
 
124
229
    vector<double>::iterator n;
125
230
    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
 
}
 
231
    for (m = _vres.begin();m != _vres.end();++m)
 
232
      for (n = m->begin();n != m->end();++n)
 
233
        *n *= RAD_TO_DEG;
 
234
  }
130
235
 
131
 
void OBRotamerList::Setup(OBMol &mol,unsigned char *ref,int nrotors)
132
 
{
 
236
  void OBRotamerList::Setup(OBMol &mol,unsigned char *ref,int nrotors)
 
237
  {
133
238
    //clear the old stuff out if necessary
134
239
    _vres.clear();
135
240
    vector<unsigned char*>::iterator j;
136
 
    for (j = _vrotamer.begin();j != _vrotamer.end();j++)
137
 
        delete [] *j;
 
241
    for (j = _vrotamer.begin();j != _vrotamer.end();++j)
 
242
      delete [] *j;
138
243
    _vrotamer.clear();
139
244
 
140
245
    vector<pair<OBAtom**,vector<int> > >::iterator k;
141
 
    for (k = _vrotor.begin();k != _vrotor.end();k++)
142
 
        delete [] k->first;
 
246
    for (k = _vrotor.begin();k != _vrotor.end();++k)
 
247
      delete [] k->first;
143
248
    _vrotor.clear();
144
249
 
145
250
    //create the new list
148
253
 
149
254
    int refatoms[4];
150
255
    OBAtom **atomlist;
151
 
    for (i = 0; i < nrotors; i++)
152
 
    {
 
256
    for (i = 0; i < nrotors; ++i)
 
257
      {
153
258
        atomlist = new OBAtom* [4];
154
259
        refatoms[0] = (int)ref[i*4  ];
155
260
        refatoms[1] = (int)ref[i*4+1];
161
266
        atomlist[2] = mol.GetAtom(refatoms[2]);
162
267
        atomlist[3] = mol.GetAtom(refatoms[3]);
163
268
        _vrotor.push_back(pair<OBAtom**,vector<int> > (atomlist,children));
164
 
    }
165
 
 
166
 
}
167
 
 
168
 
void OBRotamerList::AddRotamer(double *c)
169
 
{
 
269
      }
 
270
 
 
271
  }
 
272
 
 
273
  void OBRotamerList::AddRotamer(double *c)
 
274
  {
170
275
    int idx,size;
171
 
    double angle,res=255.0f/360.0f;
 
276
    double angle,res=255.0/360.0;
172
277
    vector3 v1,v2,v3,v4;
173
278
 
174
279
    unsigned char *rot = new unsigned char [_vrotor.size()+1];
175
280
    rot[0] = (char) 0;
176
281
 
177
282
    vector<pair<OBAtom**,vector<int> > >::iterator i;
178
 
    for (size=1,i = _vrotor.begin();i != _vrotor.end();i++,size++)
179
 
    {
 
283
    for (size=1,i = _vrotor.begin();i != _vrotor.end();++i,++size)
 
284
      {
180
285
        idx = (i->first[0])->GetCIdx();
181
286
        v1.Set(c[idx],c[idx+1],c[idx+2]);
182
287
        idx = (i->first[1])->GetCIdx();
187
292
        v4.Set(c[idx],c[idx+1],c[idx+2]);
188
293
 
189
294
        angle = CalcTorsionAngle(v1,v2,v3,v4);
190
 
        while (angle < 0.0f)
191
 
            angle += 360.0f;
192
 
        while (angle > 360.0f)
193
 
            angle -= 360.0f;
 
295
        while (angle < 0.0)
 
296
          angle += 360.0;
 
297
        while (angle > 360.0)
 
298
          angle -= 360.0;
194
299
        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
 
    {
 
300
      }
 
301
 
 
302
    _vrotamer.push_back(rot);
 
303
  }
 
304
 
 
305
  void OBRotamerList::AddRotamer(int *arr)
 
306
  {
 
307
    unsigned int i;
 
308
    double angle,res=255.0/360.0;
 
309
 
 
310
    unsigned char *rot = new unsigned char [_vrotor.size()+1];
 
311
    rot[0] = (unsigned char)arr[0];
 
312
 
 
313
    for (i = 0;i < _vrotor.size();++i)
 
314
      {
 
315
        angle = _vres[i][arr[i+1]];
 
316
        while (angle < 0.0)
 
317
          angle += 360.0;
 
318
        while (angle > 360.0)
 
319
          angle -= 360.0;
 
320
        rot[i+1] = (unsigned char)rint(angle*res);
 
321
      }
 
322
    _vrotamer.push_back(rot);
 
323
  }
 
324
 
 
325
  void OBRotamerList::AddRotamer(std::vector<int> arr)
 
326
  {
 
327
    unsigned int i;
 
328
    double angle,res=255.0/360.0;
 
329
    
 
330
    if (arr.size() != (_vrotor.size() + 1))
 
331
      return; // wrong size key
 
332
 
 
333
    unsigned char *rot = new unsigned char [_vrotor.size()+1];
 
334
    rot[0] = (unsigned char)arr[0];
 
335
 
 
336
    for (i = 0;i < _vrotor.size();++i)
 
337
      {
 
338
        angle = _vres[i][arr[i+1]];
 
339
        while (angle < 0.0)
 
340
          angle += 360.0;
 
341
        while (angle > 360.0)
 
342
          angle -= 360.0;
 
343
        rot[i+1] = (unsigned char)rint(angle*res);
 
344
      }
 
345
    _vrotamer.push_back(rot);
 
346
  }
 
347
 
 
348
  void OBRotamerList::AddRotamer(unsigned char *arr)
 
349
  {
 
350
    unsigned int i;
 
351
    double angle,res=255.0/360.0;
 
352
 
 
353
    unsigned char *rot = new unsigned char [_vrotor.size()+1];
 
354
    rot[0] = (unsigned char)arr[0];
 
355
 
 
356
    for (i = 0;i < _vrotor.size();++i)
 
357
      {
230
358
        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;
 
359
        while (angle < 0.0)
 
360
          angle += 360.0;
 
361
        while (angle > 360.0)
 
362
          angle -= 360.0;
235
363
        rot[i+1] = (unsigned char)rint(angle*res);
236
 
    }
 
364
      }
237
365
    _vrotamer.push_back(rot);
238
 
}
 
366
  }
239
367
 
240
 
void OBRotamerList::AddRotamers(unsigned char *arr,int nrotamers)
241
 
{
 
368
  void OBRotamerList::AddRotamers(unsigned char *arr,int nrotamers)
 
369
  {
242
370
    unsigned int size;
243
371
    int i;
244
372
 
245
373
    size = (unsigned int)_vrotor.size()+1;
246
 
    for (i = 0;i < nrotamers;i++)
247
 
    {
 
374
    for (i = 0;i < nrotamers;++i)
 
375
      {
248
376
        unsigned char *rot = new unsigned char [size];
249
377
        memcpy(rot,&arr[i*size],sizeof(char)*size);
250
378
        _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
 
    }
 
379
      }
 
380
  }
 
381
 
 
382
  void OBRotamerList::ExpandConformerList(OBMol &mol,vector<double*> &clist)
 
383
  {
 
384
    vector<double*> tmpclist = CreateConformerList(mol);
277
385
 
278
386
    //transfer the conf list
279
387
    vector<double*>::iterator k;
280
 
    for (k = clist.begin();k != clist.end();k++)
281
 
        delete [] *k;
 
388
    for (k = clist.begin();k != clist.end();++k)
 
389
      delete [] *k;
282
390
    clist = tmpclist;
283
 
}
 
391
  }
284
392
 
285
 
//! Create a conformer list using the internal base set of coordinates
286
 
vector<double*> OBRotamerList::CreateConformerList(OBMol& mol)
287
 
{
 
393
  //! Create a conformer list using the internal base set of coordinates
 
394
  vector<double*> OBRotamerList::CreateConformerList(OBMol& mol)
 
395
  {
288
396
    unsigned int j;
289
 
    double angle,invres=360.0f/255.0f;
 
397
    double angle,invres=360.0/255.0;
290
398
    unsigned char *conf;
291
399
    vector<double*> tmpclist;
292
400
    vector<unsigned char*>::iterator i;
293
401
 
294
 
    for (i = _vrotamer.begin();i != _vrotamer.end();i++)
295
 
    {
 
402
    for (i = _vrotamer.begin();i != _vrotamer.end();++i)
 
403
      {
296
404
        conf = *i;
297
405
        double *c = new double [mol.NumAtoms()*3];
298
406
        memcpy(c,_c[(int)conf[0]],sizeof(double)*mol.NumAtoms()*3);
299
407
 
300
 
        for (j = 0;j < _vrotor.size();j++)
301
 
        {
 
408
        for (j = 0;j < _vrotor.size();++j)
 
409
          {
302
410
            angle = invres*((double)conf[j+1]);
303
411
            if (angle > 180.0)
304
 
                angle -= 360.0;
 
412
              angle -= 360.0;
305
413
            SetRotorToAngle(c,_vrotor[j].first,angle,_vrotor[j].second);
306
 
        }
 
414
          }
307
415
        tmpclist.push_back(c);
308
 
    }
 
416
      }
309
417
 
310
418
    return tmpclist;
311
 
}
 
419
  }
 
420
  
 
421
  //! Change the current coordinate set
 
422
  //! \since version 2.2
 
423
  void OBRotamerList::SetCurrentCoordinates(OBMol &mol, std::vector<int> arr)
 
424
  {
 
425
    unsigned int i;
 
426
    double angle;
 
427
    
 
428
    if (arr.size() != (_vrotor.size() + 1))
 
429
      return; // wrong size key
 
430
    
 
431
    double *rot = new double [_vrotor.size()+1];
 
432
    rot[0] = arr[0];
 
433
    
 
434
    double *c = mol.GetCoordinates();
 
435
    for (i = 0;i < _vrotor.size();++i)
 
436
      {
 
437
                                if (arr[i+1] == -1) // skip this rotor
 
438
                                        continue;
 
439
                                else {
 
440
                angle = _vres[i][arr[i+1]];
 
441
                while (angle < 0.0)
 
442
                angle += 360.0;
 
443
                while (angle > 360.0)
 
444
                angle -= 360.0;
 
445
                SetRotorToAngle(c,_vrotor[i].first,angle,_vrotor[i].second);
 
446
                                } // set an angle
 
447
      } // for rotors
 
448
  }
312
449
 
313
 
//Copies the coordinates in bc, NOT the pointers, into the object
314
 
void OBRotamerList::SetBaseCoordinateSets(vector<double*> bc, unsigned int N)
315
 
{
 
450
  //Copies the coordinates in bc, NOT the pointers, into the object
 
451
  void OBRotamerList::SetBaseCoordinateSets(vector<double*> bc, unsigned int N)
 
452
  {
316
453
    unsigned int i,j;
317
454
 
318
455
    //Clear out old data
319
 
    for (i=0 ; i<_c.size() ; i++)
320
 
        delete [] _c[i];
 
456
    for (i=0 ; i<_c.size() ; ++i)
 
457
      delete [] _c[i];
321
458
    _c.clear();
322
459
 
323
460
    //Copy new data
324
461
    double *c = NULL;
325
462
    double *cc= NULL;
326
 
    for (i=0 ; i<bc.size() ; i++)
327
 
    {
 
463
    for (i=0 ; i<bc.size() ; ++i)
 
464
      {
328
465
        c = new double [3*N];
329
466
        cc = bc[i];
330
 
        for (j=0 ; j<3*N ; j++)
331
 
            c[j] = cc[j];
 
467
        for (j=0 ; j<3*N ; ++j)
 
468
          c[j] = cc[j];
332
469
        _c.push_back(c);
333
 
    }
 
470
      }
334
471
    _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; 
 
472
  }
 
473
 
 
474
  //! Rotate the coordinates of 'atoms'
 
475
  //! such that tor == ang.
 
476
  //! Atoms in 'tor' should be ordered such that the 3rd atom is 
 
477
  //! the pivot around which atoms rotate (ang is in degrees)
 
478
  //! \todo This code is identical to OBMol::SetTorsion() and should be combined
 
479
  void SetRotorToAngle(double *c, OBAtom **ref,double ang,vector<int> atoms)
 
480
  {
 
481
    double v1x,v1y,v1z,v2x,v2y,v2z,v3x,v3y,v3z;
 
482
    double c1x,c1y,c1z,c2x,c2y,c2z,c3x,c3y,c3z;
 
483
    double c1mag,c2mag,radang,costheta,m[9];
 
484
    double x,y,z,mag,rotang,sn,cs,t,tx,ty,tz;
 
485
 
 
486
    int tor[4];
 
487
    tor[0] = ref[0]->GetCIdx();
 
488
    tor[1] = ref[1]->GetCIdx();
 
489
    tor[2] = ref[2]->GetCIdx();
 
490
    tor[3] = ref[3]->GetCIdx();
 
491
 
 
492
    //
 
493
    //calculate the torsion angle
 
494
    //
 
495
    v1x = c[tor[0]]   - c[tor[1]];   v2x = c[tor[1]]   - c[tor[2]];
 
496
    v1y = c[tor[0]+1] - c[tor[1]+1]; v2y = c[tor[1]+1] - c[tor[2]+1];
 
497
    v1z = c[tor[0]+2] - c[tor[1]+2]; v2z = c[tor[1]+2] - c[tor[2]+2];
 
498
    v3x = c[tor[2]]   - c[tor[3]];
 
499
    v3y = c[tor[2]+1] - c[tor[3]+1];
 
500
    v3z = c[tor[2]+2] - c[tor[3]+2];
 
501
 
 
502
    c1x = v1y*v2z - v1z*v2y;   c2x = v2y*v3z - v2z*v3y;
 
503
    c1y = -v1x*v2z + v1z*v2x;  c2y = -v2x*v3z + v2z*v3x;
 
504
    c1z = v1x*v2y - v1y*v2x;   c2z = v2x*v3y - v2y*v3x;
 
505
    c3x = c1y*c2z - c1z*c2y;
 
506
    c3y = -c1x*c2z + c1z*c2x;
 
507
    c3z = c1x*c2y - c1y*c2x; 
370
508
  
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));
 
509
    c1mag = c1x*c1x + c1y*c1y + c1z*c1z;
 
510
    c2mag = c2x*c2x + c2y*c2y + c2z*c2z;
 
511
    if (c1mag*c2mag < 0.01) costheta = 1.0; //avoid div by zero error
 
512
    else costheta = (c1x*c2x + c1y*c2y + c1z*c2z)/(sqrt(c1mag*c2mag));
375
513
 
376
 
  if (costheta < -0.999999) costheta = -0.999999f;
377
 
  if (costheta >  0.999999) costheta =  0.999999f;
 
514
    if (costheta < -0.999999) costheta = -0.999999;
 
515
    if (costheta >  0.999999) costheta =  0.999999;
378
516
                              
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;
 
517
    if ((v2x*c3x + v2y*c3y + v2z*c3z) > 0.0) radang = -acos(costheta);
 
518
    else                                     radang = acos(costheta);
 
519
 
 
520
    //
 
521
    // now we have the torsion angle (radang) - set up the rot matrix
 
522
    //
 
523
 
 
524
    //find the difference between current and requested
 
525
    rotang = (DEG_TO_RAD*ang) - radang; 
 
526
 
 
527
    sn = sin(rotang); cs = cos(rotang);t = 1 - cs;
 
528
    //normalize the rotation vector
 
529
    mag = sqrt(v2x*v2x + v2y*v2y + v2z*v2z);
 
530
    if (mag < 0.1) mag = 0.1; // avoid divide by zero error
 
531
    x = v2x/mag; y = v2y/mag; z = v2z/mag;
393
532
  
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
 
{
 
533
    //set up the rotation matrix
 
534
    m[0]= t*x*x + cs;     m[1] = t*x*y + sn*z;  m[2] = t*x*z - sn*y;
 
535
    m[3] = t*x*y - sn*z;  m[4] = t*y*y + cs;    m[5] = t*y*z + sn*x;
 
536
    m[6] = t*x*z + sn*y;  m[7] = t*y*z - sn*x;  m[8] = t*z*z + cs;
 
537
 
 
538
    //
 
539
    //now the matrix is set - time to rotate the atoms
 
540
    //
 
541
    tx = c[tor[1]];ty = c[tor[1]+1];tz = c[tor[1]+2];
 
542
    vector<int>::iterator i;int j;
 
543
    for (i = atoms.begin();i != atoms.end();++i)
 
544
      {
 
545
        j = ((*i)-1)*3;
 
546
        c[j] -= tx;c[j+1] -= ty;c[j+2]-= tz;
 
547
        x = c[j]*m[0] + c[j+1]*m[1] + c[j+2]*m[2];
 
548
        y = c[j]*m[3] + c[j+1]*m[4] + c[j+2]*m[5];
 
549
        z = c[j]*m[6] + c[j+1]*m[7] + c[j+2]*m[8];
 
550
        c[j] = x; c[j+1] = y; c[j+2] = z;
 
551
        c[j] += tx;c[j+1] += ty;c[j+2] += tz;
 
552
      }
 
553
  }
 
554
 
 
555
  int PackCoordinate(double c[3],double max[3])
 
556
  {
418
557
    int tmp;
419
 
    float cf;
 
558
    double cf;
420
559
    cf = c[0];
421
560
    tmp  = ((int)(cf*max[0])) << 20;
422
561
    cf = c[1];
424
563
    cf = c[2];
425
564
    tmp |= ((int)(cf*max[2]));
426
565
    return(tmp);
427
 
}
 
566
  }
428
567
 
429
 
void UnpackCoordinate(double c[3],double max[3],int tmp)
430
 
{
431
 
    float cf;
432
 
    cf = (float)(tmp>>20);
 
568
  void UnpackCoordinate(double c[3],double max[3],int tmp)
 
569
  {
 
570
    double cf;
 
571
    cf = (double)(tmp>>20);
433
572
    c[0] = cf;
434
573
    c[0] *= max[0];
435
 
    cf = (float)((tmp&0xffc00)>>10);
 
574
    cf = (double)((tmp&0xffc00)>>10);
436
575
    c[1] = cf;
437
576
    c[1] *= max[1];
438
 
    cf = (float)(tmp&0x3ff);
 
577
    cf = (double)(tmp&0x3ff);
439
578
    c[2] = cf;
440
579
    c[2] *= max[2];
441
 
}
 
580
  }
442
581
 
443
582
} //namespace OpenBabel
444
583