#include <openbabel/obconversion.h>
#include <openbabel/mol.h>
#include <iostream>


using namespace OpenBabel;
using namespace std;

int main(int argc,char **argv) {

  OBMol mymol;
  OBAtom *atom;
  OBAtom *hatom;
  OBConversion conv;
  OBFormat *inFormat;
  OBSmartsPattern smarts;

  string infilename;
  string smpat;

  vector< vector <int> > maplist;
  vector< vector <int> >::iterator matches;
  vector< OBBond* >::iterator ib;
  vector< int >::iterator j;

  int nat, counter;
  int nmatches;
  int test;
  double cx,cy,cz,dnc;
  bool notatend;

  mymol.Clear();

  if(argc!=7) {
    cout << endl << "usage: dips filename pattern number x y z" << endl << endl;
    cout << "       filename: molecular structures in any babel supported format" << endl;
    cout << "       pattern: smarts/smiles pattern in inverted commas" << endl;
    cout << "       number: index number of atom in pattern (if <0 run test)" << endl;
    cout << "       x y z: coordinates of anchor point" << endl << endl;
    exit(1);
  }
  infilename = argv[1];
  smpat      = argv[2];
  test       = atoi(argv[3]);

  smarts.Init(smpat);

  inFormat = conv.FormatFromExt(infilename.c_str());
  conv.SetInFormat(inFormat);

  cx=atof(argv[4]);
  cy=atof(argv[5]);
  cz=atof(argv[6]);

  notatend = conv.ReadFile(&mymol, infilename);

  while (notatend) {

    if (smarts.Match(mymol)) {

      maplist = smarts.GetUMapList();

      cout << mymol.GetTitle();

      nmatches=0;

      for (matches = maplist.begin(); matches != maplist.end(); matches++) {

        nmatches++;

        if(nmatches==1){

          if(test<0) {

            counter=0;

            for (j = (*matches).begin();j != (*matches).end();++j) {

              atom = mymol.GetAtom(*j);
              cout << " (" << counter << " " << atom->GetIdx() << " " << atom->GetType() << ")";
              counter++;
            }

          } else {

            atom = mymol.GetAtom((*matches)[test]);
            dnc = sqrt( (cx - atom->x())*(cx - atom->x()) + 
                        (cy - atom->y())*(cy - atom->y()) + 
                        (cz - atom->z())*(cz - atom->z()) );
            printf("%12.3f",dnc);

          }
        }
      }

      if(nmatches>1) cout << " ... more than one match!" << endl;
      else cout << endl;

    }

    mymol.Clear();
    notatend = conv.Read(&mymol);
  }
}