/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. |