#include <iostream>
#include <fstream>
#include <string>
#include <cmath>

#include <openbabel/mol.h>
#include <openbabel/obconversion.h>
#include <openbabel/parsmart.h>
#include <openbabel/mcdlutil.h>

using namespace std;
using namespace OpenBabel;

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

  ifstream ifs_prot(argv[1]);
  if (ifs_prot.fail()) {
    cerr << "Failed to open the file: " << argv[1] << endl;
    return -1;
  }

  ifstream ifs_sdf(argv[2]);
  if (ifs_sdf.fail()) {
    cerr << "Failed to open the file: " << argv[2] << endl;
    return -1;
  }

  ofstream ofs(argv[3]);

  OBConversion conv1;
  conv1.SetInFormat("PDB");
  OBConversion conv2;
  conv2.SetInAndOutFormats("SDF", "SDF");

  OBMol prot;
  conv1.Read(&prot, &ifs_prot);
  OBAtom* zinc;
  for (auto ii = prot.BeginResidues(); ii != prot.EndResidues(); ii++) {
    if ((*ii)->GetNum() == 2036) {
      zinc = *((*ii)->BeginAtoms());
    }
  }

  OBSmartsPattern smp;
  smp.Init("[OX1]NC(=O)[#6]");

  OBMol mol;
  while (conv2.Read(&mol, &ifs_sdf)) {
    smp.Match(mol);
    vector<vector<int> > mlist = smp.GetUMapList();
    double gx = 0;
    double gy = 0;
    double gz = 0;
    for (int i = 0; i < mlist[0].size() - 1; i++) {
      OBAtom* a = mol.GetAtom(mlist[0][i]);
      gx += a->x();
      gy += a->y();
      gz += a->z();
    }
    gx /= (double)(mlist[0].size() - 1);
    gy /= (double)(mlist[0].size() - 1);
    gz /= (double)(mlist[0].size() - 1);

    double d = (zinc->x() - gx) * (zinc->x() - gx) + (zinc->y() - gy) * (zinc->y() - gy) +
               (zinc->z() - gz) * (zinc->z() - gz);
    d = sqrt(d);

    if (d < 3.5) {
      conv2.Write(&mol, &ofs);
    }
  }
}