rotor.h

Go to the documentation of this file.

/**********************************************************************
rotor.h - Rotate torsional according to rotor rules.

Copyright (C) 1998-2000 by 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.org/>

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

#ifndef OB_ROTOR_H
#define OB_ROTOR_H

#include <openbabel/parsmart.h>
#include <openbabel/typer.h>
#include <openbabel/bitvec.h>

#ifdef UNUSED
#elif (__GNUC__ == 4)
# define UNUSED(x) UNUSED_ ## x __attribute__((unused))
#elif defined(__LCLINT__)
# define UNUSED(x) /*@unused@*/ x
#else
# define UNUSED(x) x
#endif

namespace OpenBabel
{
  class OBRing;

#ifndef SQUARE
#define SQUARE(x) ((x)*(x))
#endif

  class OBAPI OBRotorRule
  {
    int                 _ref[4]; 
    double              _delta;  
    std::string         _s;      
    OBSmartsPattern*    _sp;     
    std::vector<double> _vals;   
  public:

  OBRotorRule(char *buffer,int ref[4],std::vector<double> &vals,double d):
    _delta(d), _s(buffer), _vals(vals)
    {
      _sp = new OBSmartsPattern;
      _sp->Init(buffer);
      memcpy(_ref,ref,sizeof(int)*4);
    }

    ~OBRotorRule()
      {
        if (_sp)
          {
            delete _sp;
            _sp = NULL;
          }
      }

    bool    IsValid()    {        return(_sp->IsValid());       }
    void    GetReferenceAtoms(int ref[4]) { memcpy(ref,_ref,sizeof(int)*4); }
    void    SetDelta(double d)    {       _delta = d;           }
    double  GetDelta()            {       return(_delta);       }
    std::vector<double>   &GetTorsionVals()    { return(_vals); }
    std::string  &GetSmartsString(){      return(_s);           }
    OBSmartsPattern *GetSmartsPattern() {  return(_sp);         }
  };

  class OBAPI OBRotorRules : public OBGlobalDataBase
  {
    bool                       _quiet;  
    std::vector<OBRotorRule*>  _vr;     
    std::vector<double>        _sp3sp3; 
    std::vector<double>        _sp3sp2; 
    std::vector<double>        _sp2sp2; 
  public:
    OBRotorRules();
    ~OBRotorRules();

    void ParseLine(const char*);
    size_t GetSize()                 { return _vr.size();}

    void SetFilename(std::string &s)       { _filename = s;    }

    void GetRotorIncrements(OBMol& mol,OBBond* bond,int refs[4],
                            std::vector<double> &vals,double &delta);
    void Quiet()                           { _quiet=true;      }
  };

  class OBAPI OBRotor
  {
    int _idx; 
    std::vector<int> _rotatoms; 
    double _imag, _refang; 
    OBBond *_bond; 
    std::vector<int> _ref, _torsion; 
    OBBitVec _fixedatoms,_fixedbonds, _evalatoms; 
    std::vector<double> _torsionAngles;  
    std::vector<double> _invmag; 
    std::vector<std::vector<double> > _sn,_cs,_t; 
    std::vector<OBRing *> _rings; 
  public:
    OBRotor();
    ~OBRotor()
      {
      }


    void SetBond(OBBond *bond)
    {
      _bond = bond;
      SetRings();
    }
    void SetRings();
    void SetIdx(int idx)
    {
      _idx = idx;
    }
    void SetDihedralAtoms(std::vector<int> &ref);
    void SetDihedralAtoms(int ref[4]);
    void SetRotAtoms(std::vector<int> &atoms);
    void SetTorsionValues(std::vector<double> &angles)
    {
      _torsionAngles = angles;
    }
    void SetFixedBonds(OBBitVec &bv)
    {
      _fixedbonds = bv;
    }


    inline void SetToAngle(double *coordinates, double setang)
    {
      double /*dx,dy,dz,*/ sn,cs,t,ang,mag;
      // compute the angle to rotate (radians)
      ang = setang - CalcTorsion(coordinates);
      // if the angle to rotate is too small, we're done
      if (fabs(ang) < 1e-5)
        return;

      // compute the bond length
      mag = CalcBondLength(coordinates);
      // compute some rotation matrix elements
      sn = sin(ang);
      cs = cos(ang);
      t = 1 - cs;

      // perform rotation
      Set(coordinates, sn, cs, t, 1.0 / mag);
    }
    void SetRotor(double *coordinates, int next, int prev = -1);
    void Set(double *coordinates, double sine, double cosine, double translation, double invmag);
    void Precompute(double *coordinates);
    void Set(double *coordinates, int idx);
    void Precalc(std::vector<double*> &conformers);
    void Set(double *coordinates, int conformer, int idx)
    {
      Set(coordinates, _sn[conformer][idx], _cs[conformer][idx], _t[conformer][idx], _invmag[conformer]);
    }



    OBBond *GetBond()
    {
      return(_bond);
    }
    size_t Size()
    {
      return _torsionAngles.size();
    }
    int GetIdx() const
    {
      return _idx;
    }
    void GetDihedralAtoms(int ref[4])
    {
      for (int i = 0; i < 4; ++i)
        ref[i] = _ref[i];
    }
    std::vector<int> &GetDihedralAtoms()
      {
        return _ref;
      }
    const std::vector<int>& GetRotAtoms() const
    {
      return _rotatoms;
    }
    const std::vector<double> &GetTorsionValues() const
    {
      return _torsionAngles;
    }
    OBBitVec &GetFixedBonds()
      {
        return _fixedbonds;
      }
    double CalcTorsion(double *coordinates);
    double CalcBondLength(double *coordinates);


    std::vector<double>::iterator BeginTorIncrement()
      {
        return _torsionAngles.begin();
      }
    std::vector<double>::iterator EndTorIncrement()
      {
        return _torsionAngles.end();
      }

    void RemoveSymTorsionValues(int);


    void SetDelta(double UNUSED(d)) {}
    double GetDelta() { return 10.0; }
    OBBitVec &GetFixedAtoms() { return _fixedatoms; }
    void SetFixedAtoms(OBBitVec &bv) { _fixedatoms = bv; }
    OBBitVec &GetEvalAtoms() { return _evalatoms; }
    void SetEvalAtoms(OBBitVec &bv) { _evalatoms = bv; }
    void* GetRotAtoms() { return &_rotatoms; }
    std::vector<double> &GetResolution() { return _torsionAngles; }
    void SetNumCoords(int UNUSED(nc)) {}

  };


  typedef std::vector<OBRotor*>::iterator OBRotorIterator;

  class OBAPI OBRotorList
  {
    bool _quiet;                    
    bool _removesym;                
    bool _ringRotors;               
    OBBitVec _fixedatoms, _fixedbonds; 
    OBRotorRules _rr;               
    std::vector<int> _dffv;         
    std::vector<OBRotor*> _rotor;   
    std::vector<std::pair<OBSmartsPattern*,std::pair<int,int> > > _vsym2;
    std::vector<std::pair<OBSmartsPattern*,std::pair<int,int> > > _vsym3;
  public:
    OBRotorList();
    ~OBRotorList();

    void Clear();
    size_t Size()
    {
      return _rotor.size();
    }
    bool IsFixedBond(OBBond*);
    bool HasFixedBonds()
    {
      return !_fixedbonds.IsEmpty();
    }
    void RemoveSymVals(OBMol&);

    bool HasRingRotors()
    { return _ringRotors; }


    bool Setup(OBMol &mol, bool sampleRings = false);
    void SetFixedBonds(OBBitVec &fix)
    {
      _fixedbonds = fix;
      _fixedatoms.Clear();
    }
    void Init(std::string &fname)
    {
      _rr.SetFilename(fname);
      _rr.Init();
    }
    void SetQuiet() {
      _quiet=true;
      _rr.Quiet();
    }
    bool SetRotAtoms(OBMol&);
    bool FindRotors(OBMol &mol, bool sampleRingBonds = false);
    bool SetEvalAtoms(OBMol&);
    bool AssignTorVals(OBMol &);



    OBRotor *BeginRotor(OBRotorIterator &i)
    {
      i = _rotor.begin();
      return((i ==_rotor.end()) ? NULL:*i);
    }
    OBRotor *NextRotor(OBRotorIterator &i)
    {
      ++i;
      return((i ==_rotor.end()) ? NULL:*i);
    }
    OBRotorIterator BeginRotors()   { return(_rotor.begin()); }
    OBRotorIterator EndRotors()     { return(_rotor.end());   }

    // Not declared
    bool   IdentifyEvalAtoms(OBMol &mol) { return SetEvalAtoms(mol); }
    void SetFixAtoms(OBBitVec &fix)
    {
      _fixedatoms = fix;
      _fixedbonds.Clear();
    }
    bool HasFixedAtoms()
    {
      return(!_fixedatoms.IsEmpty());
    }
    void IgnoreSymmetryRemoval()    { _removesym = false;}
    void SetRotAtomsByFix(OBMol&);

  };

  class rotor_digit {
  public:
    rotor_digit(unsigned int rs)
      {
        resolution_size = rs;
        state = 0;
      }

    rotor_digit()
      {
        resolution_size = 0;
        state = 0;
      }

    void set_size(unsigned int rs)
    {
      resolution_size = rs;
      state = 0;
    }

    void set_state(int st)
    {
      state = st;
    }

    int get_state()
    {
      return state;
    }

    unsigned int size()
    {
      return resolution_size;
    }

    bool next()
    {
      if (state < static_cast<int>(resolution_size - 1)) {
        ++state;
        return false;
      } else
        state = 0;

      return true;
    }
  private:
    unsigned int resolution_size;
    int state;
#ifndef SWIG
  } typedef rotor_digit;
#else
};
#endif

  class OBAPI OBRotorKeys
  {
    public:
      OBRotorKeys()
      {
        _vr.clear();
      }

      void Clear(){
        _vr.clear();
      }

      unsigned int NumKeys()
      {
        unsigned int numKeys = 0;

        while (Next())
          numKeys++;

        return numKeys;
      }

      void AddRotor(unsigned int size)
      {
        rotor_digit rd(size);
        _vr.push_back(rd);
      }

      bool Next()
      {
        if(_vr.size() == 0)
          return false;

        bool carry = _vr[0].next();
        unsigned int i = 1;
        while (carry) {
          if(i == _vr.size())
            return false;

          carry = _vr[i].next();
          i++;
        }
        return true;
      }

      std::vector<int> GetKey()
      {
        std::vector<int> rt;
        rt.clear();
        rt.push_back(0);
        for(unsigned int i = 0; i < _vr.size(); i++){
          rt.push_back(_vr[i].get_state());
        }

        return rt;
      }

    private:
      std::vector<rotor_digit> _vr;
  };


} // end namespace OpenBabel

#endif // OB_ROTOR_H