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

Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.

Some portions Copyright (C) 2001-2007 by Geoffrey R. Hutchison

Some portions Copyright (C) 2004 by Chris Morley

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 <openbabel/babelconfig.h>

#include <openbabel/obmolecformat.h>

#include <openbabel/mol.h>

#include <openbabel/atom.h>

#include <openbabel/bond.h>

#include <openbabel/obiter.h>

#include <openbabel/elements.h>

#include <openbabel/generic.h>

#include <openbabel/kekulize.h>

#include <openbabel/obfunctions.h>

#include <openbabel/data.h>

#include <cstdlib>

using namespace std;

namespace OpenBabel

{

  //The routine WriteSmiOrderedMol2() in the original mol2.cpp is presumably

  //another output format, but was not made available in version 100.1.2, nor

  //is it here.

  class MOL2Format : public OBMoleculeFormat

  {

  public:

    //Register this format type ID

    MOL2Format()

    {

      OBConversion::RegisterFormat("mol2",this, "chemical/x-mol2");

      OBConversion::RegisterFormat("ml2",this);

      OBConversion::RegisterFormat("sy2",this);

      OBConversion::RegisterOptionParam("c", this, 0, OBConversion::INOPTIONS);

      OBConversion::RegisterOptionParam("c", this, 0, OBConversion::OUTOPTIONS);

      OBConversion::RegisterOptionParam("l", this, 0, OBConversion::OUTOPTIONS);

      OBConversion::RegisterOptionParam("u", this, 0, OBConversion::OUTOPTIONS);

    }

    const char* Description() override  // required

    {

      return

        "Sybyl Mol2 format\n"

        "Read Options e.g. -ac\n"

        "  c               Read UCSF Dock scores saved in comments preceding molecules\n\n"

        "Write Options e.g. -xl\n"

        "  l               Output ignores residue information (only ligands)\n"

        "  c               Write UCSF Dock scores saved in comments preceding molecules\n"

        "  u               Do not write formal charge information in UNITY records\n\n";

    }

    const char* SpecificationURL() override

    {

      return "http://www.tripos.com/data/support/mol2.pdf"; // XXX dead

    }

    const char* GetMIMEType() override

    { return "chemical/x-mol2"; }

    int SkipObjects(int n, OBConversion* pConv) override;

    //*** This section identical for most OBMol conversions ***

    ////////////////////////////////////////////////////

    /// The "API" interface functions

    bool ReadMolecule(OBBase* pOb, OBConversion* pConv) override;

    bool WriteMolecule(OBBase* pOb, OBConversion* pConv) override;

  };

  //***

  //Make an instance of the format class

  MOL2Format theMOL2Format;

  // Helper function for ReadMolecule

  // \return Is this atom a sulfur in a (di)thiocarboxyl (-CS2, -COS, CS2H or COSH) group?

  static bool IsThiocarboxylSulfur(OBAtom* queryatom)

  {

    if (queryatom->GetAtomicNum() != OBElements::Sulfur)

      return(false);

    if (queryatom->GetHvyDegree() != 1)

      return(false);

    OBAtom *atom = nullptr;

    OBBond *bond;

    OBBondIterator i;

    for (bond = queryatom->BeginBond(i); bond; bond = queryatom->NextBond(i))

      if ((bond->GetNbrAtom(queryatom))->GetAtomicNum() == OBElements::Carbon)

      {

        atom = bond->GetNbrAtom(queryatom);

        break;

      }

    if (!atom)

      return(false);

    if (!(atom->CountFreeSulfurs() == 2)

      && !(atom->CountFreeOxygens() == 1 && atom->CountFreeSulfurs() == 1))

      return(false);

    //atom is connected to a carbon that has a total

    //of 2 attached free sulfurs or 1 free oxygen and 1 free sulfur

    return(true);

  }

  static bool IsOxygenOrSulfur(OBAtom *atom)

  {

    switch (atom->GetAtomicNum()) {

    case 8: case 16: return true;

    default: return false;

    }

  }

  static unsigned int GetAtomicNumAndIsotope(const char* symbol, int *isotope)

  {

    const char* p = symbol;

    switch (p[0]) {

    case 'D':

      if (p[1] == '\0') {

        *isotope = 2;

        return 1;

      }

      break;

    case 'T':

      if (p[1] == '\0') {

        *isotope = 3;

        return 1;

      }

      break;

    }

    return OBElements::GetAtomicNum(symbol);

  }

  //read from ifs until next rti is found and return it

  static string read_until_rti(istream & ifs) 

  {

      char buffer[BUFF_SIZE];

      for (;;)

      {

        if (!ifs.getline(buffer,BUFF_SIZE))

          return "";

        if (!strncmp(buffer,"@<TRIPOS>",9))

          return string(buffer);

      }      

  }

  /////////////////////////////////////////////////////////////////

  bool MOL2Format::ReadMolecule(OBBase* pOb, OBConversion* pConv)

  {

    OBMol* pmol = pOb->CastAndClear<OBMol>();

    if (pmol == nullptr)

      return false;

    //Define some references so we can use the old parameter names

    istream &ifs = *pConv->GetInStream();

    OBMol &mol = *pmol;

    //Old code follows...

    bool foundAtomLine = false;

    char buffer[BUFF_SIZE];

    char *comment = nullptr;

    string str,str1;

    vector<string> vstr;

    int len;

    // Prevent reperception

    mol.SetChainsPerceived();

    mol.BeginModify();

    for (;;)

      {

        if (!ifs.getline(buffer,BUFF_SIZE))

          return(false);

        if (pConv->IsOption("c", OBConversion::INOPTIONS) != nullptr && EQn(buffer, "###########", 10))

          {

            char attr[32], val[32];

            // CVE-2022-43607: width-limit both specifiers to the buffer size.

            sscanf(buffer, "########## %31[^:]:%31s", attr, val);

            OBPairData *dd = new OBPairData;

            dd->SetAttribute(attr);

            dd->SetValue(val);

            dd->SetOrigin(fileformatInput);

            mol.SetData(dd);

          }

        if (EQn(buffer,"@<TRIPOS>MOLECULE",17))

          break;

      }

    // OK, just read MOLECULE line

    int lcount;

    int natoms = 0, nbonds = 0;

    bool hasPartialCharges = true;

    for (lcount=0;;lcount++)

      {

        if (!ifs.getline(buffer,BUFF_SIZE))

          return(false);

        if (EQn(buffer,"@<TRIPOS>ATOM",13))

          {

            foundAtomLine = true;

            break;

          }

        if (lcount == 0)

          {

            tokenize(vstr,buffer);

            if (!vstr.empty())

              mol.SetTitle(buffer);

          }

        else if (lcount == 1)

          sscanf(buffer,"%d%d",&natoms,&nbonds);

        else if (lcount == 3) // charge descriptions

          {

            // Annotate origin of partial charges

            OBPairData *dp = new OBPairData;

            dp->SetAttribute("PartialCharges");

            dp->SetValue(buffer);

            dp->SetOrigin(fileformatInput);

            mol.SetData(dp);

            if (strncasecmp(buffer, "NO_CHARGES", 10) == 0)

              hasPartialCharges = false;

          }

        else if (lcount == 4) //energy (?)

          {

            tokenize(vstr,buffer);

            if (!vstr.empty() && vstr.size() == 3)

              if (vstr[0] == "Energy")

                mol.SetEnergy(atof(vstr[2].c_str()));

          }

        else if (lcount == 5) //comment

          {

            if ( buffer[0] )

              {

                len = (int) strlen(buffer)+1;

                //! @todo allow better multi-line comments

                // which don't allow ill-formed data to consume memory

                // Thanks to Andrew Dalke for the pointer

                delete [] comment;

                comment = new char [len];

                memcpy(comment,buffer,len);

              }

          }

      }

    if (!foundAtomLine)

      {

        mol.EndModify();

        mol.Clear();

        obErrorLog.ThrowError(__FUNCTION__, "Unable to read Mol2 format file. No atoms found.", obWarning);

        return(false);

      }

    if (natoms < 0 || natoms >= 100000000 || nbonds < 0 || nbonds >= 100000000)

      {

        mol.EndModify();

        mol.Clear();

        obErrorLog.ThrowError(__FUNCTION__, "Unable to read Mol2 format file. Invalid atom or bond count.", obWarning);

        return(false);

      }

    mol.ReserveAtoms(natoms);

    int i;

    vector3 v;

    OBAtom atom;

    double x,y,z,pcharge;

    char temp_type[BUFF_SIZE], resname[BUFF_SIZE], atmid[BUFF_SIZE];

    int elemno, resnum = -1;

    int isotope = 0;

    bool has_explicit_hydrogen = false;

    bool has_residue_information = false;

    ttab.SetFromType("SYB");

    for (i = 0;i < natoms;i++)

      {

        if (!ifs.getline(buffer,BUFF_SIZE))

          return(false);

        // Required ATOM columns are the atom name, x/y/z and the SYBYL atom

        // type.  subst_id, subst_name and charge are optional, so reset them

        // to their defaults each iteration (these locals are declared outside

        // the loop and would otherwise carry over from the previous atom).

        resnum = -1;

        resname[0] = '\0';

        pcharge = 0.0;

        isotope = 0;

        if (sscanf(buffer," %*s %1024s %lf %lf %lf %1024s %d %1024s %lf",

               atmid, &x,&y,&z, temp_type, &resnum, resname, &pcharge) < 5)

          {

            obErrorLog.ThrowError(__FUNCTION__,

                                  "Unable to read Mol2 format file. Truncated "

                                  "atom record.", obWarning);

            return(false);

          }

        atom.SetVector(x, y, z);

        atom.SetFormalCharge(0);

        // Handle "CL" and "BR" and other mis-typed atoms

        str = temp_type;

        if (strncmp(temp_type, "CL", 2) == 0) {

          str = "Cl";

        } else  if (strncmp(temp_type,"BR",2) == 0) {

          str = "Br";

        } else if (strncmp(temp_type,"S.o2", 4) == 0) {

          str = "S.O2";

        } else if (strncmp(temp_type,"S.o", 3) == 0) {

          str = "S.O";

        } else if (strncmp(temp_type,"SI", 2) == 0) {

          str = "Si";

        // The following cases are entries which are not in openbabel/data/types.txt

        // and should probably be added there

        } else if (strncmp(temp_type,"S.1", 3) == 0) {

          str = "S.2"; // no idea what the best type might be here

        } else if (strncmp(temp_type,"P.", 2) == 0) {

          str = "P.3";

        } else if (strncasecmp(temp_type,"Ti.", 3) == 0) { // e.g. Ti.th

          str = "Ti";

        } else if (strncasecmp(temp_type,"Ru.", 3) == 0) { // e.g. Ru.oh

          str = "Ru";

        // Fixes PR#3557898

        } else if (strncmp(temp_type, "N.4", 3) == 0) {

          atom.SetFormalCharge(1);

        }

        ttab.SetToType("ATN");

        ttab.Translate(str1,str);

        elemno = atoi(str1.c_str());

        ttab.SetToType("IDX");

        // We might have missed some SI or FE type things above, so here's

        // another check

        if( !elemno && isupper(temp_type[1]) )

          {

            temp_type[1] = (char)tolower(temp_type[1]);

            str = temp_type;

            ttab.SetToType("ATN");

            ttab.Translate(str1,str);

            elemno = atoi(str1.c_str());

            ttab.SetToType("IDX");

          }

        // One last check if there isn't a period in the type,

        // it's a malformed atom type, but it may be the element symbol

        // GaussView does this (PR#1739905)

        if ( !elemno ) {

          // check if it's "Du" or "Xx" and the element is in the atom name

          if (str == "Du" || str == "Xx") {

            str = atmid;

            for (unsigned int i = 0; i < str.length(); ++i)

              if (!isalpha(str[i])) {

                str.erase(i);

                break; // we've erased the end of the string

              }

          }

          std::stringstream errorMsg;

          errorMsg << "This Mol2 file is non-standard. Problem with molecule: "

                   << mol.GetTitle()

                   << " Cannot interpret atom types correctly, instead attempting to interpret atom type: "

                   << str << " as elements instead.";

          obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obWarning);

          string::size_type dotPos = str.find('.');

          if (dotPos == string::npos) {

            elemno = GetAtomicNumAndIsotope(str.c_str(), &isotope);

          }

        }

        atom.SetAtomicNum(elemno);

        if (isotope)

          atom.SetIsotope(isotope);

        else if (elemno == 1)

          has_explicit_hydrogen = true;

        ttab.SetToType("INT");

        ttab.Translate(str1,str);

        atom.SetType(str1);

        atom.SetPartialCharge(pcharge);

        // MMFF94 has different atom types for Cu(I) and Cu(II)

        // as well as for Fe(II) and Fe(III), so the correct formal

        // charge is needed for correct atom type assignment

        if (str1 == "Cu" || str1 == "Fe")

          atom.SetFormalCharge((int)pcharge);

        if (!mol.AddAtom(atom))

          return(false);

        if (fabs(pcharge) >= 2e-6)

          hasPartialCharges = true;

        // Add residue information if it exists

        if (resnum != -1 && resnum != 0 &&

            strlen(resname) != 0 && strncmp(resname,"<1>", 3) != 0)

          {

            has_residue_information = true;

            OBResidue *res  = (mol.NumResidues() > 0) ?

              mol.GetResidue(mol.NumResidues()-1) : nullptr;

            if (res == nullptr || res->GetName() != resname ||

                res->GetNum() != resnum)

              {

                vector<OBResidue*>::iterator ri;

                for (res = mol.BeginResidue(ri) ; res ; res = mol.NextResidue(ri))

                  if (res->GetName() == resname &&

                      res->GetNum() == resnum)

                    break;

                if (res == nullptr)

                  {

                    res = mol.NewResidue();

                    res->SetName(resname);

                    res->SetNum(resnum);

                  }

              }

            OBAtom *atomPtr = mol.GetAtom(mol.NumAtoms());

            res->AddAtom(atomPtr);

            res->SetAtomID(atomPtr, atmid);

          } // end adding residue info

      }

    string nextrti;

    do { nextrti = read_until_rti(ifs); } 

    while(nextrti != "@<TRIPOS>UNITY_ATOM_ATTR" && nextrti != "@<TRIPOS>BOND" && nextrti.length() > 0);

    if(nextrti == "@<TRIPOS>UNITY_ATOM_ATTR")

    { //read in formal charge information, must be done before Kekulization

        int aid = 0, num = 0;

        while (ifs.peek() != '@' && ifs.getline(buffer,BUFF_SIZE))

        {

          // Reset before each parse so a malformed header line cannot reuse

          // the previous entry's atom id / attribute count.

          aid = 0;

          num = 0;

          if (sscanf(buffer,"%d %d",&aid, &num) < 2)

            continue;

          for(int i = 0; i < num; i++)

          {

            if (!ifs.getline(buffer,BUFF_SIZE))

              return(false);

            if(strncmp(buffer, "charge", 6) == 0)

            {

              int charge = 0;

              sscanf(buffer,"%*s %d",&charge);

              if(aid >= 1 && aid <= (int)mol.NumAtoms())

              {

                OBAtom *atom = mol.GetAtom(aid);

                if (atom != nullptr)

                  atom->SetFormalCharge(charge);

              }

            }

          }

        }

    }

    while(nextrti != "@<TRIPOS>BOND" && nextrti.length() > 0)

      nextrti = read_until_rti(ifs);

    if(nextrti != "@<TRIPOS>BOND")

      return false;

    int start, end;

    bool needs_kekulization = false;

    for (i = 0; i < nbonds; i++)

      {

        if (!ifs.getline(buffer,BUFF_SIZE))

          return(false);

        // origin/target atom ids and the bond type are all required.

        temp_type[0] = '\0';

        if (sscanf(buffer,"%*d %d %d %1024s",&start,&end,temp_type) < 3)

          {

            obErrorLog.ThrowError(__FUNCTION__,

                                  "Unable to read Mol2 format file. Truncated "

                                  "bond record.", obWarning);

            return(false);

          }

        str = temp_type;

        unsigned int flags = 0;

        int order;

        if (str == "ar" || str == "AR" || str == "Ar") {

          order = 1;

          flags = OB_AROMATIC_BOND;

          needs_kekulization = true;

        }

        else if (str == "AM" || str == "am" || str == "Am")

          order = 1;

        else

          order = atoi(str.c_str());

        mol.AddBond(start, end, order, flags);

      }

    // TODO: Add a test case for the statement below of Paolo Tosco

    //       - I am currently assuming that is not a problem for the

    //         the current kekulization code, but it needs to be

    //         checked

    //   "Make a pass to ensure that there are no double bonds

    //    between atoms which are also involved in aromatic bonds

    //    as that may ill-condition kekulization (fixes potential

    //    issues with molecules like CEWYIM30 (MMFF94 validation suite)"

    mol.SetAromaticPerceived(); // don't trigger reperception

    if (has_explicit_hydrogen) {

      FOR_ATOMS_OF_MOL(atom, mol) {

        unsigned int total_valence = atom->GetTotalDegree();

        switch (atom->GetAtomicNum()) {

        case 8:

          if (total_valence != 1) continue;

          if (strcmp(atom->GetType(), "O2") != 0) continue; // TODO: the O.co2 type is lost by this point

          {

            OBAtomBondIter bit(&*atom);

            if (!bit->IsAromatic() && bit->GetBondOrder() == 1)

              atom->SetFormalCharge(-1); // set -1 charge on dangling O.co2

          }

          break;

        case 17: // Cl

          if (total_valence == 0)

            atom->SetFormalCharge(-1);

          break;

        }

      }

    }

    // Kekulization is necessary if an aromatic bond is present

    if (needs_kekulization) {

      // "de-aromatize" carboxylates and (di)thiocarboxylates

      // The typical case (in our test suite anyway) is a carboxylate binding to

      // a metal ion. The two O's have charges of -0.5 and the bond orders are aromatic

      FOR_ATOMS_OF_MOL(atom, mol) {

        OBAtom* oxygenOrSulfur = &*atom;

        // Look first for a terminal O/S

        if (!IsOxygenOrSulfur(oxygenOrSulfur) || oxygenOrSulfur->GetTotalDegree() != 1) continue;

        OBAtomBondIter bitA(oxygenOrSulfur);

        OBBond *bondA = &*bitA;

        if (!bondA->IsAromatic()) continue;

        // Look for the carbon

        OBAtom *carbon = bondA->GetNbrAtom(oxygenOrSulfur);

        if (carbon->GetAtomicNum() != 6) continue;

        // Look for the other oxygen or sulfur

        OBAtom* otherOxygenOrSulfur = nullptr;

        OBBond* bondB = nullptr;

        FOR_BONDS_OF_ATOM(bitB, carbon) {

          if (&*bitB == bondA || !bitB->IsAromatic()) continue;

          OBAtom* nbr = bitB->GetNbrAtom(carbon);

          if (IsOxygenOrSulfur(nbr) && nbr->GetTotalDegree() == 1) {

            otherOxygenOrSulfur = nbr;

            bondB = &*bitB;

          }

        }

        if (!otherOxygenOrSulfur) continue;

        if(otherOxygenOrSulfur->GetFormalCharge() != 0) continue; //formal charge already set on one

        // Now set as C(=O)O

        bondA->SetAromatic(false);

        oxygenOrSulfur->SetFormalCharge(-1);

        bondB->SetAromatic(false);

        bondB->SetBondOrder(2);

      }

      // First of all, set the atoms at the ends of the aromatic bonds to also

      // be aromatic. This information is required for OBKekulize.

      FOR_BONDS_OF_MOL(bond, mol) {

        if (bond->IsAromatic()) {

          bond->GetBeginAtom()->SetAromatic();

          bond->GetEndAtom()->SetAromatic();

        }

      }

      bool ok = OBKekulize(&mol);

      if (!ok) {

        stringstream errorMsg;

        errorMsg << "Failed to kekulize aromatic bonds in MOL2 file";

        std::string title = mol.GetTitle();

        if (!title.empty())

          errorMsg << " (title is " << title << ")";

        errorMsg << endl;

        obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obWarning);

        // return false; Should we return false for a kekulization failure?

      }

    }

    mol.EndModify();

    // Suggestion by Liu Zhiguo 2008-01-26

    // Mol2 files define atom types -- there is no need to re-perceive

    mol.SetAtomTypesPerceived();

    if (has_residue_information)

      mol.SetChainsPerceived();

    if (!has_explicit_hydrogen) {

      // Guess how many hydrogens are present on each atom based on typical valencies

      // TODO: implement the MOL2 valence model (if it exists)

      FOR_ATOMS_OF_MOL(matom, mol) {

        if (matom->GetImplicitHCount() == 0)

          OBAtomAssignTypicalImplicitHydrogens(&*matom);

      }

    }

    //must add generic data after end modify - otherwise it will be blown away

    if (comment)

      {

        OBCommentData *cd = new OBCommentData;

        cd->SetData(comment);

        cd->SetOrigin(fileformatInput);

        mol.SetData(cd);

        delete [] comment;

        comment = nullptr;

      }

    if (hasPartialCharges)

      mol.SetPartialChargesPerceived();

    /* Disabled due to PR#3048758 -- seekg is very slow with gzipped mol2

    // continue untill EOF or untill next molecule record

    streampos pos;

    for(;;)

      {

        pos = ifs.tellg();

        if (!ifs.getline(buffer,BUFF_SIZE))

          break;

        if (EQn(buffer,"@<TRIPOS>MOLECULE",17))

          break;

      }

    ifs.seekg(pos); // go back to the end of the molecule

    */

    return(true);

  }

  ////////////////////////////////////////////////////////////////

  bool MOL2Format::WriteMolecule(OBBase* pOb, OBConversion* pConv)

  {

    OBMol* pmol = dynamic_cast<OBMol*>(pOb);

    if (pmol == nullptr)

      return false;

    //Define some references so we can use the old parameter names

    ostream &ofs = *pConv->GetOutStream();

    OBMol &mol = *pmol;

    bool ligandsOnly = pConv->IsOption("l", OBConversion::OUTOPTIONS) != nullptr;

    bool skipFormalCharge = pConv->IsOption("u", OBConversion::OUTOPTIONS) != nullptr;

    //The old code follows....

    string str,str1;

    char buffer[BUFF_SIZE],label[BUFF_SIZE];

    char rnum[BUFF_SIZE],rlabel[BUFF_SIZE];

    //Check if UCSF Dock style coments are on

    if (pConv->IsOption("c", OBConversion::OUTOPTIONS) != nullptr) {

        vector<OBGenericData*>::iterator k;

        vector<OBGenericData*> vdata = mol.GetData();

        ofs << endl;

        for (k = vdata.begin();k != vdata.end();++k) {

            if ((*k)->GetDataType() == OBGenericDataType::PairData

            && (*k)->GetOrigin()!=local //internal OBPairData is not written

            && (*k)->GetAttribute()!="PartialCharges")

            {

                ofs << "##########\t" << (*k)->GetAttribute() << ":\t" << ((OBPairData*)(*k))->GetValue() << endl;

            }

        }

        ofs << endl;

    }

    ofs << "@<TRIPOS>MOLECULE" << endl;

    str = mol.GetTitle();

    if (str.empty())

      ofs << "*****" << endl;

    else

      ofs << str << endl;

    snprintf(buffer, BUFF_SIZE," %d %d 0 0 0", mol.NumAtoms(),mol.NumBonds());

    ofs << buffer << endl;

    ofs << "SMALL" << endl; // TODO: detect if we have protein, biopolymer, etc.

    OBPairData *dp = (OBPairData*)mol.GetData("PartialCharges");

    if (dp != nullptr) {

        // Tripos spec says:

        // NO_CHARGES, DEL_RE, GASTEIGER, GAST_HUCK, HUCKEL, PULLMAN,

        // GAUSS80_CHARGES, AMPAC_CHARGES, MULLIKEN_CHARGES, DICT_ CHARGES,

        // MMFF94_CHARGES, USER_CHARGES

      if (strcasecmp(dp->GetValue().c_str(),"Mulliken") == 0)

        ofs << "MULLIKEN_CHARGES" << endl;

      else if (strcasecmp(dp->GetValue().c_str(),"MMFF94") == 0)

        ofs << "MMFF94_CHARGES" << endl;

      else if (strcasecmp(dp->GetValue().c_str(),"ESP") == 0)

        ofs << "USER_CHARGES" << endl;

      else if (strcasecmp(dp->GetValue().c_str(),"Gasteiger") == 0)

        ofs << "GASTEIGER" << endl;

      else // ideally, code should pick from the Tripos types

        ofs << "USER_CHARGES" << endl;

    }

    else { // No idea what these charges are... all our code sets "PartialCharges"

        ofs << "GASTEIGER" << endl;

    }

    //    ofs << "Energy = " << mol.GetEnergy() << endl;

    if (mol.HasData(OBGenericDataType::CommentData))

      {

        ofs << "****\n"; // comment line printed, so we need to add "no status bits set"

        OBCommentData *cd = (OBCommentData*)mol.GetData(OBGenericDataType::CommentData);

        ofs << cd->GetData();

      }

    ofs << endl;

    ofs << "@<TRIPOS>ATOM" << endl;

    OBAtom *atom;

    OBResidue *res;

    vector<OBAtom*>::iterator i;

    std::map<int, int> labelcount;

    ttab.SetFromType("INT");

    ttab.SetToType("SYB");

    bool hasFormalCharges = false;

    for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))

      {

        //

        //  Use sequentially numbered atom names if no residues

        //

        snprintf(label,BUFF_SIZE, "%s%d",

                 OBElements::GetSymbol(atom->GetAtomicNum()),

                 ++labelcount[atom->GetAtomicNum()]);

        strcpy(rlabel,"<1>");

        strcpy(rnum,"1");

        str = atom->GetType();

        ttab.Translate(str1,str);

        if (atom->GetFormalCharge() != 0) hasFormalCharges = true;

        //

        //  Use original atom names if there are residues

        //

        if (!ligandsOnly && (res = atom->GetResidue()) )

          {

            // use original atom names defined by residue

            snprintf(label,BUFF_SIZE,"%s",(char*)res->GetAtomID(atom).c_str());

            // make sure that residue name includes its number

            snprintf(rlabel,BUFF_SIZE,"%s%d",res->GetName().c_str(), res->GetNum());

            snprintf(rnum,BUFF_SIZE,"%d",res->GetNum());

          }

        snprintf(buffer,BUFF_SIZE,"%7d %-6s   %9.4f %9.4f %9.4f %-5s %3s  %-8s %9.4f",

                 atom->GetIdx(),label,

                 atom->GetX(),atom->GetY(),atom->GetZ(),

                 str1.c_str(),

                 rnum,rlabel,

                 atom->GetPartialCharge());

        ofs << buffer << endl;

      }

    //store formal charge info; put before bonds so we don't have

    //to read past the end of the molecule to realize it is there

    if(hasFormalCharges && !skipFormalCharge) {

      //dkoes - to enable roundtriping of charges

      ofs << "@<TRIPOS>UNITY_ATOM_ATTR\n";

      for (atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i))

      {

        int charge = atom->GetFormalCharge();

        if (charge != 0) 

        {

          ofs << atom->GetIdx() << " 1\n"; //one attribute

          ofs << "charge " << charge << "\n"; //namely charge

        }

      }

    }

    ofs << "@<TRIPOS>BOND" << endl;

    OBBond *bond;

    vector<OBBond*>::iterator j;

    string s1, s2;

    for (bond = mol.BeginBond(j);bond;bond = mol.NextBond(j))

      {

        s1 = bond->GetBeginAtom()->GetType();

        s2 = bond->GetEndAtom()->GetType();

        if (bond->IsAromatic() || s1 == "O.co2" || s2 == "O.co2")

          strcpy(label,"ar");

        else if (bond->IsAmide())

          strcpy(label,"am");

        else

          snprintf(label,BUFF_SIZE,"%d",bond->GetBondOrder());

        snprintf(buffer, BUFF_SIZE,"%6d %5d %5d   %2s",

                 bond->GetIdx()+1,bond->GetBeginAtomIdx(),bond->GetEndAtomIdx(),

                 label);

        ofs << buffer << endl;

      }

    // NO trailing blank line (PR#1868929).

    //    ofs << endl;

    return(true);

  }

  int MOL2Format::SkipObjects(int n, OBConversion* pConv)

  {

    const char txt[] = "@<TRIPOS>MOLECULE";

    istream& ifs = *pConv->GetInStream();

    if(!ifs)

      return -1;

    if(n>0 && ifs.peek()==txt[0])

      ifs.ignore(); // move past '@' so that next mol will be found

    do {

      ignore(ifs, txt);

    } while(ifs && (--n)>0);

    if(!ifs.eof())

      ifs.seekg(1-sizeof(txt), ios::cur);//1 for '/0'

    char ch = ifs.peek();

   return 1;

  }

} //namespace OpenBabel