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 <openbabel/babelconfig.h>
#include <openbabel/mol.h> // TODO: Move OBMol code out of the header (use OBMol*)
#include <openbabel/atom.h> // TODO: Move OBAtom code out of the header
#include <openbabel/plugin.h>
#include <openbabel/bitvec.h>
#include <float.h>

namespace OpenBabel
{
 class OBGridData;

 // log levels
#define OBFF_LOGLVL_NONE 0 
#define OBFF_LOGLVL_LOW 1 
#define OBFF_LOGLVL_MEDIUM 2 
#define OBFF_LOGLVL_HIGH 3 

 // terms
#define OBFF_ENERGY (1 << 0) 
#define OBFF_EBOND (1 << 1) 
#define OBFF_EANGLE (1 << 2) 
#define OBFF_ESTRBND (1 << 3) 
#define OBFF_ETORSION (1 << 4) 
#define OBFF_EOOP (1 << 5) 
#define OBFF_EVDW (1 << 6) 
#define OBFF_EELECTROSTATIC (1 << 7) 

 // constraint types
#define OBFF_CONST_IGNORE (1 << 0) 
#define OBFF_CONST_ATOM (1 << 1) 
#define OBFF_CONST_ATOM_X (1 << 2) 
#define OBFF_CONST_ATOM_Y (1 << 3) 
#define OBFF_CONST_ATOM_Z (1 << 4) 
#define OBFF_CONST_DISTANCE (1 << 5) 
#define OBFF_CONST_ANGLE (1 << 6) 
#define OBFF_CONST_TORSION (1 << 7) 
#define OBFF_CONST_CHIRAL (1 << 8) 

 // mode arguments for SteepestDescent, ConjugateGradients, ...
#define OBFF_NUMERICAL_GRADIENT (1 << 0) 
#define OBFF_ANALYTICAL_GRADIENT (1 << 1) 

const double KCAL_TO_KJ = 4.1868;
const double GAS_CONSTANT = 8.31446261815324e-3 / KCAL_TO_KJ; 

 // 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 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, _gconv, _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 unsigned int _fixAtom; 
 static unsigned int _ignoreAtom; 
 // cut-off variables
 bool _cutoff; 
 double _rvdw; 
 double _rele; 
 double _epsilon; 
 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 SetDielectricConstant(double epsilon)
 {
 _epsilon = epsilon;
 }
 /* Get the dielectric permittivity used for electrostatic calculations
 * \rreturn The current relative permittivity
 */
 double GetDielectricConstant()
 {
 return _epsilon;
 }
 void SetUpdateFrequency(int f)
 {
 _pairfreq = f;
 }
 int GetUpdateFrequency()
 {
 return _pairfreq;
 }
 void UpdatePairsSimple();

 //void UpdatePairsGroup(); TODO

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

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

 double* GetGradientPtr()
 {
 return _gradientPtr;
 }

 // 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, bool sampleRingBonds = false);
 int SystematicRotorSearchInitialize(unsigned int geomSteps = 2500, bool sampleRingBonds = false);
 bool SystematicRotorSearchNextConformer(unsigned int geomSteps = 2500);
 void RandomRotorSearch(unsigned int conformers, unsigned int geomSteps = 2500,
 bool sampleRingBonds = false);
 void RandomRotorSearchInitialize(unsigned int conformers, unsigned int geomSteps = 2500,
 bool sampleRingBonds = false);
 bool RandomRotorSearchNextConformer(unsigned int geomSteps = 2500);
 void WeightedRotorSearch(unsigned int conformers, unsigned int geomSteps,
 bool sampleRingBonds = false);
 int FastRotorSearch(bool permute = true);

#ifdef HAVE_EIGEN
 int DiverseConfGen(double rmsd, unsigned int nconfs = 0, double energy_gap = 50, bool verbose = false);
#endif

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