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

Copyright (C) 2011 by Geoffrey Hutchison

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 <cstdlib>

using namespace std;

namespace OpenBabel

{

#define BOHR_TO_ANGSTROM 0.5291772108

#define ANGSTROM_TO_BOHR 1.0 / BOHR_TO_ANGSTROM

  class ABINITFormat : public OBMoleculeFormat

  {

  public:

    //Register this format type ID

    ABINITFormat()

    {

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

    }

    const char* Description() override  // required

    {

      return

        "ABINIT Output Format\n"

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

        "  s  Output single bonds only\n"

        "  b  Disable bonding entirely\n\n";

    }

    const char* SpecificationURL() override

    { return "http://abinit.org/" ; }  // optional

    //Flags() can return be any the following combined by | or be omitted if none apply

    // NOTREADABLE  READONEONLY  NOTWRITABLE  WRITEONEONLY

    unsigned int Flags() override

    {

      return READONEONLY | NOTWRITABLE;

    }

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

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

    /// The "API" interface functions

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

  };

  //***

  //Make an instance of the format class

  ABINITFormat theABINITFormat;

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

  bool ABINITFormat::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;

    const char* title = pConv->GetTitle();

    char buffer[BUFF_SIZE];

    string str;

    vector<string> vs;

    OBAtom *atom;

    int natom;

    vector<int> atomicNumbers, atomTypes;

    double x, y, z;

    vector<vector3> atomPositions;

    vector<double> energies;

    // Translation vectors (if present)

    // aka rprim

    vector3 translationVectors[3];

    double acell[3]; // scale of lattice vectors

    int numTranslationVectors = 0;

    int symmetryCode = 0;

    bool last = false;

    mol.BeginModify();

    while (ifs.getline(buffer,BUFF_SIZE))

      {

  if(strstr(buffer,"outvars")){

       last = true;

       continue;

     }

  if(! last){

    continue;

  }

        // tokens are listed in alphabetical order for code clarity

        if (strstr(buffer, "acell")) {

          tokenize(vs, buffer);

          if (vs.size() < 4)

            continue; // invalid line

          // acell=  7.6967369631E+00  7.6967369631E+00  7.6967369631E+00

          for (int i = 0; i < 3; ++i) {

            acell[i] = atof(vs[i+1].c_str());

          }

        }

        else if (strstr(buffer, " xcart ")) {

          double unit = BOHR_TO_ANGSTROM;

          if (strstr(buffer, "ngstrom"))

            unit = 1.0; // no conversion needed

    //          ifs.getline(buffer,BUFF_SIZE);

          tokenize(vs, buffer);

          // first line, rprim takes up a token

          x = atof((char*)vs[1].c_str()) * unit;

          y = atof((char*)vs[2].c_str()) * unit;

    z = atof((char*)vs[3].c_str()) * unit;

    atomPositions.push_back(vector3(x, y, z));

    // get next line

    ifs.getline(buffer,BUFF_SIZE);

    tokenize(vs, buffer);

    //rest of lines

    while(vs.size() == 3) {

      x = atof(vs[0].c_str()) * unit;

      y = atof(vs[1].c_str()) * unit;

      z = atof(vs[2].c_str()) * unit;

            atomPositions.push_back(vector3(x, y, z));

            // get next line

            ifs.getline(buffer,BUFF_SIZE);

            tokenize(vs, buffer);

          }

      }

        else if (strstr(buffer, "natom")) {

          tokenize(vs, buffer);

          if (vs.size() != 2)

            continue;

          natom = atoi(vs[1].c_str());

        }

        else if (strstr(buffer, "rprim")) {

          numTranslationVectors = 0;

          int column;

    ifs.getline(buffer,BUFF_SIZE);

          for (int i = 1; i < 4; ++i) {

            tokenize(vs, buffer);

            if (vs.size() < 3)

              break;

            // first line, rprim takes up a token

      //            if (i == 0)

      // column = 1;

            //else

              column = 0;

            x = atof((char*)vs[column].c_str()) * BOHR_TO_ANGSTROM;

            y = atof((char*)vs[column+1].c_str()) * BOHR_TO_ANGSTROM;

            z = atof((char*)vs[column+2].c_str()) * BOHR_TO_ANGSTROM;

            translationVectors[numTranslationVectors++].Set(x, y,z);

            ifs.getline(buffer,BUFF_SIZE);

          }

        }

        else if (strstr(buffer, "Symmetries")) {

          tokenize(vs, buffer, "()");

          // Should be something like (#160)

          symmetryCode = atoi(vs[1].substr(1).c_str());

        }

        else if (strstr(buffer, "typat")) {

          atomTypes.clear();

    int n=0;

      while( n<=natom){

        tokenize(vs, buffer);

        for (unsigned int i = 1; i < vs.size(); ++i) {

    atomTypes.push_back(atoi(vs[i].c_str()));

        }

        n+=vs.size();

        ifs.getline(buffer,BUFF_SIZE);

      }

        }

        else if (strstr(buffer, "znucl")) {

          tokenize(vs, buffer);

          // make sure znucl is first token

          if (vs[0] != "znucl")

            continue;

          // push back the remaining tokens into atomicNumbers

          atomicNumbers.clear();

          atomicNumbers.push_back(0); // abinit starts typat with 1

          for (unsigned int i = 1; i < vs.size(); ++i)

            atomicNumbers.push_back(int(atof(vs[i].c_str())));

        }

        // xangst

        // forces

      }

    for (int i = 0; i < natom; ++i) {

      atom = mol.NewAtom();

      //set atomic number

      int idx = atom->GetIdx();

      int type = atomTypes[idx - 1];

      atom->SetAtomicNum(atomicNumbers[type]);

      // we set the coordinates by conformers in another loop

    }

    mol.EndModify();

    if (natom == 0)

      return false;

    int numConformers = atomPositions.size() / natom;

    for (int i = 0; i < numConformers; ++i) {

      double *coordinates = new double[natom * 3];

      for (int j = 0; j < natom; ++j) {

        vector3 currentPosition = atomPositions[i*natom + j];

        coordinates[j*3] = currentPosition.x();

        coordinates[j*3 + 1] = currentPosition.y();

        coordinates[j*3 + 2] = currentPosition.z();

      }

      mol.AddConformer(coordinates);

    }

    // Delete first conformer, created by EndModify, bunch of 0s

    mol.DeleteConformer(0);

    // Set geometry to last one

    mol.SetConformer(mol.NumConformers() - 1);

    if (!pConv->IsOption("b",OBConversion::INOPTIONS))

      mol.ConnectTheDots();

    if (!pConv->IsOption("s",OBConversion::INOPTIONS) && !pConv->IsOption("b",OBConversion::INOPTIONS))

      mol.PerceiveBondOrders();

    // Attach unit cell translation vectors if found

    if (numTranslationVectors > 0) {

      OBUnitCell* uc = new OBUnitCell;

      //      uc->SetData(acell[0] * translationVectors[0], acell[1] * translationVectors[1], acell[2] * translationVectors[2]);

      uc->SetData(translationVectors[0], translationVectors[1], translationVectors[2]);

      uc->SetOrigin(fileformatInput);

      if (symmetryCode)

        uc->SetSpaceGroup(symmetryCode);

      mol.SetData(uc);

    }

    mol.SetTitle(title);

    return(true);

  }

} //namespace OpenBabel