/* 
 *            Copyright 2009-2018 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

#include <iostream>
#include <fstream>
#include <stdexcept>

#include "votca/tools/application.h"
#include <votca/csg/trajectorywriter.h>
#include <votca/csg/trajectoryreader.h>
#include <votca/csg/topologyreader.h>
#include <votca/xtp/statesaversqlite.h>
#include <votca/xtp/version.h>
#include <votca/tools/globals.h>
#include <votca/xtp/Md2QmEngine.h>

#include <votca/xtp/topology.h>
#include <votca/xtp/molecule.h>
#include <votca/xtp/segment.h>
#include <votca/xtp/fragment.h>
#include <votca/xtp/atom.h>
#include <votca/xtp/segmenttype.h>

using namespace std;

namespace CSG = votca::csg;
namespace XTP = votca::xtp;
namespace TOOLS = votca::tools;

class XtpMap : public TOOLS::Application
{

public:
    string ProgramName()  { return "xtp_map"; }
    void   HelpText(ostream &out) {out << "Generates QM|MD topology" << endl;}
    void ShowHelpText(std::ostream &out);

    void Initialize();
    bool EvaluateOptions();
    void Run();
    void Save(string mode);

    void BeginEvaluate() { ; }
    bool DoTrajectory() { return 0; }
    bool DoMapping() { return 0; }


protected:
    TOOLS::Property               _options;
    CSG::Topology          _mdtopol;
    XTP::Topology          _qmtopol;

    XTP::Md2QmEngine            _md2qm;
    XTP::StateSaverSQLite  _statsav;
    string                 _outdb;

};

namespace propt = boost::program_options;

void XtpMap::Initialize() {

    TOOLS::Application::Initialize();
    
    CSG::TrajectoryWriter::RegisterPlugins();
    CSG::TrajectoryReader::RegisterPlugins();
    CSG::TopologyReader::RegisterPlugins();

    AddProgramOptions() ("topology,t", propt::value<string> (),
                         "  topology");
    AddProgramOptions() ("coordinates,c", propt::value<string> (),
                         "  coordinates or trajectory");
    AddProgramOptions() ("segments,s",  propt::value<string> (),
                         "  definition of segments and fragments");
    AddProgramOptions() ("file,f", propt::value<string> (),
                         "  state file");
}

bool XtpMap::EvaluateOptions() {

    CheckRequired("topology", "Missing topology file");
    CheckRequired("segments", "Missing segment definition file");
    CheckRequired("coordinates", "Missing trajectory input");
    CheckRequired("file", "Missing state file");

    return 1;
}

void XtpMap::Run() {
  
    std::string name = ProgramName();
    if (VersionString() != "") name = name + ", version " + VersionString();
    XTP::HelpTextHeader(name);

    // +++++++++++++++++++++++++++++++++++++ //
    // Initialize MD2QM Engine and SQLite Db //
    // +++++++++++++++++++++++++++++++++++++ //

    bool abort = false;
    _outdb = _op_vm["file"].as<string> ();
    _statsav.Open(_qmtopol, _outdb, false);
    int frames_in_db = _statsav.FramesInDatabase();
    if (frames_in_db > 0) {
        cout << endl << "ERROR <xtp_map> : state file '" 
             << _outdb << "' already in use. Abort." << endl;
        abort = true;
    }
    _statsav.Close();
    if (abort) return;

    string cgfile = _op_vm["segments"].as<string> ();
    _md2qm.Initialize(cgfile);

    
    // ++++++++++++++++++++++++++++ //
    // Create MD topology from file //
    // ++++++++++++++++++++++++++++ //

    // Create topology reader
    string topfile = _op_vm["topology"].as<string> ();
    CSG::TopologyReader *topread;
    topread = CSG::TopReaderFactory().Create(topfile);

    if (topread == NULL) {
        throw runtime_error( string("Input format not supported: ")
                           + _op_vm["topology"].as<string> () );
    }

    topread->ReadTopology(topfile, this->_mdtopol);
    if (TOOLS::globals::verbose) {
        cout << "Read MD topology from " << topfile << ": Found "
             << _mdtopol.BeadCount() << " atoms in "
             << _mdtopol.MoleculeCount() << " molecules. "
             << endl;
    }

    // ++++++++++++++++++++++++++++++ //
    // Create MD trajectory from file //
    // ++++++++++++++++++++++++++++++ //

    // Create trajectory reader and initialize
    string trjfile =  _op_vm["coordinates"].as<string> ();
    CSG::TrajectoryReader *trjread;
    trjread = CSG::TrjReaderFactory().Create(trjfile);

    if (trjread == NULL) {
        throw runtime_error( string("Input format not supported: ")
                           + _op_vm["coordinates"].as<string> () );
    }
    trjread->Open(trjfile);
    trjread->FirstFrame(this->_mdtopol);

    int    firstFrame = 1;
    int    nFrames    = 1;
    bool   beginAt    = 0;
    double startTime  = _mdtopol.getTime();

    if (_op_vm.count("nframes")) {
        nFrames = _op_vm["nframes"].as<int> ();
    }
    if (_op_vm.count("first-frame")) {
        firstFrame = _op_vm["first-frame"].as<int> ();
    }    
    if (_op_vm.count("begin")) {
        beginAt = true;
        startTime = _op_vm["begin"].as<double> ();
    }

    // Extract first frame specified
    bool hasFrame;

    for (hasFrame = true; hasFrame == true;
         hasFrame = trjread->NextFrame(this->_mdtopol)) {
         if (  ((_mdtopol.getTime() < startTime) && beginAt )
               || firstFrame > 1 ) {
             firstFrame--;
             continue;
         }
         break;
    }
    if ( ! hasFrame) {
        trjread->Close();
        delete trjread;

        throw runtime_error("Time or frame number exceeds trajectory length");
    }
    
    // +++++++++++++++++++++++++ //
    // Convert MD to QM Topology //
    // +++++++++++++++++++++++++ //

    for (int saved = 0; hasFrame && saved < nFrames;
         hasFrame = trjread->NextFrame(this->_mdtopol), saved++) {

        _md2qm.Md2Qm(&_mdtopol, &_qmtopol);

    // +++++++++++++++++++++++++ //
    // Save to SQLite State File //
    // +++++++++++++++++++++++++ //

        this->Save("");
    }

    // trjread->Close();
    // delete trjread;

}

void XtpMap::Save(string mode) {    
    
    _statsav.Open(_qmtopol, _outdb);

    _statsav.WriteFrame();

    _statsav.Close();

}

void XtpMap::ShowHelpText(std::ostream &out) {
    string name = ProgramName();
    if (VersionString() != "") name = name + ", version " + VersionString();
    XTP::HelpTextHeader(name);
    HelpText(out);
    //out << "\n\n" << OptionsDesc() << endl;
    out << "\n\n" << VisibleOptions() << endl;
}


int main(int argc, char** argv)
{
    XtpMap xtpmap;
    return xtpmap.Exec(argc, argv);
}