/**********************************************************************

forcefieldmmff94.h - MMFF94

Copyright (C) 2006 by Tim Vandermeersch <tim.vandermeersch@gmail.com>

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.

***********************************************************************/

#include <vector>

#include <string>

#include <map>

#include <openbabel/parsmart.h>

#include <openbabel/forcefield.h>

#include <openbabel/base.h>

#include <openbabel/mol.h>

namespace OpenBabel

{

  class OBFFBondCalculationMMFF94 : public OBFFCalculation2

  {

    public:

      int bt; // bondtype (BTIJ)

      double kb, r0, rab, delta;

      template<bool> void Compute();

  };

  class OBFFAngleCalculationMMFF94 : public OBFFCalculation3

  {

    public:

      int at; //angletype (ATIJK)

      bool linear;

      double ka, theta, theta0, delta;

      template<bool> void Compute();

  };

  class OBFFStrBndCalculationMMFF94 : public OBFFCalculation3

  {

    public:

      int sbt; //strbndtype (SBTIJK)

      double kbaABC, kbaCBA, theta0, rab0, rbc0, delta_theta, delta_rab, delta_rbc;

      double theta, rab, rbc;

      double force_ab_a[3], force_ab_b[3], force_bc_b[3], force_bc_c[3];

      double force_abc_a[3], force_abc_b[3], force_abc_c[3];

      template<bool> void Compute();

  };

  class OBFFTorsionCalculationMMFF94 : public OBFFCalculation4

  {

    public:

      int tt; //torsiontype (TTIJKL)

      double v1, v2, v3, tor, cosine;

      template<bool> void Compute();

  };

  class OBFFOOPCalculationMMFF94 : public OBFFCalculation4

  {

    public:

      double koop, angle;

      template<bool> void Compute();

  };

  class OBFFVDWCalculationMMFF94 : public OBFFCalculation2

  {

    public:

      int aDA, bDA; // hydrogen donor/acceptor (A=1, D=2, neither=0)

      double rab, epsilon, alpha_a, alpha_b, Na, Nb, Aa, Ab, Ga, Gb;

      double R_AB, R_AB7/*, erep, erep7, eattr*/;

      int pairIndex; // index into iteration using FOR_PAIRS_OF_MOL(..., _mol)

      template<bool> void Compute();

  };

  class OBFFElectrostaticCalculationMMFF94 : public OBFFCalculation2

  {

    public:

      double qq, rab;

      int pairIndex; // index into iteration using FOR_PAIRS_OF_MOL(..., _mol)

      template<bool> void Compute();

  };

  // Class OBForceFieldMMFF94

  // class introduction in forcefieldmmff94.cpp

  class OBForceFieldMMFF94: public OBForceField

  {

    protected:

      //! \return Parses the parameter file

      bool ParseParamFile() override;

      bool ParseParamProp(std::string &filename);

      bool ParseParamDef(std::string &filename);

      bool ParseParamBond(std::string &filename);

      bool ParseParamBndk(std::string &filename);

      bool ParseParamAngle(std::string &filename);

      bool ParseParamStrBnd(std::string &filename);

      bool ParseParamDfsb(std::string &filename);

      bool ParseParamOOP(std::string &filename);

      bool ParseParamTorsion(std::string &filename);

      bool ParseParamVDW(std::string &filename);

      bool ParseParamCharge(std::string &filename);

      bool ParseParamPbci(std::string &filename);

      //! detect which rings are aromatic

      bool PerceiveAromatic();

      //! \return Get the MMFF94 atom type for atom

      int GetType(OBAtom *atom);

      //! \return Sets atomtypes to MMFF94 in _mol

      bool SetTypes() override;

      //! fill OBFFXXXCalculation vectors

      bool SetupCalculations() override;

      //! Setup pointers in OBFFXXXCalculation vectors

      bool SetupPointers() override;

      //!  Sets formal charges

      bool SetFormalCharges() override;

      //!  Sets partial charges

      bool SetPartialCharges() override;

      //! \return The row of the element atom in the periodic table

      int GetElementRow(OBAtom *atom);

      //! \return The bond type (BTIJ)

      int GetBondType(OBAtom* a, OBAtom* b);

      //! \return The angle type (ATIJK)

      int GetAngleType(OBAtom* a, OBAtom* b, OBAtom *c);

      //! \return The stretch-bend type (SBTIJK)

      int GetStrBndType(OBAtom* a, OBAtom* b, OBAtom *c);

      //! \return The torsion type (TTIJKL)

      int GetTorsionType(OBAtom* a, OBAtom* b, OBAtom *c, OBAtom *d);

      //! \return true if atomtype has sbmb set in mmffprop.par

      bool HasSbmbSet(int atomtype);

      //! \return true if atomtype has pilp set in mmffprop.par

      bool HasPilpSet(int atomtype);

      //! \return true if atomtype has arom set in mmffprop.par

      bool HasAromSet(int atomtype);

      //! \return true if atomtype has lin set in mmffprop.par

      bool HasLinSet(int atomtype);

      //! \return the crd value for the atomtype in mmffprop.par

      int GetCrd(int atomtype);

      //! \return the val value for the atomtype in mmffprop.par

      int GetVal(int atomtype);

      //! \return the mltb value for the atomtype in mmffprop.par

      int GetMltb(int atomtype);

      //! \return the level 2 equivalent atom type for type (mmffdef.par)

      int EqLvl2(int type);

      //! \return the level 3 equivalent atom type for type (mmffdef.par)

      int EqLvl3(int type);

      //! \return the level 4 equivalent atom type for type (mmffdef.par)

      int EqLvl4(int type);

      //! \return the level 5 equivalent atom type for type (mmffdef.par)

      int EqLvl5(int type);

      //! \return the canonical bond index

      unsigned int GetCXB(int type, int a, int b);

      //! \return the canonical angle index

      unsigned int GetCXA(int type, int a, int b, int c);

      //! \return the canonical stretch-bend index

      unsigned int GetCXS(int type, int a, int b, int c);

      //! \return the canonical out-of-plane index

      unsigned int GetCXO(int a, int b, int c, int d);

      //! \return the canonical torsion index

      unsigned int GetCXT(int type, int a, int b, int c, int d);

      //! \return the canonical bond-charge-increment index

      unsigned int GetCXQ(int type, int a, int b);

      //! \return the U value for the atom from table X page 631

      double GetUParam(OBAtom* atom);

      //! \return the Z value for the atom from table VI page 628

      double GetZParam(OBAtom* atom);

      //! \return the C value for the atom from table VI page 628

      double GetCParam(OBAtom* atom);

      //! \return the V value for the atom from table X page 631

      double GetVParam(OBAtom* atom);

      //! return the covalent radius from Blom and Haaland, value from etab if not available

      double GetCovalentRadius(OBAtom* a);

      //! return the bond length calculated with a modified version of the Schomaker-Stevenson rule

      double GetRuleBondLength(OBAtom* a, OBAtom* b);

      //! return the bond length from mmffbond.par, if not found, one is calculated with a modified version of the Schomaker-Stevenson rule

      double GetBondLength(OBAtom* a, OBAtom* b);

      //! Same as OBForceField::GetParameter, but takes (bond/angle/torsion) type in account and takes 0 as wildcart.

      OBFFParameter* GetParameter1Atom(int a, std::vector<OBFFParameter> &parameter);

      OBFFParameter* GetParameter2Atom(int a, int b, std::vector<OBFFParameter> &parameter);

      OBFFParameter* GetParameter3Atom(int a, int b, int c, std::vector<OBFFParameter> &parameter);

      //! Same as OBForceField::GetParameter, but takes (bond/angle/torsion) type in account and takes 0 as wildcart.

      OBFFParameter* GetTypedParameter2Atom(int ffclass, int a, int b, std::vector<OBFFParameter> &parameter);

      OBFFParameter* GetTypedParameter3Atom(int ffclass, int a, int b, int c, std::vector<OBFFParameter> &parameter);

      OBFFParameter* GetTypedParameter4Atom(int ffclass, int a, int b, int c, int d, std::vector<OBFFParameter> &parameter);

      // OBFFParameter vectors to contain the parameters

      std::vector<OBFFParameter> _ffbondparams;

      std::vector<OBFFParameter> _ffbndkparams;

      std::vector<OBFFParameter> _ffangleparams;

      std::vector<OBFFParameter> _ffstrbndparams;

      std::vector<OBFFParameter> _ffdfsbparams;

      std::vector<OBFFParameter> _fftorsionparams;

      std::vector<OBFFParameter> _ffoopparams;

      std::vector<OBFFParameter> _ffvdwparams;

      std::vector<OBFFParameter> _ffchgparams;

      std::vector<OBFFParameter> _ffpbciparams;

      std::vector<OBFFParameter> _ffdefparams;

      std::vector<OBFFParameter> _ffpropparams;

      OBBitVec       _ffpropPilp;

      OBBitVec       _ffpropArom;

      OBBitVec       _ffpropLin;

      OBBitVec       _ffpropSbmb;

      // OBFFXXXCalculationYYY vectors to contain the calculations

      std::vector<OBFFBondCalculationMMFF94>          _bondcalculations;

      std::vector<OBFFAngleCalculationMMFF94>         _anglecalculations;

      std::vector<OBFFStrBndCalculationMMFF94>        _strbndcalculations;

      std::vector<OBFFTorsionCalculationMMFF94>       _torsioncalculations;

      std::vector<OBFFOOPCalculationMMFF94>           _oopcalculations;

      std::vector<OBFFVDWCalculationMMFF94>           _vdwcalculations;

      std::vector<OBFFElectrostaticCalculationMMFF94> _electrostaticcalculations;

      bool mmff94s;

    public:

      //! Constructor

      explicit OBForceFieldMMFF94(const char* ID, bool IsDefault=true) : OBForceField(ID, IsDefault)

      {

        _validSetup = false;

        _init = false;

        _rvdw = 7.0;

        _rele = 15.0;

        _epsilon = 1.0; // default electrostatics

        _pairfreq = 15;

        _cutoff = false;

        _linesearch = LineSearchType::Newton2Num;

        _gradientPtr = nullptr;

        _grad1 = nullptr;

  if (!strncmp(ID, "MMFF94s", 7)) {

          mmff94s = true;

          _parFile = std::string("mmff94s.ff");

  } else {

          mmff94s = false;

          _parFile = std::string("mmff94.ff");

  }

      }

      //! Destructor

      virtual ~OBForceFieldMMFF94();

      //! Assignment

      OBForceFieldMMFF94 &operator = (OBForceFieldMMFF94 &);

      //!Clone the current instance. May be desirable in multithreaded environments

      OBForceFieldMMFF94* MakeNewInstance() override

      {

        return new OBForceFieldMMFF94(_id, false);

      }

      //! Get the description for this force field

      const char* Description() override

      {

        if (mmff94s)

          return "MMFF94s force field.";

  else

          return "MMFF94 force field.";

      }

      //! Get the unit in which the energy is expressed

      std::string GetUnit() override

      {

        return std::string("kcal/mol");

      }

      //! \return that analytical gradients are implemented for MMFF94

      bool HasAnalyticalGradients() override { return true; }

      //! Returns total energy

      double Energy(bool gradients = true) override;

      //! Returns the bond stretching energy

      template<bool> double E_Bond();

      double E_Bond(bool gradients = true) override

      {

        return gradients ? E_Bond<true>() : E_Bond<false>();

      }

      //! Returns the angle bending energy

      template<bool> double E_Angle();

      double E_Angle(bool gradients = true) override

      {

        return gradients ? E_Angle<true>() : E_Angle<false>();

      }

      //! Returns the stretch-bend energy

      template<bool> double E_StrBnd();

      double E_StrBnd(bool gradients = true) override

      {

        return gradients ? E_StrBnd<true>() : E_StrBnd<false>();

      }

      //! Returns the torsional energy

      template<bool> double E_Torsion();

      double E_Torsion(bool gradients = true) override

      {

        return gradients ? E_Torsion<true>() : E_Torsion<false>();

      }

      //! Returns the out-of-plane bending energy

      template<bool> double E_OOP();

      double E_OOP(bool gradients = true) override

      {

        return gradients ? E_OOP<true>() : E_OOP<false>();

      }

      //! Returns the Van der Waals energy (Buckingham potential)

      template<bool> double E_VDW();

      double E_VDW(bool gradients = true) override

      {

        return gradients ? E_VDW<true>() : E_VDW<false>();

      }

      //! Returns the dipole-dipole interaction energy

      template<bool> double E_Electrostatic();

      double E_Electrostatic(bool gradients = true) override

      {

        return gradients ? E_Electrostatic<true>() : E_Electrostatic<false>();

      }

      //! Validate MMFF94 using validation suite

      bool Validate() override;

      //! Compare and print the numerical and analytical gradients

      bool ValidateGradients() override;

  }; // class OBForceFieldMM2

}// namespace OpenBabel

//! \file forcefieldmmff94.h

//! \brief MMFF94 force field