src/openbabel/rotamer.cpp

/**********************************************************************
rotamer.cpp - Handle rotamer list data.
 
Copyright (C) 1998, 1999, 2000-2002 OpenEye Scientific Software, Inc.
Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
 
This file is part of the Open Babel project.
For more information, see <http://openbabel.sourceforge.net/>
 
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation version 2 of the License.
 
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
***********************************************************************/

#include "rotamer.hpp"
#include "mol.hpp"
#include "obutil.hpp"
#include "rotor.hpp"

#define OB_TITLE_SIZE     254
#define OB_BINARY_SETWORD 32

using namespace std;

namespace OpenBabel
{

//test byte ordering
static int SINT = 0x00000001;
static unsigned char *STPTR = (unsigned char*)&SINT;
const bool SwabInt = (STPTR[0]!=0);

#if !HAVE_RINT
inline double rint(double x)
{
    return ( (x < 0.0) ? ceil(x-0.5) : floor(x+0.5));
}
#endif

void SetRotorToAngle(double *c,OBAtom **ref,double ang,vector<int> atoms);

int Swab(int i)
{
    unsigned char tmp[4],c;
    memcpy(tmp,(char*)&i,sizeof(int));
    c = tmp[0];
    tmp[0] = tmp[3];
    tmp[3] = c;
    c = tmp[1];
    tmp[1] = tmp[2];
    tmp[2] = c;
    memcpy((char*)&i,tmp,sizeof(int));
    return(i);
}

OBRotamerList::~OBRotamerList()
{
    vector<unsigned char*>::iterator i;
    for (i = _vrotamer.begin();i != _vrotamer.end();i++)
        delete [] *i;

    vector<pair<OBAtom**,vector<int> > >::iterator j;
    for (j = _vrotor.begin();j != _vrotor.end();j++)
        delete [] j->first;

    //Delete the interal base coordinate list
    unsigned int k;
    for (k=0 ; k<_c.size() ; k++)
        delete [] _c[k];
}

void OBRotamerList::GetReferenceArray(unsigned char *ref)
{
    int j;
    vector<pair<OBAtom**,vector<int> > >::iterator i;
    for (j=0,i = _vrotor.begin();i != _vrotor.end();i++)
    {
        ref[j++] = (unsigned char)(i->first[0])->GetIdx();
        ref[j++] = (unsigned char)(i->first[1])->GetIdx();
        ref[j++] = (unsigned char)(i->first[2])->GetIdx();
        ref[j++] = (unsigned char)(i->first[3])->GetIdx();
    }
}

void OBRotamerList::Setup(OBMol &mol,OBRotorList &rl)
{
    //clear the old stuff out if necessary
    _vres.clear();
    vector<unsigned char*>::iterator j;
    for (j = _vrotamer.begin();j != _vrotamer.end();j++)
        delete [] *j;
    _vrotamer.clear();

    vector<pair<OBAtom**,vector<int> > >::iterator k;
    for (k = _vrotor.begin();k != _vrotor.end();k++)
        delete [] k->first;
    _vrotor.clear();

    //create the new list
    OBRotor *rotor;
    vector<OBRotor*>::iterator i;
    vector<int> children;

    int ref[4];
    OBAtom **atomlist;
    for (rotor = rl.BeginRotor(i);rotor;rotor = rl.NextRotor(i))
    {
        atomlist = new OBAtom* [4];
        rotor->GetDihedralAtoms(ref);
        atomlist[0] = mol.GetAtom(ref[0]);
        atomlist[1] = mol.GetAtom(ref[1]);
        atomlist[2] = mol.GetAtom(ref[2]);
        atomlist[3] = mol.GetAtom(ref[3]);
        mol.FindChildren(children,ref[1],ref[2]);
        _vrotor.push_back(pair<OBAtom**,vector<int> > (atomlist,children));
        _vres.push_back(rotor->GetResolution());
    }

    vector<double>::iterator n;
    vector<vector<double> >::iterator m;
    for (m = _vres.begin();m != _vres.end();m++)
        for (n = m->begin();n != m->end();n++)
            *n *= RAD_TO_DEG;
}

void OBRotamerList::Setup(OBMol &mol,unsigned char *ref,int nrotors)
{
    //clear the old stuff out if necessary
    _vres.clear();
    vector<unsigned char*>::iterator j;
    for (j = _vrotamer.begin();j != _vrotamer.end();j++)
        delete [] *j;
    _vrotamer.clear();

    vector<pair<OBAtom**,vector<int> > >::iterator k;
    for (k = _vrotor.begin();k != _vrotor.end();k++)
        delete [] k->first;
    _vrotor.clear();

    //create the new list
    int i;
    vector<int> children;

    int refatoms[4];
    OBAtom **atomlist;
    for (i = 0; i < nrotors; i++)
    {
        atomlist = new OBAtom* [4];
        refatoms[0] = (int)ref[i*4  ];
        refatoms[1] = (int)ref[i*4+1];
        refatoms[2] = (int)ref[i*4+2];
        refatoms[3] = (int)ref[i*4+3];
        mol.FindChildren(children,refatoms[1],refatoms[2]);
        atomlist[0] = mol.GetAtom(refatoms[0]);
        atomlist[1] = mol.GetAtom(refatoms[1]);
        atomlist[2] = mol.GetAtom(refatoms[2]);
        atomlist[3] = mol.GetAtom(refatoms[3]);
        _vrotor.push_back(pair<OBAtom**,vector<int> > (atomlist,children));
    }

}

void OBRotamerList::AddRotamer(double *c)
{
    int idx,size;
    double angle,res=255.0f/360.0f;
    vector3 v1,v2,v3,v4;

    unsigned char *rot = new unsigned char [_vrotor.size()+1];
    rot[0] = (char) 0;

    vector<pair<OBAtom**,vector<int> > >::iterator i;
    for (size=1,i = _vrotor.begin();i != _vrotor.end();i++,size++)
    {
        idx = (i->first[0])->GetCIdx();
        v1.Set(c[idx],c[idx+1],c[idx+2]);
        idx = (i->first[1])->GetCIdx();
        v2.Set(c[idx],c[idx+1],c[idx+2]);
        idx = (i->first[2])->GetCIdx();
        v3.Set(c[idx],c[idx+1],c[idx+2]);
        idx = (i->first[3])->GetCIdx();
        v4.Set(c[idx],c[idx+1],c[idx+2]);

        angle = CalcTorsionAngle(v1,v2,v3,v4);
        while (angle < 0.0f)
            angle += 360.0f;
        while (angle > 360.0f)
            angle -= 360.0f;
        rot[size] = (unsigned char)rint(angle*res);
    }

    _vrotamer.push_back(rot);
}

void OBRotamerList::AddRotamer(int *arr)
{
    unsigned int i;
    double angle,res=255.0f/360.0f;

    unsigned char *rot = new unsigned char [_vrotor.size()+1];
    rot[0] = (unsigned char)arr[0];

    for (i = 0;i < _vrotor.size();i++)
    {
        angle = _vres[i][arr[i+1]];
        while (angle < 0.0f)
            angle += 360.0f;
        while (angle > 360.0f)
            angle -= 360.0f;
        rot[i+1] = (unsigned char)rint(angle*res);
    }
    _vrotamer.push_back(rot);
}

void OBRotamerList::AddRotamer(unsigned char *arr)
{
    unsigned int i;
    double angle,res=255.0f/360.0f;

    unsigned char *rot = new unsigned char [_vrotor.size()+1];
    rot[0] = (unsigned char)arr[0];

    for (i = 0;i < _vrotor.size();i++)
    {
        angle = _vres[i][(int)arr[i+1]];
        while (angle < 0.0f)
            angle += 360.0f;
        while (angle > 360.0f)
            angle -= 360.0f;
        rot[i+1] = (unsigned char)rint(angle*res);
    }
    _vrotamer.push_back(rot);
}

void OBRotamerList::AddRotamers(unsigned char *arr,int nrotamers)
{
    unsigned int size;
    int i;

    size = (unsigned int)_vrotor.size()+1;
    for (i = 0;i < nrotamers;i++)
    {
        unsigned char *rot = new unsigned char [size];
        memcpy(rot,&arr[i*size],sizeof(char)*size);
        _vrotamer.push_back(rot);
    }
}

void OBRotamerList::ExpandConformerList(OBMol &mol,vector<double*> &clist)
{
    unsigned int j;
    double angle,invres=360.0f/255.0f;
    unsigned char *conf;
    vector<double*> tmpclist;
    vector<unsigned char*>::iterator i;

    for (i = _vrotamer.begin();i != _vrotamer.end();i++)
    {
        conf = *i;
        double *c = new double [mol.NumAtoms()*3];
        memcpy(c,clist[(int)conf[0]],sizeof(double)*mol.NumAtoms()*3);

        for (j = 0;j < _vrotor.size();j++)
        {
            angle = invres*((double)conf[j+1]);
            if (angle > 180.0)
                angle -= 360.0;
            SetRotorToAngle(c,_vrotor[j].first,angle,_vrotor[j].second);
        }
        tmpclist.push_back(c);
    }

    //transfer the conf list
    vector<double*>::iterator k;
    for (k = clist.begin();k != clist.end();k++)
        delete [] *k;
    clist = tmpclist;
}

//! Create a conformer list using the internal base set of coordinates
vector<double*> OBRotamerList::CreateConformerList(OBMol& mol)
{
    unsigned int j;
    double angle,invres=360.0f/255.0f;
    unsigned char *conf;
    vector<double*> tmpclist;
    vector<unsigned char*>::iterator i;

    for (i = _vrotamer.begin();i != _vrotamer.end();i++)
    {
        conf = *i;
        double *c = new double [mol.NumAtoms()*3];
        memcpy(c,_c[(int)conf[0]],sizeof(double)*mol.NumAtoms()*3);

        for (j = 0;j < _vrotor.size();j++)
        {
            angle = invres*((double)conf[j+1]);
            if (angle > 180.0)
                angle -= 360.0;
            SetRotorToAngle(c,_vrotor[j].first,angle,_vrotor[j].second);
        }
        tmpclist.push_back(c);
    }

    return tmpclist;
}

//Copies the coordinates in bc, NOT the pointers, into the object
void OBRotamerList::SetBaseCoordinateSets(vector<double*> bc, unsigned int N)
{
    unsigned int i,j;

    //Clear out old data
    for (i=0 ; i<_c.size() ; i++)
        delete [] _c[i];
    _c.clear();

    //Copy new data
    double *c = NULL;
    double *cc= NULL;
    for (i=0 ; i<bc.size() ; i++)
    {
        c = new double [3*N];
        cc = bc[i];
        for (j=0 ; j<3*N ; j++)
            c[j] = cc[j];
        _c.push_back(c);
    }
    _NBaseCoords = N;
}

//! Rotate the coordinates of 'atoms'
//! such that tor == ang - atoms in 'tor' should be ordered such 
//! that the 3rd atom is the pivot around which atoms rotate
//! ang is in degrees
void SetRotorToAngle(double *c, OBAtom **ref,double ang,vector<int> atoms)
{
  double v1x,v1y,v1z,v2x,v2y,v2z,v3x,v3y,v3z;
  double c1x,c1y,c1z,c2x,c2y,c2z,c3x,c3y,c3z;
  double c1mag,c2mag,radang,costheta,m[9];
  double x,y,z,mag,rotang,sn,cs,t,tx,ty,tz;

  int tor[4];
  tor[0] = ref[0]->GetCIdx();
  tor[1] = ref[1]->GetCIdx();
  tor[2] = ref[2]->GetCIdx();
  tor[3] = ref[3]->GetCIdx();

  //
  //calculate the torsion angle
  //
  v1x = c[tor[0]]   - c[tor[1]];   v2x = c[tor[1]]   - c[tor[2]];
  v1y = c[tor[0]+1] - c[tor[1]+1]; v2y = c[tor[1]+1] - c[tor[2]+1];
  v1z = c[tor[0]+2] - c[tor[1]+2]; v2z = c[tor[1]+2] - c[tor[2]+2];
  v3x = c[tor[2]]   - c[tor[3]];
  v3y = c[tor[2]+1] - c[tor[3]+1];
  v3z = c[tor[2]+2] - c[tor[3]+2];

  c1x = v1y*v2z - v1z*v2y;   c2x = v2y*v3z - v2z*v3y;
  c1y = -v1x*v2z + v1z*v2x;  c2y = -v2x*v3z + v2z*v3x;
  c1z = v1x*v2y - v1y*v2x;   c2z = v2x*v3y - v2y*v3x;
  c3x = c1y*c2z - c1z*c2y;
  c3y = -c1x*c2z + c1z*c2x;
  c3z = c1x*c2y - c1y*c2x; 
  
  c1mag = c1x*c1x + c1y*c1y + c1z*c1z;
  c2mag = c2x*c2x + c2y*c2y + c2z*c2z;
  if (c1mag*c2mag < 0.01) costheta = 1.0; //avoid div by zero error
  else costheta = (c1x*c2x + c1y*c2y + c1z*c2z)/(sqrt(c1mag*c2mag));

  if (costheta < -0.999999) costheta = -0.999999f;
  if (costheta >  0.999999) costheta =  0.999999f;
                              
  if ((v2x*c3x + v2y*c3y + v2z*c3z) > 0.0) radang = -acos(costheta);
  else                                     radang = acos(costheta);

  //
  // now we have the torsion angle (radang) - set up the rot matrix
  //

  //find the difference between current and requested
  rotang = (DEG_TO_RAD*ang) - radang; 

  sn = sin(rotang); cs = cos(rotang);t = 1 - cs;
  //normalize the rotation vector
  mag = sqrt(v2x*v2x + v2y*v2y + v2z*v2z);
  x = v2x/mag; y = v2y/mag; z = v2z/mag;
  
  //set up the rotation matrix
  m[0]= t*x*x + cs;     m[1] = t*x*y + sn*z;  m[2] = t*x*z - sn*y;
  m[3] = t*x*y - sn*z;  m[4] = t*y*y + cs;    m[5] = t*y*z + sn*x;
  m[6] = t*x*z + sn*y;  m[7] = t*y*z - sn*x;  m[8] = t*z*z + cs;

  //
  //now the matrix is set - time to rotate the atoms
  //
  tx = c[tor[1]];ty = c[tor[1]+1];tz = c[tor[1]+2];
  vector<int>::iterator i;int j;
  for (i = atoms.begin();i != atoms.end();i++)
    {
      j = ((*i)-1)*3;
      c[j] -= tx;c[j+1] -= ty;c[j+2]-= tz;
      x = c[j]*m[0] + c[j+1]*m[1] + c[j+2]*m[2];
      y = c[j]*m[3] + c[j+1]*m[4] + c[j+2]*m[5];
      z = c[j]*m[6] + c[j+1]*m[7] + c[j+2]*m[8];
      c[j] = x; c[j+1] = y; c[j+2] = z;
      c[j] += tx;c[j+1] += ty;c[j+2] += tz;
    }
}

int PackCoordinate(double c[3],double max[3])
{
    int tmp;
    float cf;
    cf = c[0];
    tmp  = ((int)(cf*max[0])) << 20;
    cf = c[1];
    tmp |= ((int)(cf*max[1])) << 10;
    cf = c[2];
    tmp |= ((int)(cf*max[2]));
    return(tmp);
}

void UnpackCoordinate(double c[3],double max[3],int tmp)
{
    float cf;
    cf = (float)(tmp>>20);
    c[0] = cf;
    c[0] *= max[0];
    cf = (float)((tmp&0xffc00)>>10);
    c[1] = cf;
    c[1] *= max[1];
    cf = (float)(tmp&0x3ff);
    c[2] = cf;
    c[2] *= max[2];
}

} //namespace OpenBabel

//! \file rotamer.cpp
//! \brief Handle rotamer list data.
Revision:	741
Committed:	Wed Nov 16 19:42:11 2005 UTC (19 years, 5 months ago) by tim
File size:	12652 byte(s)
Log Message:	adding openbabel
#	User	Rev	Content
1	tim	741	/**********************************************************************
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[isize],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 = v1yv2z - v1zv2y; c2x = v2yv3z - v2zv3y;
365			c1y = -v1xv2z + v1zv2x; c2y = -v2xv3z + v2zv3x;
366			c1z = v1xv2y - v1yv2x; c2z = v2xv3y - v2yv3x;
367			c3x = c1yc2z - c1zc2y;
368			c3y = -c1xc2z + c1zc2x;
369			c3z = c1xc2y - c1yc2x;
370
371			c1mag = c1xc1x + c1yc1y + c1z*c1z;
372			c2mag = c2xc2x + c2yc2y + c2z*c2z;
373			if (c1mag*c2mag < 0.01) costheta = 1.0; //avoid div by zero error
374			else costheta = (c1xc2x + c1yc2y + c1zc2z)/(sqrt(c1magc2mag));
375
376			if (costheta < -0.999999) costheta = -0.999999f;
377			if (costheta > 0.999999) costheta = 0.999999f;
378
379			if ((v2xc3x + v2yc3y + 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(v2xv2x + v2yv2y + v2z*v2z);
392			x = v2x/mag; y = v2y/mag; z = v2z/mag;
393
394			//set up the rotation matrix
395			m[0]= txx + cs; m[1] = txy + snz; m[2] = txz - sny;
396			m[3] = txy - snz; m[4] = tyy + cs; m[5] = tyz + snx;
397			m[6] = txz + sny; m[7] = tyz - snx; m[8] = tzz + 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.