/////////////////////////////////////////////////////////////////////////////
/// molgeom.cpp - calculate some molecular geometry.
///    - bond lengths
///    - bond angles
///    - torsion angles
///    - "ring pucker", i.e. non-planarity (avg abs torsion angle)
///
/// Author: Jeremy Yang
///  10 Jul 2003
/////////////////////////////////////////////////////////////////////////////
#include "openeye.h"
#include <iostream>
#include <vector>
#include <math.h>
#include "oesystem.h"
#include "oechem.h"

using namespace OESystem;
using namespace OEChem;
using namespace std;

static oemolistream ifs;
static int ifmt=0;
static int wegot_ifile=0;
static int verbose=0;
static int minring=5;
static int maxring=7;


string Sprintformats();   /// from formats.cpp
int ParseCommandLine(int argc,char *argv[]);
void findcenter(OEMolBase *m,double xyzavg[3]);
double *rmat(char axis,double a);

void help(string msg)
{
cerr << msg << endl
     <<"\t************************************************************"<<endl
     <<"\t| molgeom [options] [<infile]                              |"<<endl
     <<"\t|   options:                                               |"<<endl
     <<"\t|     -i INFILE                                            |"<<endl
     <<"\t|     -ifmt FMT (lowcase, e.g., smi,mdl)                   |"<<endl
     <<"\t|     -minring MINRING [5]                                 |"<<endl
     <<"\t|     -maxring MAXRING [7]                                 |"<<endl
     <<"\t|     -v .......  verbose                                  |"<<endl
     <<"\t|     -help                                                |"<<endl
     <<"\t************************************************************"<<endl;
  cerr << endl << Sprintformats() << endl;
  exit(1);
}


int main(int argc,char *argv[])
{
  OEIter<OEMolBase> m;
  OEGraphMol matchmol;                // matched submol
  OESubSearch pat;
  OEIter<OEMatchBase> match;
  OEIter<OEMatchPair<OEAtomBase> > mp;
  OEAtomBase *atom1=NULL,*atom2=NULL;
  OEAtomBase *ablen[2];
  OEAtomBase *abang[3];
  OEAtomBase *ator[4];
  OEAtomBase *aring[7];
  OEIter<OEAtomBase> a1,a2;
  OEIter<OEBondBase> b;
  string smi;
  string smarts;
  double *rot;
  double xyz[3];
  double maxl2,l2;
  double t;
  double sumtor;
  int dim,i,j;
  int k;

  ParseCommandLine(argc,argv);

  if (!wegot_ifile) {
    if (ifmt) ifs.open();
    else help("Input fmt required w/ stdin.");
  }

  for (m=ifs.GetMolBases();m;++m)
  {
    OECreateCanSmiString(smi,m);
    cout << "cansmi: " << smi << endl;
    OETriposAtomNames(m);
    dim=m->GetDimension();
    cout << "formal dimension of mol: " << dim << endl;
    OESetDimensionFromCoords(m);
    dim=m->GetDimension();
    cout << "perceived dimension of mol: " << dim << endl;
    if (dim<1) 
    {
      cerr << "No coordinates -- skipping mol." << endl;
      continue;
    }

    /// Center mol (avg position), and verify.
    findcenter(m,xyz);
    printf("center: %.2f,%.2f,%.2f\n",xyz[0],xyz[1],xyz[2]);
    OECenter(m,xyz);  // OECenter() centers and fills xyz.
    printf("center: %.2f,%.2f,%.2f\n",xyz[0],xyz[1],xyz[2]);
    findcenter(m,xyz);

    /// Now move back to original position.
    OETranslate(m,xyz);
    findcenter(m,xyz);
    printf("center: %.2f,%.2f,%.2f\n",xyz[0],xyz[1],xyz[2]);

    /// Find max atom-atom length.
    maxl2=0.0;
    for(a1=m->GetAtoms();a1;++a1)
    {
      for(a2=m->GetAtoms();a2;++a2)
      {
         if (a1->GetIdx()>=a2->GetIdx()) continue;
         l2=OEGetDistance2(m,a1,a2);
         if (l2>maxl2)
         {
           maxl2=l2;
           atom1=a1;
           atom2=a2;
         }
       }
    }
    m->GetCoords(atom1,xyz);
    printf("%s: %.2f,%.2f,%.2f\n",atom1->GetName(),xyz[0],
      xyz[1],xyz[2]);
    m->GetCoords(atom2,xyz);
    printf("%s: %.2f,%.2f,%.2f\n",atom2->GetName(),xyz[0],
      xyz[1],xyz[2]);
    printf("length of mol (%s <-> %s): %.2f\n",
      atom1->GetName(),atom2->GetName(),sqrt(maxl2));

    /// Rotate molecule.
    rot = rmat('z',M_PI/2.0);  // 90degrees, CW, z-axis
    OERotate(m,rot);
    m->GetCoords(atom1,xyz);
    printf("%s: %.2f,%.2f,%.2f\n",atom1->GetName(),xyz[0],
      xyz[1],xyz[2]);
    m->GetCoords(atom2,xyz);
    printf("%s: %.2f,%.2f,%.2f\n",atom2->GetName(),xyz[0],
      xyz[1],xyz[2]);
    printf("length of mol (%s <-> %s): %.2f\n",
      atom1->GetName(),atom2->GetName(),sqrt(maxl2));

    /// Find bond lengths.
    pat.Init("*~*");  // smarts for any 2 bonded atoms
    for (i=0,match=pat.Match(m,1);match;++match)
    {
      for (mp=match->GetAtoms(),j=0;mp;++mp)
        ablen[j++]=mp->target;
      if (j!=2) cerr << "ERROR: bad bondlength atom count." << endl;
      OESubsetMol(matchmol,match);
      OEAssignMDLHydrogens(matchmol);
      OECreateCanSmiString(smi,matchmol);
      printf("%d.  bond length: %s  %.2f\n",++i,smi.c_str(),
        OEGetDistance(m,ablen[0],ablen[1]));
      matchmol.Clear();
    }

    /// Find bond angles.
    pat.Init("*~*~*");  // smarts for any 3 bonded atoms
    for (i=0,match=pat.Match(m,1);match;++match)
    {
      for (mp=match->GetAtoms(),j=0;mp;++mp)
        abang[j++]=mp->target;
      if (j!=3) cerr << "ERROR: bad bondangle atom count." << endl;
      OESubsetMol(matchmol,match);
      OEAssignMDLHydrogens(matchmol);
      OECreateCanSmiString(smi,matchmol);
      printf("%d.  bond angle: %s  %.1f\n",++i,smi.c_str(),
        OEGetAngle(m,abang[0],abang[1],abang[2])*OEMath::Rad2Deg);
      matchmol.Clear();
    }

    /// Find torsion angles.
    pat.Init("*~[!$(*#*)&!D1]-!@[!$(*#*)&!D1]~*");
    for (i=0,match=pat.Match(m,1);match;++match)
    {
      for (mp=match->GetAtoms(),j=0;mp;++mp)
        ator[j++]=mp->target;
      if (j!=4) cerr << "ERROR: bad torsion atom count." << endl;
      OESubsetMol(matchmol,match);
      OEAssignMDLHydrogens(matchmol);
      OECreateCanSmiString(smi,matchmol);
      printf("%d.  torsion angle: %s  %.1f\n",++i,smi.c_str(),
           OEGetTorsion(m,ator[0],ator[1],ator[2],ator[3])*OEMath::Rad2Deg);
      matchmol.Clear();
    }

    /// Calculate pucker (avg absolute torsion angle) for rings.
    for (i=minring,j=0;i<=maxring;++i)
    {
      smarts="*~1";
      for (k=1;k<i;++k) smarts+="~*";
      smarts+="1";
      pat.Init(smarts.c_str());
      for (match=pat.Match(m,1);match;++match)
      {
        OESubsetMol(matchmol,match);
        OEAssignMDLHydrogens(matchmol);
        OECreateCanSmiString(smi,matchmol);
	printf("%d.  %d member ring: %s",++j,i,smi.c_str());
        for (mp=match->GetAtoms(),k=0;mp;++mp)
          aring[k++]=mp->target;
	sumtor=0.0;
	for (k=0;k<i;++k)
	{
          for (int kk=0;kk<4;++kk) ator[kk]=aring[(k+kk)%i];
	  t=OEGetAbsTorsion(m,ator[0],ator[1],ator[2],ator[3]);
          sumtor+=t;
	}
        printf(" avg torsion (pucker): %.3f\n",(sumtor/i)*OEMath::Rad2Deg);
        matchmol.Clear();
      }
    }
  }
  return 0;
}

/////////////////////////////////////////////////////////////////////////////
int sameword(char *p, char *q)
{ while (*p||*q) {if (((*p++)-(*q++))%('a'-'A')) return 0;} return 1; }

/////////////////////////////////////////////////////////////////////////////
int ParseCommandLine(int argc,char *argv[])
{
  for (++argv; --argc; ++argv)
  {
    if (sameword(*argv,"-i"))
    {
    if (--argc==0)
      help("ERROR: -i requires value.");
      if (!ifs.open(*++argv))
        OEThrow.Error("Unable to open \"%s\"",*argv);
      wegot_ifile=1;
    }
    else if (sameword(*argv,"-ifmt"))
    {
      if (--argc==0) help("ERROR: -ifmt requires value.");
      if (!(ifmt=OEGetFileType(*++argv)))
        OEThrow.Error("Unknown file type \"%s\"",*argv);
    }
    else if (sameword(*argv,"-minring"))
    {
      if (--argc==0) help("ERROR: -minring requires value.");
      minring=atoi(*++argv);
    }
    else if (sameword(*argv,"-maxring"))
    {
      if (--argc==0) help("ERROR: -maxring requires value.");
      maxring=atoi(*++argv);
    }
    else if (sameword(*argv,"-v"))
      verbose=1;
    else if (sameword(*argv,"-help"))
      help("");
    else
      help("ERROR: Bad option");
  }
  return 1;
}
  
/////////////////////////////////////////////////////////////////////////////
void findcenter(OEMolBase *m,double xyzavg[3])
{
  double xyz[3],xyzsum[3]={0.0,0.0,0.0};
  OEIter<OEAtomBase> a;
  for(a=m->GetAtoms();a;++a)
  {
    m->GetCoords(a,xyz);
    xyzsum[0]+=xyz[0];
    xyzsum[1]+=xyz[1];
    xyzsum[2]+=xyz[2];
  }
  xyzavg[0]=xyzsum[0]/m->NumAtoms();
  xyzavg[1]=xyzsum[1]/m->NumAtoms();
  xyzavg[2]=xyzsum[2]/m->NumAtoms();
  return;
}

/////////////////////////////////////////////////////////////////////////////
///  rmat - returns axial rotation matrices (clockwise looking
///  into axis).
/// 
/// Refs: http://mathworld.wolfram.com/RotationMatrix.html,
/// Goldstein 1980, pp. 146-147 and 608; Arfken 1985, pp. 199-200
/////////////////////////////////////////////////////////////////////////////
double *rmat(char axis,double a)
{
  static double mat[9];

  switch (axis) 
  {
    case 'x':
      mat[0]=1.0;      mat[1]=1.0;       mat[2]=0.0;
      mat[3]=0.0;      mat[4]=cos(a);    mat[5]=sin(a);
      mat[6]=0.0;      mat[7]=-sin(a);   mat[8]=cos(a);
      break;
    case 'y':
      mat[0]=cos(a);   mat[1]=0.0;       mat[2]=-sin(a);
      mat[3]=0.0;      mat[4]=1.0;       mat[5]=0.0;    
      mat[6]=sin(a);   mat[7]=0.0;       mat[8]=cos(a);
      break;
    default:
      mat[0]=cos(a);   mat[1]=sin(a);    mat[2]=0.0;    
      mat[3]=-sin(a);  mat[4]=cos(a);    mat[5]=0.0;    
      mat[6]=0.0;      mat[7]=0.0;       mat[8]=1.0;   
  }
  return mat;
}