forcefield.h

Go to the documentation of this file.

/**********************************************************************
forcefield.h - Handle OBForceField class.

Copyright (C) 2006-2007 by Tim Vandermeersch <[email protected]>

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_FORCEFIELD_H
#define OB_FORCEFIELD_H

#include <vector>
#include <string>
#include <map>

#include <list>
#include <set>
#include <openbabel/babelconfig.h>
#include <openbabel/base.h>
#include <openbabel/mol.h>
#include <openbabel/plugin.h>
#include <openbabel/grid.h>
#include <openbabel/griddata.h>
#include <float.h>

namespace OpenBabel
{
  // log levels
#define OBFF_LOGLVL_NONE        0   //!< no output
#define OBFF_LOGLVL_LOW         1   //!< SteepestDescent progress... (no output from Energy())
#define OBFF_LOGLVL_MEDIUM      2   //!< individual energy terms
#define OBFF_LOGLVL_HIGH        3   //!< individual calculations and parameters

  // terms
#define OBFF_ENERGY             (1 << 0)   //!< all terms
#define OBFF_EBOND              (1 << 1)   //!< bond term
#define OBFF_EANGLE             (1 << 2)   //!< angle term
#define OBFF_ESTRBND            (1 << 3)   //!< strbnd term
#define OBFF_ETORSION           (1 << 4)   //!< torsion term
#define OBFF_EOOP               (1 << 5)   //!< oop term
#define OBFF_EVDW               (1 << 6)   //!< vdw term
#define OBFF_EELECTROSTATIC     (1 << 7)   //!< electrostatic term

  // constraint types
#define OBFF_CONST_IGNORE       (1 << 0)   //!< ignore the atom while setting up calculations
#define OBFF_CONST_ATOM         (1 << 1)   //!< fix the atom position
#define OBFF_CONST_ATOM_X       (1 << 2)   //!< fix the x coordinate of the atom position
#define OBFF_CONST_ATOM_Y       (1 << 3)   //!< fix the y coordinate of the atom position
#define OBFF_CONST_ATOM_Z       (1 << 4)   //!< fix the z coordinate of the atom position
#define OBFF_CONST_DISTANCE     (1 << 5)   //!< constrain distance length
#define OBFF_CONST_ANGLE        (1 << 6)   //!< constrain angle
#define OBFF_CONST_TORSION      (1 << 7)   //!< constrain torsion
#define OBFF_CONST_CHIRAL       (1 << 8)   //!< constrain chiral volume

  // mode arguments for SteepestDescent, ConjugateGradients, ...
#define OBFF_NUMERICAL_GRADIENT         (1 << 0)  //!< use numerical gradients
#define OBFF_ANALYTICAL_GRADIENT        (1 << 1)  //!< use analytical gradients

#define KCAL_TO_KJ      4.1868

  // inline if statements for logging.
#define IF_OBFF_LOGLVL_LOW    if(_loglvl >= OBFF_LOGLVL_LOW)
#define IF_OBFF_LOGLVL_MEDIUM if(_loglvl >= OBFF_LOGLVL_MEDIUM)
#define IF_OBFF_LOGLVL_HIGH   if(_loglvl >= OBFF_LOGLVL_HIGH)

  struct LineSearchType
  {
    enum {
      Simple, Newton2Num
    };
  };
  /*
  struct ConstraintType
  {
    enum {
      Ignore, Atom, AtomX, AtomY, AtomZ, Distance, Angle, Torsion, Chiral
    };
  };
  */

  class OBFPRT OBFFParameter {
  public:
    int         a, b, c, d;
    std::string _a, _b, _c, _d;

    std::vector<int>    _ipar;
    std::vector<double> _dpar;

    OBFFParameter& operator=(const OBFFParameter &ai)
      {
        if (this != &ai) {
          a = ai.a;
          b = ai.b;
          c = ai.c;
          d = ai.d;
          _a = ai._a;
          _b = ai._b;
          _c = ai._c;
          _d = ai._d;
          _ipar = ai._ipar;
          _dpar = ai._dpar;
        }

        return *this;
      }

    void clear ()
    {
      a = b = c = d = 0;
      _ipar.clear();
      _dpar.clear();
    }
  }; // class OBFFParameter

  // specific class introductions in forcefieldYYYY.cpp (for YYYY calculations)

  class OBFPRT OBFFCalculation2
  {
  public:
    double energy;
    OBAtom *a, *b;
    int idx_a, idx_b;
    double *pos_a, *pos_b;
    double force_a[3], force_b[3];
    virtual ~OBFFCalculation2()
    {
    }
    virtual void SetupPointers()
    {
      if (!a || !b) return;
      pos_a = a->GetCoordinate();
      idx_a = a->GetIdx();
      pos_b = b->GetCoordinate();
      idx_b = b->GetIdx();
    }
  };

  class OBFPRT OBFFCalculation3: public OBFFCalculation2
  {
  public:
    OBAtom *c;
    int idx_c;
    double *pos_c;
    double force_c[3];
    virtual ~OBFFCalculation3()
    {
    }
    virtual void SetupPointers()
    {
      if (!a || !b || !c) return;
      pos_a = a->GetCoordinate();
      idx_a = a->GetIdx();
      pos_b = b->GetCoordinate();
      idx_b = b->GetIdx();
      pos_c = c->GetCoordinate();
      idx_c = c->GetIdx();
    }
  };

  class OBFPRT OBFFCalculation4: public OBFFCalculation3
  {
  public:
    OBAtom *d;
    int idx_d;
    double *pos_d;
    double force_d[3];
    virtual ~OBFFCalculation4()
    {
    }
    void SetupPointers()
    {
      if (!a || !b || !c || !d) return;
      pos_a = a->GetCoordinate();
      idx_a = a->GetIdx();
      pos_b = b->GetCoordinate();
      idx_b = b->GetIdx();
      pos_c = c->GetCoordinate();
      idx_c = c->GetIdx();
      pos_d = d->GetCoordinate();
      idx_d = d->GetIdx();
    }
  };

  class OBFPRT OBFFConstraint
  {
  public:
    double factor, constraint_value;
    double rab0, rbc0;
    int type, ia, ib, ic, id;
    OBAtom *a, *b, *c, *d;
    vector3 grada, gradb, gradc, gradd;

    OBFFConstraint()
      {
        a = b = c = d = NULL;
        ia = ib = ic = id = 0;
        constraint_value = 0.0;
        factor = 0.0;
      }
    ~OBFFConstraint()
      {
      }

    vector3 GetGradient(int a)
    {
      if (a == ia)
        return grada;
      else if (a == ib)
        return gradb;
      else if (a == ic)
        return gradc;
      else if (a == id)
        return gradd;
      else
        return  VZero;
    }
  };

  class OBFPRT OBFFConstraints
  {
  public:
    OBFFConstraints();
    ~OBFFConstraints()
      {
        _constraints.clear();
        _ignored.Clear();
        _fixed.Clear();
        _Xfixed.Clear();
        _Yfixed.Clear();
        _Zfixed.Clear();
      }
    void Clear();
    double GetConstraintEnergy();
    vector3 GetGradient(int a);
    OBFFConstraints& operator=(const OBFFConstraints &ai)
      {
        if (this != &ai) {
          _constraints = ai._constraints;
          _ignored = ai._ignored;
          _fixed = ai._fixed;
          _Xfixed = ai._Xfixed;
          _Yfixed = ai._Yfixed;
          _Zfixed = ai._Zfixed;
        }
        return *this;
      }

    void Setup(OBMol &mol);

    // Set Constraints                                                     //


    void SetFactor(double factor);
    void AddIgnore(int a);
    void AddAtomConstraint(int a);
    void AddAtomXConstraint(int a);
    void AddAtomYConstraint(int a);
    void AddAtomZConstraint(int a);
    void AddDistanceConstraint(int a, int b, double length);
    void AddAngleConstraint(int a, int b, int c, double angle);
    void AddTorsionConstraint(int a, int b, int c, int d, double torsion);
    void DeleteConstraint(int index);

    // Get Constraints                                                     //


    double GetFactor();
    int Size() const;
    int GetConstraintType(int index) const;
    double GetConstraintValue(int index) const;
    int GetConstraintAtomA(int index) const;
    int GetConstraintAtomB(int index) const;
    int GetConstraintAtomC(int index) const;
    int GetConstraintAtomD(int index) const;
    bool IsIgnored(int a);
    bool IsFixed(int a);
    bool IsXFixed(int a);
    bool IsYFixed(int a);
    bool IsZFixed(int a);
    OBBitVec GetIgnoredBitVec() { return _ignored; }
    OBBitVec GetFixedBitVec() { return _fixed; }

  private:
    std::vector<OBFFConstraint> _constraints;
    OBBitVec    _ignored;
    OBBitVec    _fixed;
    OBBitVec    _Xfixed;
    OBBitVec    _Yfixed;
    OBBitVec    _Zfixed;
    double _factor;
  };

  // Class OBForceField
  // class introduction in forcefield.cpp
  class OBFPRT OBForceField : public OBPlugin
  {
    MAKE_PLUGIN(OBForceField)

    protected:

    OBFFParameter* GetParameter(int a, int b, int c, int d, std::vector<OBFFParameter> &parameter);
    OBFFParameter* GetParameter(const char* a, const char* b, const char* c, const char* d,
        std::vector<OBFFParameter> &parameter);
    int GetParameterIdx(int a, int b, int c, int d, std::vector<OBFFParameter> &parameter);

    vector3 NumericalDerivative(OBAtom *a, int terms = OBFF_ENERGY);
    vector3 NumericalSecondDerivative(OBAtom *a, int terms = OBFF_ENERGY);

    /*
     *   NEW gradients functions
     */

    void SetGradient(double *grad, int idx)
    {
      const int coordIdx = (idx - 1) * 3;
      for (unsigned int i = 0; i < 3; ++i) {
        _gradientPtr[coordIdx + i] = grad[i];
      }
    }

    void AddGradient(double *grad, int idx)
    {
      const int coordIdx = (idx - 1) * 3;
      for (unsigned int i = 0; i < 3; ++i) {
        _gradientPtr[coordIdx + i] += grad[i];
      }
    }

    virtual vector3 GetGradient(OBAtom *a, int /*terms*/ = OBFF_ENERGY)
    {
      const int coordIdx = (a->GetIdx() - 1) * 3;
      return _gradientPtr + coordIdx;
    }

    double* GetGradientPtr()
    {
      return _gradientPtr;
    }

    virtual void ClearGradients()
    {
      // We cannot use memset because IEEE floating point representations
      // are not guaranteed by C/C++ standard, but this loop can be
      // unrolled or vectorized by compilers
      for (unsigned int i = 0; i < _ncoords; ++i)
        _gradientPtr[i] = 0.0;
      //      memset(_gradientPtr, '\0', sizeof(double)*_ncoords);
    }

    bool IsInSameRing(OBAtom* a, OBAtom* b);

    // general variables
    OBMol       _mol; 
    bool        _init; 
    std::string _parFile; 
    bool        _validSetup; 
    double      *_gradientPtr; 
    // logging variables
    std::ostream* _logos; 
    char        _logbuf[BUFF_SIZE+1]; 
    int         _loglvl; 
    int         _origLogLevel;
    // conformer genereation (rotor search) variables
    int         _current_conformer; 
    std::vector<double> _energies; 
    // minimization variables
    double      _econv, _e_n1; 
    int         _cstep, _nsteps; 
    double      *_grad1; 
    unsigned int _ncoords; 
    int         _linesearch; 
    // molecular dynamics variables
    double      _timestep; 
    double      _temp; 
    double      *_velocityPtr; 
    // contraint varibles
    static OBFFConstraints _constraints; 
    static int _fixAtom; 
    static int _ignoreAtom; 
    // cut-off variables
    bool        _cutoff; 
    double      _rvdw; 
    double      _rele; 
    OBBitVec    _vdwpairs; 
    OBBitVec    _elepairs; 
    int         _pairfreq; 
    // group variables
    std::vector<OBBitVec> _intraGroup; 
    std::vector<OBBitVec> _interGroup; 
    std::vector<std::pair<OBBitVec, OBBitVec> > _interGroups; 

  public:
    virtual OBForceField* MakeNewInstance()=0;

    virtual ~OBForceField()
    {
      if (_grad1 != NULL) {
        delete [] _grad1;
        _grad1 = NULL;
      }
      if (_gradientPtr != NULL) {
        delete [] _gradientPtr;
        _gradientPtr = NULL;
      }
    }

    const char* TypeID()
    {
      return "forcefields";
    }

    static OBForceField* FindForceField(const std::string& ID)
    {
      return FindType(ID.c_str());
    }
    static OBForceField* FindForceField(const char *ID)
    {
      return FindType(ID);
    }
    /*
     *
     */
    void SetParameterFile(const std::string &filename)
    {
      _parFile = filename;
      _init = false;
    }
    virtual std::string GetUnit() { return std::string("au"); }
    /* Does this force field have analytical gradients defined for all
     * calculation components (bonds, angles, non-bonded, etc.)
     * If this is true, code should default to using OBFF_ANALYTICAL_GRADIENT
     * for SteepestDescent() or ConjugateGradients().
     * \return True if all analytical gradients are implemented.
     */
    virtual bool HasAnalyticalGradients() { return false; }
    bool Setup(OBMol &mol);
    bool Setup(OBMol &mol, OBFFConstraints &constraints);
    // move to protected in future version
    virtual bool ParseParamFile() { return false; }
    // move to protected in future version
    virtual bool SetTypes() { return false; }
    // move to protected in future version
    virtual bool SetFormalCharges() { return false; }
    // move to protected in future version
    virtual bool SetPartialCharges() { return false; }
    // move to protected in future version
    virtual bool SetupCalculations() { return false; }
    // move to protected in future version
    virtual bool SetupPointers() { return false; }
    bool IsSetupNeeded(OBMol &mol);
    bool GetAtomTypes(OBMol &mol);
    bool GetPartialCharges(OBMol &mol);



    bool GetCoordinates(OBMol &mol);
    bool UpdateCoordinates(OBMol &mol) {return GetCoordinates(mol); }
    bool GetConformers(OBMol &mol);
    bool UpdateConformers(OBMol &mol) { return GetConformers(mol); }
    bool SetCoordinates(OBMol &mol);
    bool SetConformers(OBMol &mol);
    OBGridData *GetGrid(double step, double padding, const char *type, double pchg);

    // Interacting groups                                                  //



    void AddIntraGroup(OBBitVec &group);
    void AddInterGroup(OBBitVec &group);
    void AddInterGroups(OBBitVec &group1, OBBitVec &group2);
    void ClearGroups();
    bool HasGroups();

    // Cut-off                                                             //



    void EnableCutOff(bool enable)
    {
      _cutoff = enable;
    }
    bool IsCutOffEnabled()
    {
      return _cutoff;
    }
    void SetVDWCutOff(double r)
    {
      _rvdw = r;
    }
    double GetVDWCutOff()
    {
      return _rvdw;
    }
    void SetElectrostaticCutOff(double r)
    {
      _rele = r;
    }
    double GetElectrostaticCutOff()
    {
      return _rele;
    }
    void SetUpdateFrequency(int f)
    {
      _pairfreq = f;
    }
    int GetUpdateFrequency()
    {
      return _pairfreq;
    }
    void UpdatePairsSimple();

    //void UpdatePairsGroup(); TODO

    unsigned int GetNumPairs();
    void EnableAllPairs()
    {
      // TODO: OBBitVec doesn't seem to be allocating it's memory correctly
      //_vdwpairs.SetRangeOn(0, _numpairs-1);
      //_elepairs.SetRangeOn(0, _numpairs-1);
    }

    // Energy Evaluation                                                   //



    virtual double Energy(bool UNUSED(gradients) = true) { return 0.0f; }
    virtual double E_Bond(bool UNUSED(gradients) = true) { return 0.0f; }
    virtual double E_Angle(bool UNUSED(gradients) = true) { return 0.0f; }
    virtual double E_StrBnd(bool UNUSED(gradients) = true) { return 0.0f; }
    virtual double E_Torsion(bool UNUSED(gradients) = true) { return 0.0f; }
    virtual double E_OOP(bool UNUSED(gradients) = true) { return 0.0f; }
    virtual double E_VDW(bool UNUSED(gradients) = true) { return 0.0f; }
    virtual double E_Electrostatic(bool UNUSED(gradients) = true) { return 0.0f; }

    // Logging                                                             //



    void PrintTypes();
    void PrintFormalCharges();
    void PrintPartialCharges();
    void PrintVelocities();
    bool SetLogFile(std::ostream *pos);
    bool SetLogLevel(int level);
    int GetLogLevel() { return _loglvl; }
    void OBFFLog(std::string msg)
    {
      if (!_logos)
        return;

      *_logos << msg;
    }
    void OBFFLog(const char *msg)
    {
      if (!_logos)
        return;

      *_logos << msg;
    }

    // Structure Generation                                                //



    void DistanceGeometry();
    void SystematicRotorSearch(unsigned int geomSteps = 2500);
    int SystematicRotorSearchInitialize(unsigned int geomSteps = 2500);
    bool SystematicRotorSearchNextConformer(unsigned int geomSteps = 2500);
    void RandomRotorSearch(unsigned int conformers, unsigned int geomSteps = 2500);
    void RandomRotorSearchInitialize(unsigned int conformers, unsigned int geomSteps = 2500);
    bool RandomRotorSearchNextConformer(unsigned int geomSteps = 2500);
    void WeightedRotorSearch(unsigned int conformers, unsigned int geomSteps);

    // Energy Minimization                                                 //



    void SetLineSearchType(int type)
    {
      _linesearch = type;
    }
    int GetLineSearchType()
    {
      return _linesearch;
    }
    vector3 LineSearch(OBAtom *atom, vector3 &direction);
    double LineSearch(double *currentCoords, double *direction);
    double Newton2NumLineSearch(double *direction);
    void   LineSearchTakeStep(double *origCoords, double *direction, double step);
    void SteepestDescent(int steps, double econv = 1e-6f, int method = OBFF_ANALYTICAL_GRADIENT);
    void SteepestDescentInitialize(int steps = 1000, double econv = 1e-6f, int method = OBFF_ANALYTICAL_GRADIENT);
    bool SteepestDescentTakeNSteps(int n);
    void ConjugateGradients(int steps, double econv = 1e-6f, int method = OBFF_ANALYTICAL_GRADIENT);
    void ConjugateGradientsInitialize(int steps = 1000, double econv = 1e-6f, int method = OBFF_ANALYTICAL_GRADIENT);
    bool ConjugateGradientsTakeNSteps(int n);

    // Molecular Dynamics                                                  //



    void GenerateVelocities();
    void CorrectVelocities();
    void MolecularDynamicsTakeNSteps(int n, double T, double timestep = 0.001, int method = OBFF_ANALYTICAL_GRADIENT);

    // Constraints                                                         //



    OBFFConstraints& GetConstraints();
    void SetConstraints(OBFFConstraints& constraints);
    void SetFixAtom(int index);
    void UnsetFixAtom();
    void SetIgnoreAtom(int index);
    void UnsetIgnoreAtom();

    static bool IgnoreCalculation(int a, int b);
    static bool IgnoreCalculation(int a, int b, int c);
    static bool IgnoreCalculation(int a, int b, int c, int d);


    // Validation                                                          //



    bool DetectExplosion();
    vector3 ValidateLineSearch(OBAtom *atom, vector3 &direction);
    void ValidateSteepestDescent(int steps);
    void ValidateConjugateGradients(int steps);
    virtual bool Validate() { return false; }
    virtual bool ValidateGradients() { return false; }
    vector3 ValidateGradientError(vector3 &numgrad, vector3 &anagrad);

    // Vector Analysis                                                     //



    static double VectorBondDerivative(double *pos_a, double *pos_b,
                                       double *force_a, double *force_b);
    static double VectorDistanceDerivative(const double* const pos_i, const double* const pos_j,
                                           double *force_i, double *force_j);
    static double VectorLengthDerivative(vector3 &a, vector3 &b);

    static double VectorAngleDerivative(double *pos_a, double *pos_b, double *pos_c,
                                        double *force_a, double *force_b, double *force_c);
    static double VectorAngleDerivative(vector3 &a, vector3 &b, vector3 &c);
    static double VectorOOPDerivative(double *pos_a, double *pos_b, double *pos_c, double *pos_d,
                                      double *force_a, double *force_b, double *force_c, double *force_d);
    static double VectorOOPDerivative(vector3 &a, vector3 &b, vector3 &c, vector3 &d);
    static double VectorTorsionDerivative(double *pos_a, double *pos_b, double *pos_c, double *pos_d,
                                          double *force_a, double *force_b, double *force_c, double *force_d);
    static double VectorTorsionDerivative(vector3 &a, vector3 &b, vector3 &c, vector3 &d);

    static void VectorSubtract(double *i, double *j, double *result)
    {
      for (unsigned int c = 0; c < 3; ++c)
        result[c] = i[c] - j[c];
    }

    static void VectorSubtract(const double* const i, const double* const j, double *result)
    {
      for (unsigned int c = 0; c < 3; ++c)
        result[c] = i[c] - j[c];
    }

    static void VectorAdd(double *i, double *j, double *result)
    {
      for (unsigned int c = 0; c < 3; ++c)
        result[c] = i[c] + j[c];
    }

    static void VectorDivide(double *i, double n, double *result)
    {
      for (unsigned int c = 0; c < 3; ++c)
        result[c] = i[c] / n;
    }

    static void VectorMultiply(double *i, double n, double *result)
    {
      for (unsigned int c = 0; c < 3; ++c)
        result[c] = i[c] * n;
    }

    static void VectorMultiply(const double* const i, const double n, double *result)
    {
      for (unsigned int c = 0; c < 3; ++c)
        result[c] = i[c] * n;
    }

    static void VectorSelfMultiply(double *i, double n)
    {
      for (unsigned int c = 0; c < 3; ++c)
        i[c] *= n;
    }

    static void VectorNormalize(double *i)
    {
      double length = VectorLength(i);
      for (unsigned int c = 0; c < 3; ++c)
        i[c] /= length;
    }

    static void VectorCopy(double *from, double *to)
    {
      for (unsigned int c = 0; c < 3; ++c)
        to[c] = from[c];
    }

    static double VectorLength(double *i)
    {
      return sqrt( i[0]*i[0] + i[1]*i[1] + i[2]*i[2] );
    }

    static double VectorDistance(double *pos_i, double *pos_j)
    {
      double ij[3];
      VectorSubtract(pos_i, pos_j, ij);
      const double rij = VectorLength(ij);
      return rij;
    }

    static double VectorAngle(double *i, double *j, double *k);

    static double VectorTorsion(double *i, double *j, double *k, double *l);

    static double VectorOOP(double *i, double *j, double *k, double *l);

    static void VectorClear(double *i)
    {
      for (unsigned int c = 0; c < 3; ++c)
        i[c] = 0.0;
    }

    static double VectorDot(double *i, double *j)
    {
      double result = 0.0;
      // Written as a loop for vectorization
      // Loop will be unrolled by compiler otherwise
      for (unsigned int c = 0; c < 3; ++c)
        result += i[c]*j[c];
      return result;
    }

    static void VectorCross(double *i, double *j, double *result)
    {
      result[0] =   i[1]*j[2] - i[2]*j[1];
      result[1] = - i[0]*j[2] + i[2]*j[0];
      result[2] =   i[0]*j[1] - i[1]*j[0];
    }

    static void PrintVector(double *i)
    {
      std::cout << "<" << i[0] << ", " << i[1] << ", " << i[2] << ">" << std::endl;
    }

  }; // class OBForceField

}// namespace OpenBabel

#endif   // OB_FORCEFIELD_H