<commit_before>/*
 *            Copyright 2009-2017 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 "orca.h"
#include <votca/xtp/qminterface.h>

#include <votca/ctp/segment.h>

#include <votca/xtp/elements.h>
#include <boost/algorithm/string.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/format.hpp>
#include <boost/filesystem.hpp>
#include <stdio.h>
#include <iomanip>
#include <sys/stat.h>
#include <vector>



namespace votca { namespace xtp {
    namespace ub = boost::numeric::ublas;
    

void Orca::Initialize( Property *options ) {
    
    //good luck

     // Orca file names
    std::string fileName = "system";

    _xyz_file_name = fileName + ".xyz";
    _input_file_name = fileName + ".inp";
    _log_file_name = fileName + ".log"; 
    _shell_file_name = fileName + ".sh"; 
    _orb_file_name=fileName+".gbw";

    std::string key = "package";
    std::string _name = options->get(key+".name").as<std::string> ();
    
    if ( _name != "orca" ) {
        cerr << "Tried to use " << _name << " package. ";
        throw std::runtime_error( "Wrong options file");
    }
    
    _executable =       options->get(key + ".executable").as<std::string> ();
    _charge =           options->get(key + ".charge").as<int> ();
    _spin =             options->get(key + ".spin").as<int> ();
    _options =          options->get(key + ".options").as<std::string> ();
    _memory =           options->get(key + ".memory").as<std::string> ();
    _threads =          options->get(key + ".threads").as<int> ();
    _scratch_dir =      options->get(key + ".scratch").as<std::string> ();
    _basisset_name =    options->get(key + ".basisset").as<std::string> ();
    _cleanup =          options->get(key + ".cleanup").as<std::string> ();

    
    
    if (options->exists(key + ".outputVxc")) {
                _output_Vxc = options->get(key + ".outputVxc").as<bool> ();   
            }
             else _output_Vxc=false;
    if (_output_Vxc){
        throw std::runtime_error( "Sorry "+_name+" does not support Vxc output");
    }
  

    _write_basis_set=   options->get(key + ".writebasisset").as<bool> ();
    _write_pseudopotentials= options->get(key + ".writepseudopotentials").as<bool> ();
    

    // check if the optimize keyword is present, if yes, read updated coords
    std::string::size_type iop_pos = _options.find(" Opt"); /*optimization word in orca*/
    if (iop_pos != std::string::npos) {
        _is_optimization = true;
    } else
    {
        _is_optimization = false;
    }

    // check if the esp keyword is present, if yes, get the charges and save them
    iop_pos = _options.find(" chelpg");  /*electrostatic potential related to partial atomic charges I guess is chelpg in orca but check */
    if (iop_pos != std::string::npos) {
        _get_charges = true;
    } else
    {
        _get_charges = false;
    }

    // check if the charge keyword is present, if yes, get the self energy and save it
    iop_pos = _options.find("pointcharges");  /*??*/
    if (iop_pos != std::string::npos) {
        _get_self_energy = true;
        _write_charges = true;
    } else
    {
        _get_self_energy = false;
        _write_charges = false;
    }
    
    // check if the guess should be prepared, if yes, append the guess later
    _write_guess = false;
    iop_pos = _options.find("iterations 1 ");
    if (iop_pos != std::string::npos) _write_guess = true;
    iop_pos = _options.find("iterations 1\n");
    if (iop_pos != std::string::npos) _write_guess = true;
}    

/**
 * Prepares the *.inp file from a vector of segments
 * Appends a guess constructed from monomer orbitals if supplied
 */
bool Orca::WriteInputFile( std::vector<ctp::Segment* > segments, Orbitals* orbitals_guess ){
   
    std::vector<std::string> results;
   // int qmatoms = 0;
    std::string temp_suffix = "/id";
    std::string scratch_dir_backup = _scratch_dir;
    std::ofstream _com_file;
    std::ofstream _crg_file;

    std::string _com_file_name_full = _run_dir + "/" + _input_file_name;
    std::string _crg_file_name_full = _run_dir + "/background.crg" ;
    

    _com_file.open ( _com_file_name_full.c_str() );
    // header 
   //_com_file << "* xyzfile  "  <<  _charge << " " << _spin << " " << _xyz_file_name << "\n" << endl;
    _com_file << "* xyz  "  <<  _charge << " " << _spin << endl;
    
    std::vector< ctp::QMAtom* > qmatoms;
    

       // This is needed for the QM/MM scheme, since only orbitals have 
            // updated positions of the QM region, hence vector<Segments*> is 
            // NULL in the QMMachine and the QM region is also printed here
    if(_write_charges){
        qmatoms= orbitals_guess->QMAtoms();
    }
    else{
        QMMInterface qmmface;
        qmatoms=qmmface.Convert(segments);
    }
    
   
    std::vector< ctp::QMAtom* >::iterator it;

         
    for (it = qmatoms.begin(); it < qmatoms.end(); it++) {
        if ((*it)->from_environment) {continue;}
        
        _com_file << setw(3) << (*it)->type.c_str()
                  << setw(12) << setiosflags(ios::fixed) << setprecision(5) << (*it)->x
                  << setw(12) << setiosflags(ios::fixed) << setprecision(5) << (*it)->y
                  << setw(12) << setiosflags(ios::fixed) << setprecision(5) << (*it)->z
                  << endl;
    }
    _com_file << "* \n" << endl;
    _com_file << "%pal\n "  <<  "nprocs " <<  _threads  << "\nend" << "\n" << endl;
       
    // basis set info
    if ( _write_basis_set){
        Elements _elements;

        list<std::string> elements;
        BasisSet bs;
        bs.LoadBasisSet(_basisset_name);
        CTP_LOG(ctp::logDEBUG, *_pLog) << "Loaded Basis Set " << _basisset_name << flush;
        ofstream _el_file;
        std::string _el_file_name = _run_dir + "/" +  "system.bas";
        _el_file.open(_el_file_name.c_str());
        //_com_file << "@" << "system.bas" << endl;
        _el_file << "$DATA" << endl;

        for (it = qmatoms.begin(); it < qmatoms.end(); it++) {
            if ((*it)->from_environment) { continue;}
            
            std::string element_name = (*it)->type;
            list<std::string>::iterator ite;
            ite = find(elements.begin(), elements.end(), element_name);
            if (ite == elements.end()) {
                elements.push_back(element_name);
                Element* element = bs.getElement(element_name);
                _el_file << _elements.getEleFull(element_name) << endl; 
                for (Element::ShellIterator its = element->firstShell(); its != element->lastShell(); its++) {
                    Shell* shell = (*its);

                    string type=shell->getType();
                    // check combined shells

                       for(unsigned i=0;i<type.size();++i){
                           string subtype = string( type, i, 1 );
                          _el_file  << subtype << " " << shell->getSize() << endl; 
                          int _sh_idx =0;
                          for (Shell::GaussianIterator itg = shell->firstGaussian(); itg != shell->lastGaussian(); itg++) {
                              GaussianPrimitive* gaussian = *itg;
                              _sh_idx++;
                              _el_file << " " << _sh_idx << " " << indent(gaussian->decay);
                              _el_file << " " << indent(gaussian->contraction[FindLmax(subtype)]);

                              }
                              _el_file << endl;
                          }
                    }

                 }
             }
        
        _el_file << "STOP\n";
        _el_file.close();


        _com_file << "%basis\n "  << endl;
        _com_file << "GTOName" << " " << "=" << "\"system.bas\";"  << endl;
    } // write_basis set
            
    // ECPs
   /* WRITING ECP INTO system.inp FILE for ORCA**/
    if ( _write_pseudopotentials ){
        std::string pseudopotential_name("ecp");
        _com_file << endl;

        list<std::string> elements;
        elements.push_back("H");
        elements.push_back("He");
        BasisSet ecp;
        ecp.LoadPseudopotentialSet(pseudopotential_name);
        CTP_LOG(ctp::logDEBUG, *_pLog) << "Loaded Pseudopotentials " << pseudopotential_name << flush;



         for (it = qmatoms.begin(); it < qmatoms.end(); it++) {
            if (!(*it)->from_environment) {
                std::string element_name = (*it)->type;


                list<std::string>::iterator ite;
                ite = find(elements.begin(), elements.end(), element_name);
                if (ite == elements.end()) {
                    elements.push_back(element_name);
                    Element* element = ecp.getElement(element_name);
                    _com_file << "\n" << "NewECP" << " " <<  element_name  << endl;
                    _com_file << "N_core" << " " << element->getNcore()  << endl;
                    //lmaxnum2lmaxname 
                    _com_file << "lmax" << " " << getLName(element->getLmax()) << endl;

                    //For Orca the order doesn't matter but let's write it in ascending order
                    // write remaining shells in ascending order s,p,d...
                    for (int i = 0; i < element->getLmax(); i++) {
                        for (Element::ShellIterator its = element->firstShell(); its != element->lastShell(); its++) {
                            Shell* shell = (*its);
                            if (shell->getLmax() == i) {
                            // shell type, number primitives, scale factor
                            _com_file << shell->getType() << " " <<  shell->getSize() << endl;
                            int _sh_idx = 0;
                            for (Shell::GaussianIterator itg = shell->firstGaussian(); itg != shell->lastGaussian(); itg++) {
                                GaussianPrimitive* gaussian = *itg;
                                _sh_idx++;
                                _com_file << _sh_idx << " " << gaussian->decay << " " << gaussian->contraction[0] << " " << gaussian->power << endl;
                            }
                        }
                    }
                }

                for (Element::ShellIterator its = element->firstShell(); its != element->lastShell(); its++) {
                    Shell* shell = (*its);
                    // shell type, number primitives, scale factor
                    if (shell->getLmax() == element->getLmax()) {
                        _com_file << shell->getType() << " " <<  shell->getSize() << endl;
                        // _com_file << shell->getSize() << endl;
                        int _sh_idx = 0;
                        for (Shell::GaussianIterator itg = shell->firstGaussian(); itg != shell->lastGaussian(); itg++) {
                            GaussianPrimitive* gaussian = *itg;
                            _sh_idx++;
                            _com_file << _sh_idx << " " << gaussian->decay << " " << gaussian->contraction[0] << " " << gaussian->power << endl;
                            }
                        }
                    }

                   _com_file << "end\n "  <<  "\n" << endl; 
                }

            }
        }
    } // write pseudopotentials
                
                
    /* END   OF WRITING ECP INTO system.inp FILE for ORCA*************/                
    _com_file << "end\n "  <<  "\n" << endl;   //This end is for the basis set block
             
            
     if(_write_charges){      
        // actual writing of charges    
        _crg_file.open( _crg_file_name_full.c_str() ); 
        int _total_background = 0;
        for (it = qmatoms.begin(); it < qmatoms.end(); it++) {
            if ((*it)->from_environment) {
                if ((*it)->charge != 0.0) _total_background++;
            }
        }

        _crg_file << _total_background << endl;
        for (it = qmatoms.begin(); it < qmatoms.end(); it++) {
            if ((*it)->from_environment) {
                boost::format fmt("%1$+1.7f %2$+1.7f %3$+1.7f %4$+1.7f");
                fmt % (*it)->charge % (*it)->x % (*it)->y % (*it)->z ;
                if ((*it)->charge != 0.0) _crg_file << fmt << endl;
            }
        }
     }            
     
    // writing scratch_dir info
    if ( _scratch_dir != "" ) {
        
        CTP_LOG(ctp::logDEBUG,*_pLog) << "Setting the scratch dir to " << _scratch_dir + temp_suffix << flush;

        // boost::filesystem::create_directories( _scratch_dir + temp_suffix );
        std::string _temp( "scratch_dir " + _scratch_dir + temp_suffix + "\n" );
        //_com_file << _temp ;
    }

    _com_file << _options << "\n";

    
    if ( _write_guess ) { 
       
            throw std::runtime_error( "Not implemented in orca");
    }   

    _com_file << endl;
    _com_file.close();
    
      // and now generate a shell script to run both jobs, if neccessary
    CTP_LOG(ctp::logDEBUG, *_pLog) << "Setting the scratch dir to " << _scratch_dir + temp_suffix << flush;

    _scratch_dir = scratch_dir_backup + temp_suffix;
            
            //boost::filesystem::create_directories(_scratch_dir + temp_suffix);
            //std::string _temp("scratch_dir " + _scratch_dir + temp_suffix + "\n");
            //_com_file << _temp;
    WriteShellScript();
    _scratch_dir = scratch_dir_backup;
   
    return true;
}



bool Orca::WriteShellScript() {
    ofstream _shell_file;
    
    std::string _shell_file_name_full = _run_dir + "/" + _shell_file_name;
            
    _shell_file.open ( _shell_file_name_full.c_str() );

    _shell_file << "#!/bin/bash" << endl ;
    _shell_file << "mkdir -p " << _scratch_dir << endl;
    
    if ( _threads == 1 ){ 
        _shell_file << _executable << " " << _input_file_name << " > " <<  _log_file_name  << endl;//" 2> run.error" << endl;    
    } else {
        _shell_file << _executable << " "<< _input_file_name << " > " <<  _log_file_name <<endl;// " 2> run.error" << endl;    
    }
    _shell_file.close();
    
    return true;   
}

/**
 * Runs the Orca job. 
 */
bool Orca::Run()
{

    CTP_LOG(ctp::logDEBUG,*_pLog) << "Running Orca job" << flush;
    
    if (std::system(NULL)) {
        
        // Orca overrides input information, if *.db and *.movecs files are present
        // better trash the old version
        std::string file_name = _run_dir + "/system.db";
        remove ( file_name.c_str() );
        file_name = _run_dir + "/" + _log_file_name;
        remove ( file_name.c_str() );
        file_name = _run_dir + "/" + _orb_file_name;
        //remove ( file_name.c_str() );
               
        std::string _command;
        if ( _threads == 1 ) {
            _command = "cd " + _run_dir + "; sh " + _shell_file_name;
        } else {
            _command = "cd " + _run_dir + "; sh " + _shell_file_name;
        }
        //CTP_LOG(logDEBUG,*_pLog) << _command << flush;
       int check=std::system(_command.c_str());
        if (check==-1){
            CTP_LOG(ctp::logERROR, *_pLog) << _input_file_name << " failed to start" << flush;
            return false;
        }
       
        if ( CheckLogFile() ) {
            CTP_LOG(ctp::logDEBUG,*_pLog) << "Finished Orca job" << flush;
            return true;
        } else {
            CTP_LOG(ctp::logDEBUG,*_pLog) << "Orca job failed" << flush;
        }
    }
    else {
        CTP_LOG(ctp::logERROR,*_pLog) << _input_file_name << " failed to start" << flush; 
        return false;
    }
    
    return true;
}

/**
 * Cleans up after the Orca job
 */
void Orca::CleanUp() {
    
    // cleaning up the generated files
    if ( _cleanup.size() != 0 ) {
        Tokenizer tok_cleanup(_cleanup, ",");
        std::vector <std::string> _cleanup_info;
        tok_cleanup.ToVector(_cleanup_info);
        
        std::vector<std::string> ::iterator it;
               
        for (it = _cleanup_info.begin(); it != _cleanup_info.end(); ++it) {
            if ( *it == "inp" ) {
                std::string file_name = _run_dir + "/" + _input_file_name;
                remove ( file_name.c_str() );
            }
            
            if ( *it == "bas" ) {
                std::string file_name = _run_dir +"/system.bas";
                remove ( file_name.c_str() );
            }
            
            if ( *it == "log" ) {
                std::string file_name = _run_dir + "/" + _log_file_name;
                remove ( file_name.c_str() );
            }

           if ( *it == "gbw" ) {
                std::string file_name = _run_dir + "/" + _orb_file_name;
                remove ( file_name.c_str() );
            }
            
            if ( *it == "ges" ) {
                std::string file_name = _run_dir + "/system.ges" ;
                remove ( file_name.c_str() );
            }    
            if ( *it == "prop" ) {
                std::string file_name = _run_dir + "/system.prop" ;
                remove ( file_name.c_str() );
            }     
        }
    }
return; 
}


bool Orca::ParseLogFile( Orbitals* _orbitals )
{
    const double _conv_Hrt_eV = tools::conv::hrt2ev;
    
    _orbitals->setQMpackage("orca");
    _orbitals->setDFTbasis( _basisset_name);
    CTP_LOG(ctp::logDEBUG,*_pLog) << "Parsing " << _log_file_name << flush;
      // return true;
  std::string _log_file_name_full =  _run_dir + "/" + _log_file_name;
   // check if LOG file is complete
   if ( !CheckLogFile() ) return false;
    
    //std::map <int, std::std::vector<double> > _coefficients;
    std::map <int, double> _energies;
    std::map <int, double> _occ;
    
    std::string _line;
    unsigned _levels = 0;
    //unsigned _level;
    //unsigned _basis_size = 0;
    int _number_of_electrons = 0;
    //bool _has_basis_dim = false;
               
    std::vector<std::string> results;    
    
    std::ifstream _input_file( _log_file_name_full.c_str() );
    
    if (_input_file.fail()) {
        CTP_LOG( ctp::logERROR, *_pLog ) << "File " << _log_file_name_full << " not found " << flush;
        return false;
    } else {
        CTP_LOG(ctp::logDEBUG, *_pLog) << "Reading Coordinates and occupationnumbers and energies from " << _log_file_name_full << flush;
    }
    
      
          //Coordinates of the final configuration depending on whether it is an optimization or not
         
        
        
    while (_input_file) {
        getline(_input_file,_line);
        boost::trim(_line);
        
        
        
        if ( _is_optimization ){
              throw runtime_error("Not implemented yet!");
        }
        bool _found_optimization=true;
        
        std::string::size_type coordinates_pos = _line.find("CARTESIAN COORDINATES (ANGSTROEM)");
        
        if ( _found_optimization && coordinates_pos != std::string::npos) {
            CTP_LOG(ctp::logDEBUG,*_pLog) << "Getting the coordinates" << flush;
            
            //_has_coordinates = true;
            bool _has_QMAtoms = _orbitals->hasQMAtoms();

            // three garbage lines
            getline(_input_file, _line);
            // now starts the data in format
            // _id type Qnuc x y z 
            std::vector<std::string> _row;
            getline(_input_file, _line);
            boost::trim( _line );

            boost::algorithm::split( _row , _line, boost::is_any_of("\t "), boost::algorithm::token_compress_on); 
            int nfields =  _row.size();

            int atom_id = 0;   
            while ( nfields == 4 ) {
                //int atom_id = boost::lexical_cast< int >( _row.at(0) );
                //int atom_number = boost::lexical_cast< int >( _row.at(0) );
                std::string _atom_type = _row.at(0);
                double _x =  boost::lexical_cast<double>( _row.at(1) );
                double _y =  boost::lexical_cast<double>( _row.at(2) );
                double _z =  boost::lexical_cast<double>( _row.at(3) );
                //if ( tools::globals::verbose ) cout << "... ... " << atom_id << " " << atom_type << " " << atom_charge << endl;
                getline(_input_file, _line);
                boost::trim( _line );
                boost::algorithm::split( _row , _line, boost::is_any_of("\t "), boost::algorithm::token_compress_on);  
                nfields =  _row.size();

                if ( _has_QMAtoms == false ) {
                    _orbitals->AddAtom( _atom_type, _x, _y, _z );
                } else {
                                       
                    ctp::QMAtom* pAtom = _orbitals->_atoms.at( atom_id  );
                    pAtom->type = _atom_type;
                    pAtom->x = _x;
                    pAtom->y = _y;
                    pAtom->z = _z;
                    atom_id++;
                }
                    
            }
      
        }
              
        std::string::size_type energy_pos = _line.find("Total Energy");
       if (energy_pos != std::string::npos) {
            //cout << _line << endl;
            boost::algorithm::split(results, _line, boost::is_any_of(" "), boost::algorithm::token_compress_on);
            std::string _energy = results[3];
            boost::trim( _energy );
            //cout << _energy << endl; 
            _orbitals->setQMEnergy ( _conv_Hrt_eV * boost::lexical_cast<double>(_energy) );
            CTP_LOG(ctp::logDEBUG, *_pLog) << (boost::format("QM energy[eV]: %4.6f ") % _orbitals->getQMEnergy()).str() << flush;
            // _orbitals->_has_qm_energy = true;
            }
        
        /* Check for ScaHFX = factor of HF exchange included in functional */
              std::string::size_type HFX_pos = _line.find("Fraction HF Exchange ScalHFX");
              if (HFX_pos != std::string::npos) {
                  boost::algorithm::split(results, _line, boost::is_any_of(" "), boost::algorithm::token_compress_on);
                  double _ScaHFX = boost::lexical_cast<double>(results.back());
                  _orbitals->setScaHFX(_ScaHFX);
                  CTP_LOG(ctp::logDEBUG, *_pLog) << "DFT with " << _ScaHFX << " of HF exchange!" << flush;
              }
     
    //Finding Basis Dimension, the number of energy levels
        std::string::size_type dim_pos = _line.find("Basis Dimension");
        if (dim_pos != std::string::npos) {
                
            boost::algorithm::split(results, _line, boost::is_any_of(" "), boost::algorithm::token_compress_on);
            //_has_basis_dim = true;
            std::string _dim = results[4];  //The 4th element of results vector is the Basis Dim
            boost::trim( _dim );
            _levels = boost::lexical_cast<int>(_dim);
            //cout <<  boost::lexical_cast<int>(_dim) << endl;
            //_basis_size = _levels;
              CTP_LOG(ctp::logDEBUG,*_pLog) << "Basis Dimension: " << _levels << flush;
              CTP_LOG( ctp::logDEBUG, *_pLog ) << "Energy levels: " << _levels << flush;
        }
        /********************************************************/
         
        
        std::string::size_type OE_pos = _line.find("ORBITAL ENERGIES");
        if (OE_pos != std::string::npos) {
        getline(_input_file,_line);
        getline(_input_file,_line);
        getline(_input_file,_line);
        if (_line.find("E(Eh)")==std::string::npos){
            CTP_LOG(ctp::logDEBUG,*_pLog) << "Warning: Orbital Energies not found in log file" << flush;
        }
        for (unsigned i=0; i<_levels ; i++){
            getline(_input_file,_line);
            boost::trim( _line );
            boost::algorithm::split(results, _line, boost::is_any_of(" "), boost::algorithm::token_compress_on);
           
                    
            std::string _no =results[0];
            
            boost::trim( _no );
            unsigned levelnumber= boost::lexical_cast<unsigned>(_no);
            if (levelnumber!=i){
                CTP_LOG(ctp::logDEBUG,*_pLog) << "Have a look at the orbital energies something weird is going on" << flush;
            }
            std::string _oc = results[1];
            boost::trim( _oc );
            double occ =boost::lexical_cast<double>(_oc);
            // We only count alpha electrons, each orbital must be empty or doubly occupied
            if (occ==2){
                _number_of_electrons++;
                _occ[i]=occ;
            }
            else if (occ==0) {
                _occ[i]=occ;
            }
            else {
                throw runtime_error("Only empty or doubly occupied orbitals are allowed not running the right kind of DFT calculation");
            }
            
            std::string _e = results[2];
            boost::trim( _e );
            _energies [i] = boost::lexical_cast<double>(_e) ;
        }
    }
            /*
                 *  Partial charges from the input file
                 */
        std::string::size_type charge_pos = _line.find("CHELPG Charges");

        if (charge_pos != std::string::npos && _get_charges) {
                CTP_LOG(ctp::logDEBUG, *_pLog) << "Getting charges" << flush;
                //_has_charges = true;
                getline(_input_file, _line);
                //getline(_input_file, _line);

                bool _has_atoms = _orbitals->hasQMAtoms();

                std::vector<std::string> _row;
                getline(_input_file, _line);
                boost::trim(_line);
                //cout << _line << endl;
                boost::algorithm::split(_row, _line, boost::is_any_of("\t "), boost::algorithm::token_compress_on);
                int nfields = _row.size();
                //cout << _row.size() << endl;

                while (nfields == 4) {
                    int atom_id = boost::lexical_cast< int >(_row.at(0));
                    atom_id++;
                    //int atom_number = boost::lexical_cast< int >(_row.at(0));
                    std::string atom_type = _row.at(1);
                    double atom_charge = boost::lexical_cast< double >(_row.at(3));
                    if ( tools::globals::verbose ) cout << "... ... " << atom_id << " " << atom_type << " " << atom_charge << endl;
                    getline(_input_file, _line);
                    boost::trim(_line);
                    boost::algorithm::split(_row, _line, boost::is_any_of("\t "), boost::algorithm::token_compress_on);
                    nfields = _row.size();

                    if (_has_atoms == false) {
                        _orbitals->AddAtom(atom_type, 0, 0, 0, atom_charge);
                    } if (charge_pos != std::string::npos && _get_charges) {
                CTP_LOG(ctp::logDEBUG, *_pLog) << "Getting charges" << flush;
                //_has_charges = true;
                getline(_input_file, _line);
                //getline(_input_file, _line);

                bool _has_atoms = _orbitals->hasQMAtoms();

                std::vector<std::string> _row;
                getline(_input_file, _line);
                boost::trim(_line);
                //cout << _line << endl;
                boost::algorithm::split(_row, _line, boost::is_any_of("\t "), boost::algorithm::token_compress_on);
                int nfields = _row.size();
                //cout << _row.size() << endl;

                while (nfields == 4) {
                    int atom_id = boost::lexical_cast< int >(_row.at(0));
                    atom_id++;
                    //int atom_number = boost::lexical_cast< int >(_row.at(0));
                    std::string atom_type = _row.at(1);
                    double atom_charge = boost::lexical_cast< double >(_row.at(3));
                    if ( tools::globals::verbose ) cout << "... ... " << atom_id << " " << atom_type << " " << atom_charge << endl;
                    getline(_input_file, _line);
                    boost::trim(_line);
                    boost::algorithm::split(_row, _line, boost::is_any_of("\t "), boost::algorithm::token_compress_on);
                    nfields = _row.size();

                    if (_has_atoms == false) {
                        _orbitals->AddAtom(atom_type, 0, 0, 0, atom_charge);
                    } else {
                        ctp::QMAtom* pAtom = _orbitals->_atoms.at(atom_id - 1);
                        pAtom->type = atom_type;
                        pAtom->charge = atom_charge;
                    }

                }
                //_orbitals->_has_atoms = true;
            }else {
                        ctp::QMAtom* pAtom = _orbitals->_atoms.at(atom_id - 1);
                        pAtom->type = atom_type;
                        pAtom->charge = atom_charge;
                    }<commit_msg>fixed double reading of CHELPG charges from ORCA logfile, probably bogus merge issue<commit_after>/*
 *            Copyright 2009-2017 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 "orca.h"
#include <votca/xtp/qminterface.h>

#include <votca/ctp/segment.h>

#include <votca/xtp/elements.h>
#include <boost/algorithm/string.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/format.hpp>
#include <boost/filesystem.hpp>
#include <stdio.h>
#include <iomanip>
#include <sys/stat.h>
#include <vector>



namespace votca { namespace xtp {
    namespace ub = boost::numeric::ublas;
    

void Orca::Initialize( Property *options ) {
    
    //good luck

     // Orca file names
    std::string fileName = "system";

    _xyz_file_name = fileName + ".xyz";
    _input_file_name = fileName + ".inp";
    _log_file_name = fileName + ".log"; 
    _shell_file_name = fileName + ".sh"; 
    _orb_file_name=fileName+".gbw";

    std::string key = "package";
    std::string _name = options->get(key+".name").as<std::string> ();
    
    if ( _name != "orca" ) {
        cerr << "Tried to use " << _name << " package. ";
        throw std::runtime_error( "Wrong options file");
    }
    
    _executable =       options->get(key + ".executable").as<std::string> ();
    _charge =           options->get(key + ".charge").as<int> ();
    _spin =             options->get(key + ".spin").as<int> ();
    _options =          options->get(key + ".options").as<std::string> ();
    _memory =           options->get(key + ".memory").as<std::string> ();
    _threads =          options->get(key + ".threads").as<int> ();
    _scratch_dir =      options->get(key + ".scratch").as<std::string> ();
    _basisset_name =    options->get(key + ".basisset").as<std::string> ();
    _cleanup =          options->get(key + ".cleanup").as<std::string> ();

    
    
    if (options->exists(key + ".outputVxc")) {
                _output_Vxc = options->get(key + ".outputVxc").as<bool> ();   
            }
             else _output_Vxc=false;
    if (_output_Vxc){
        throw std::runtime_error( "Sorry "+_name+" does not support Vxc output");
    }
  

    _write_basis_set=   options->get(key + ".writebasisset").as<bool> ();
    _write_pseudopotentials= options->get(key + ".writepseudopotentials").as<bool> ();
    

    // check if the optimize keyword is present, if yes, read updated coords
    std::string::size_type iop_pos = _options.find(" Opt"); /*optimization word in orca*/
    if (iop_pos != std::string::npos) {
        _is_optimization = true;
    } else
    {
        _is_optimization = false;
    }

    // check if the esp keyword is present, if yes, get the charges and save them
    iop_pos = _options.find(" chelpg");  /*electrostatic potential related to partial atomic charges I guess is chelpg in orca but check */
    if (iop_pos != std::string::npos) {
        _get_charges = true;
    } else
    {
        _get_charges = false;
    }

    // check if the charge keyword is present, if yes, get the self energy and save it
    iop_pos = _options.find("pointcharges");  /*??*/
    if (iop_pos != std::string::npos) {
        _get_self_energy = true;
        _write_charges = true;
    } else
    {
        _get_self_energy = false;
        _write_charges = false;
    }
    
    // check if the guess should be prepared, if yes, append the guess later
    _write_guess = false;
    iop_pos = _options.find("iterations 1 ");
    if (iop_pos != std::string::npos) _write_guess = true;
    iop_pos = _options.find("iterations 1\n");
    if (iop_pos != std::string::npos) _write_guess = true;
}    

/**
 * Prepares the *.inp file from a vector of segments
 * Appends a guess constructed from monomer orbitals if supplied
 */
bool Orca::WriteInputFile( std::vector<ctp::Segment* > segments, Orbitals* orbitals_guess ){
   
    std::vector<std::string> results;
   // int qmatoms = 0;
    std::string temp_suffix = "/id";
    std::string scratch_dir_backup = _scratch_dir;
    std::ofstream _com_file;
    std::ofstream _crg_file;

    std::string _com_file_name_full = _run_dir + "/" + _input_file_name;
    std::string _crg_file_name_full = _run_dir + "/background.crg" ;
    

    _com_file.open ( _com_file_name_full.c_str() );
    // header 
   //_com_file << "* xyzfile  "  <<  _charge << " " << _spin << " " << _xyz_file_name << "\n" << endl;
    _com_file << "* xyz  "  <<  _charge << " " << _spin << endl;
    
    std::vector< ctp::QMAtom* > qmatoms;
    

       // This is needed for the QM/MM scheme, since only orbitals have 
            // updated positions of the QM region, hence vector<Segments*> is 
            // NULL in the QMMachine and the QM region is also printed here
    if(_write_charges){
        qmatoms= orbitals_guess->QMAtoms();
    }
    else{
        QMMInterface qmmface;
        qmatoms=qmmface.Convert(segments);
    }
    
   
    std::vector< ctp::QMAtom* >::iterator it;

         
    for (it = qmatoms.begin(); it < qmatoms.end(); it++) {
        if ((*it)->from_environment) {continue;}
        
        _com_file << setw(3) << (*it)->type.c_str()
                  << setw(12) << setiosflags(ios::fixed) << setprecision(5) << (*it)->x
                  << setw(12) << setiosflags(ios::fixed) << setprecision(5) << (*it)->y
                  << setw(12) << setiosflags(ios::fixed) << setprecision(5) << (*it)->z
                  << endl;
    }
    _com_file << "* \n" << endl;
    _com_file << "%pal\n "  <<  "nprocs " <<  _threads  << "\nend" << "\n" << endl;
       
    // basis set info
    if ( _write_basis_set){
        Elements _elements;

        list<std::string> elements;
        BasisSet bs;
        bs.LoadBasisSet(_basisset_name);
        CTP_LOG(ctp::logDEBUG, *_pLog) << "Loaded Basis Set " << _basisset_name << flush;
        ofstream _el_file;
        std::string _el_file_name = _run_dir + "/" +  "system.bas";
        _el_file.open(_el_file_name.c_str());
        //_com_file << "@" << "system.bas" << endl;
        _el_file << "$DATA" << endl;

        for (it = qmatoms.begin(); it < qmatoms.end(); it++) {
            if ((*it)->from_environment) { continue;}
            
            std::string element_name = (*it)->type;
            list<std::string>::iterator ite;
            ite = find(elements.begin(), elements.end(), element_name);
            if (ite == elements.end()) {
                elements.push_back(element_name);
                Element* element = bs.getElement(element_name);
                _el_file << _elements.getEleFull(element_name) << endl; 
                for (Element::ShellIterator its = element->firstShell(); its != element->lastShell(); its++) {
                    Shell* shell = (*its);

                    string type=shell->getType();
                    // check combined shells

                       for(unsigned i=0;i<type.size();++i){
                           string subtype = string( type, i, 1 );
                          _el_file  << subtype << " " << shell->getSize() << endl; 
                          int _sh_idx =0;
                          for (Shell::GaussianIterator itg = shell->firstGaussian(); itg != shell->lastGaussian(); itg++) {
                              GaussianPrimitive* gaussian = *itg;
                              _sh_idx++;
                              _el_file << " " << _sh_idx << " " << indent(gaussian->decay);
                              _el_file << " " << indent(gaussian->contraction[FindLmax(subtype)]);

                              }
                              _el_file << endl;
                          }
                    }

                 }
             }
        
        _el_file << "STOP\n";
        _el_file.close();


        _com_file << "%basis\n "  << endl;
        _com_file << "GTOName" << " " << "=" << "\"system.bas\";"  << endl;
    } // write_basis set
            
    // ECPs
   /* WRITING ECP INTO system.inp FILE for ORCA**/
    if ( _write_pseudopotentials ){
        std::string pseudopotential_name("ecp");
        _com_file << endl;

        list<std::string> elements;
        elements.push_back("H");
        elements.push_back("He");
        BasisSet ecp;
        ecp.LoadPseudopotentialSet(pseudopotential_name);
        CTP_LOG(ctp::logDEBUG, *_pLog) << "Loaded Pseudopotentials " << pseudopotential_name << flush;



         for (it = qmatoms.begin(); it < qmatoms.end(); it++) {
            if (!(*it)->from_environment) {
                std::string element_name = (*it)->type;


                list<std::string>::iterator ite;
                ite = find(elements.begin(), elements.end(), element_name);
                if (ite == elements.end()) {
                    elements.push_back(element_name);
                    Element* element = ecp.getElement(element_name);
                    _com_file << "\n" << "NewECP" << " " <<  element_name  << endl;
                    _com_file << "N_core" << " " << element->getNcore()  << endl;
                    //lmaxnum2lmaxname 
                    _com_file << "lmax" << " " << getLName(element->getLmax()) << endl;

                    //For Orca the order doesn't matter but let's write it in ascending order
                    // write remaining shells in ascending order s,p,d...
                    for (int i = 0; i < element->getLmax(); i++) {
                        for (Element::ShellIterator its = element->firstShell(); its != element->lastShell(); its++) {
                            Shell* shell = (*its);
                            if (shell->getLmax() == i) {
                            // shell type, number primitives, scale factor
                            _com_file << shell->getType() << " " <<  shell->getSize() << endl;
                            int _sh_idx = 0;
                            for (Shell::GaussianIterator itg = shell->firstGaussian(); itg != shell->lastGaussian(); itg++) {
                                GaussianPrimitive* gaussian = *itg;
                                _sh_idx++;
                                _com_file << _sh_idx << " " << gaussian->decay << " " << gaussian->contraction[0] << " " << gaussian->power << endl;
                            }
                        }
                    }
                }

                for (Element::ShellIterator its = element->firstShell(); its != element->lastShell(); its++) {
                    Shell* shell = (*its);
                    // shell type, number primitives, scale factor
                    if (shell->getLmax() == element->getLmax()) {
                        _com_file << shell->getType() << " " <<  shell->getSize() << endl;
                        // _com_file << shell->getSize() << endl;
                        int _sh_idx = 0;
                        for (Shell::GaussianIterator itg = shell->firstGaussian(); itg != shell->lastGaussian(); itg++) {
                            GaussianPrimitive* gaussian = *itg;
                            _sh_idx++;
                            _com_file << _sh_idx << " " << gaussian->decay << " " << gaussian->contraction[0] << " " << gaussian->power << endl;
                            }
                        }
                    }

                   _com_file << "end\n "  <<  "\n" << endl; 
                }

            }
        }
    } // write pseudopotentials
                
                
    /* END   OF WRITING ECP INTO system.inp FILE for ORCA*************/                
    _com_file << "end\n "  <<  "\n" << endl;   //This end is for the basis set block
             
            
     if(_write_charges){      
        // actual writing of charges    
        _crg_file.open( _crg_file_name_full.c_str() ); 
        int _total_background = 0;
        for (it = qmatoms.begin(); it < qmatoms.end(); it++) {
            if ((*it)->from_environment) {
                if ((*it)->charge != 0.0) _total_background++;
            }
        }

        _crg_file << _total_background << endl;
        for (it = qmatoms.begin(); it < qmatoms.end(); it++) {
            if ((*it)->from_environment) {
                boost::format fmt("%1$+1.7f %2$+1.7f %3$+1.7f %4$+1.7f");
                fmt % (*it)->charge % (*it)->x % (*it)->y % (*it)->z ;
                if ((*it)->charge != 0.0) _crg_file << fmt << endl;
            }
        }
     }            
     
    // writing scratch_dir info
    if ( _scratch_dir != "" ) {
        
        CTP_LOG(ctp::logDEBUG,*_pLog) << "Setting the scratch dir to " << _scratch_dir + temp_suffix << flush;

        // boost::filesystem::create_directories( _scratch_dir + temp_suffix );
        std::string _temp( "scratch_dir " + _scratch_dir + temp_suffix + "\n" );
        //_com_file << _temp ;
    }

    _com_file << _options << "\n";

    
    if ( _write_guess ) { 
       
            throw std::runtime_error( "Not implemented in orca");
    }   

    _com_file << endl;
    _com_file.close();
    
      // and now generate a shell script to run both jobs, if neccessary
    CTP_LOG(ctp::logDEBUG, *_pLog) << "Setting the scratch dir to " << _scratch_dir + temp_suffix << flush;

    _scratch_dir = scratch_dir_backup + temp_suffix;
            
            //boost::filesystem::create_directories(_scratch_dir + temp_suffix);
            //std::string _temp("scratch_dir " + _scratch_dir + temp_suffix + "\n");
            //_com_file << _temp;
    WriteShellScript();
    _scratch_dir = scratch_dir_backup;
   
    return true;
}



bool Orca::WriteShellScript() {
    ofstream _shell_file;
    
    std::string _shell_file_name_full = _run_dir + "/" + _shell_file_name;
            
    _shell_file.open ( _shell_file_name_full.c_str() );

    _shell_file << "#!/bin/bash" << endl ;
    _shell_file << "mkdir -p " << _scratch_dir << endl;
    
    if ( _threads == 1 ){ 
        _shell_file << _executable << " " << _input_file_name << " > " <<  _log_file_name  << endl;//" 2> run.error" << endl;    
    } else {
        _shell_file << _executable << " "<< _input_file_name << " > " <<  _log_file_name <<endl;// " 2> run.error" << endl;    
    }
    _shell_file.close();
    
    return true;   
}

/**
 * Runs the Orca job. 
 */
bool Orca::Run()
{

    CTP_LOG(ctp::logDEBUG,*_pLog) << "Running Orca job" << flush;
    
    if (std::system(NULL)) {
        
        // Orca overrides input information, if *.db and *.movecs files are present
        // better trash the old version
        std::string file_name = _run_dir + "/system.db";
        remove ( file_name.c_str() );
        file_name = _run_dir + "/" + _log_file_name;
        remove ( file_name.c_str() );
        file_name = _run_dir + "/" + _orb_file_name;
        //remove ( file_name.c_str() );
               
        std::string _command;
        if ( _threads == 1 ) {
            _command = "cd " + _run_dir + "; sh " + _shell_file_name;
        } else {
            _command = "cd " + _run_dir + "; sh " + _shell_file_name;
        }
        //CTP_LOG(logDEBUG,*_pLog) << _command << flush;
       int check=std::system(_command.c_str());
        if (check==-1){
            CTP_LOG(ctp::logERROR, *_pLog) << _input_file_name << " failed to start" << flush;
            return false;
        }
       
        if ( CheckLogFile() ) {
            CTP_LOG(ctp::logDEBUG,*_pLog) << "Finished Orca job" << flush;
            return true;
        } else {
            CTP_LOG(ctp::logDEBUG,*_pLog) << "Orca job failed" << flush;
        }
    }
    else {
        CTP_LOG(ctp::logERROR,*_pLog) << _input_file_name << " failed to start" << flush; 
        return false;
    }
    
    return true;
}

/**
 * Cleans up after the Orca job
 */
void Orca::CleanUp() {
    
    // cleaning up the generated files
    if ( _cleanup.size() != 0 ) {
        Tokenizer tok_cleanup(_cleanup, ",");
        std::vector <std::string> _cleanup_info;
        tok_cleanup.ToVector(_cleanup_info);
        
        std::vector<std::string> ::iterator it;
               
        for (it = _cleanup_info.begin(); it != _cleanup_info.end(); ++it) {
            if ( *it == "inp" ) {
                std::string file_name = _run_dir + "/" + _input_file_name;
                remove ( file_name.c_str() );
            }
            
            if ( *it == "bas" ) {
                std::string file_name = _run_dir +"/system.bas";
                remove ( file_name.c_str() );
            }
            
            if ( *it == "log" ) {
                std::string file_name = _run_dir + "/" + _log_file_name;
                remove ( file_name.c_str() );
            }

           if ( *it == "gbw" ) {
                std::string file_name = _run_dir + "/" + _orb_file_name;
                remove ( file_name.c_str() );
            }
            
            if ( *it == "ges" ) {
                std::string file_name = _run_dir + "/system.ges" ;
                remove ( file_name.c_str() );
            }    
            if ( *it == "prop" ) {
                std::string file_name = _run_dir + "/system.prop" ;
                remove ( file_name.c_str() );
            }     
        }
    }
return; 
}


bool Orca::ParseLogFile( Orbitals* _orbitals )
{
    const double _conv_Hrt_eV = tools::conv::hrt2ev;
    
    _orbitals->setQMpackage("orca");
    _orbitals->setDFTbasis( _basisset_name);
    CTP_LOG(ctp::logDEBUG,*_pLog) << "Parsing " << _log_file_name << flush;
      // return true;
  std::string _log_file_name_full =  _run_dir + "/" + _log_file_name;
   // check if LOG file is complete
   if ( !CheckLogFile() ) return false;
    
    //std::map <int, std::std::vector<double> > _coefficients;
    std::map <int, double> _energies;
    std::map <int, double> _occ;
    
    std::string _line;
    unsigned _levels = 0;
    //unsigned _level;
    //unsigned _basis_size = 0;
    int _number_of_electrons = 0;
    //bool _has_basis_dim = false;
               
    std::vector<std::string> results;    
    
    std::ifstream _input_file( _log_file_name_full.c_str() );
    
    if (_input_file.fail()) {
        CTP_LOG( ctp::logERROR, *_pLog ) << "File " << _log_file_name_full << " not found " << flush;
        return false;
    } else {
        CTP_LOG(ctp::logDEBUG, *_pLog) << "Reading Coordinates and occupationnumbers and energies from " << _log_file_name_full << flush;
    }
    
      
          //Coordinates of the final configuration depending on whether it is an optimization or not
         
        
        
    while (_input_file) {
        getline(_input_file,_line);
        boost::trim(_line);
        
        
        
        if ( _is_optimization ){
              throw runtime_error("Not implemented yet!");
        }
        bool _found_optimization=true;
        
        std::string::size_type coordinates_pos = _line.find("CARTESIAN COORDINATES (ANGSTROEM)");
        
        if ( _found_optimization && coordinates_pos != std::string::npos) {
            CTP_LOG(ctp::logDEBUG,*_pLog) << "Getting the coordinates" << flush;
            
            //_has_coordinates = true;
            bool _has_QMAtoms = _orbitals->hasQMAtoms();

            // three garbage lines
            getline(_input_file, _line);
            // now starts the data in format
            // _id type Qnuc x y z 
            std::vector<std::string> _row;
            getline(_input_file, _line);
            boost::trim( _line );

            boost::algorithm::split( _row , _line, boost::is_any_of("\t "), boost::algorithm::token_compress_on); 
            int nfields =  _row.size();

            int atom_id = 0;   
            while ( nfields == 4 ) {
                //int atom_id = boost::lexical_cast< int >( _row.at(0) );
                //int atom_number = boost::lexical_cast< int >( _row.at(0) );
                std::string _atom_type = _row.at(0);
                double _x =  boost::lexical_cast<double>( _row.at(1) );
                double _y =  boost::lexical_cast<double>( _row.at(2) );
                double _z =  boost::lexical_cast<double>( _row.at(3) );
                //if ( tools::globals::verbose ) cout << "... ... " << atom_id << " " << atom_type << " " << atom_charge << endl;
                getline(_input_file, _line);
                boost::trim( _line );
                boost::algorithm::split( _row , _line, boost::is_any_of("\t "), boost::algorithm::token_compress_on);  
                nfields =  _row.size();

                if ( _has_QMAtoms == false ) {
                    _orbitals->AddAtom( _atom_type, _x, _y, _z );
                } else {
                                       
                    ctp::QMAtom* pAtom = _orbitals->_atoms.at( atom_id  );
                    pAtom->type = _atom_type;
                    pAtom->x = _x;
                    pAtom->y = _y;
                    pAtom->z = _z;
                    atom_id++;
                }
                    
            }
      
        }
              
        std::string::size_type energy_pos = _line.find("Total Energy");
       if (energy_pos != std::string::npos) {
            //cout << _line << endl;
            boost::algorithm::split(results, _line, boost::is_any_of(" "), boost::algorithm::token_compress_on);
            std::string _energy = results[3];
            boost::trim( _energy );
            //cout << _energy << endl; 
            _orbitals->setQMEnergy ( _conv_Hrt_eV * boost::lexical_cast<double>(_energy) );
            CTP_LOG(ctp::logDEBUG, *_pLog) << (boost::format("QM energy[eV]: %4.6f ") % _orbitals->getQMEnergy()).str() << flush;
            // _orbitals->_has_qm_energy = true;
            }
        
        /* Check for ScaHFX = factor of HF exchange included in functional */
              std::string::size_type HFX_pos = _line.find("Fraction HF Exchange ScalHFX");
              if (HFX_pos != std::string::npos) {
                  boost::algorithm::split(results, _line, boost::is_any_of(" "), boost::algorithm::token_compress_on);
                  double _ScaHFX = boost::lexical_cast<double>(results.back());
                  _orbitals->setScaHFX(_ScaHFX);
                  CTP_LOG(ctp::logDEBUG, *_pLog) << "DFT with " << _ScaHFX << " of HF exchange!" << flush;
              }
     
    //Finding Basis Dimension, the number of energy levels
        std::string::size_type dim_pos = _line.find("Basis Dimension");
        if (dim_pos != std::string::npos) {
                
            boost::algorithm::split(results, _line, boost::is_any_of(" "), boost::algorithm::token_compress_on);
            //_has_basis_dim = true;
            std::string _dim = results[4];  //The 4th element of results vector is the Basis Dim
            boost::trim( _dim );
            _levels = boost::lexical_cast<int>(_dim);
            //cout <<  boost::lexical_cast<int>(_dim) << endl;
            //_basis_size = _levels;
              CTP_LOG(ctp::logDEBUG,*_pLog) << "Basis Dimension: " << _levels << flush;
              CTP_LOG( ctp::logDEBUG, *_pLog ) << "Energy levels: " << _levels << flush;
        }
        /********************************************************/
         
        
        std::string::size_type OE_pos = _line.find("ORBITAL ENERGIES");
        if (OE_pos != std::string::npos) {
        getline(_input_file,_line);
        getline(_input_file,_line);
        getline(_input_file,_line);
        if (_line.find("E(Eh)")==std::string::npos){
            CTP_LOG(ctp::logDEBUG,*_pLog) << "Warning: Orbital Energies not found in log file" << flush;
        }
        for (unsigned i=0; i<_levels ; i++){
            getline(_input_file,_line);
            boost::trim( _line );
            boost::algorithm::split(results, _line, boost::is_any_of(" "), boost::algorithm::token_compress_on);
           
                    
            std::string _no =results[0];
            
            boost::trim( _no );
            unsigned levelnumber= boost::lexical_cast<unsigned>(_no);
            if (levelnumber!=i){
                CTP_LOG(ctp::logDEBUG,*_pLog) << "Have a look at the orbital energies something weird is going on" << flush;
            }
            std::string _oc = results[1];
            boost::trim( _oc );
            double occ =boost::lexical_cast<double>(_oc);
            // We only count alpha electrons, each orbital must be empty or doubly occupied
            if (occ==2){
                _number_of_electrons++;
                _occ[i]=occ;
            }
            else if (occ==0) {
                _occ[i]=occ;
            }
            else {
                throw runtime_error("Only empty or doubly occupied orbitals are allowed not running the right kind of DFT calculation");
            }
            
            std::string _e = results[2];
            boost::trim( _e );
            _energies [i] = boost::lexical_cast<double>(_e) ;
        }
    }
            /*
                 *  Partial charges from the input file
                 */
        std::string::size_type charge_pos = _line.find("CHELPG Charges");

        if (charge_pos != std::string::npos && _get_charges) {
                CTP_LOG(ctp::logDEBUG, *_pLog) << "Getting charges" << flush;
                //_has_charges = true;
                getline(_input_file, _line);
                //getline(_input_file, _line);

                bool _has_atoms = _orbitals->hasQMAtoms();

                std::vector<std::string> _row;
                getline(_input_file, _line);
                boost::trim(_line);
                //cout << _line << endl;
                boost::algorithm::split(_row, _line, boost::is_any_of("\t "), boost::algorithm::token_compress_on);
                int nfields = _row.size();
                //cout << _row.size() << endl;

                while (nfields == 4) {
                    int atom_id = boost::lexical_cast< int >(_row.at(0));
                    atom_id++;
                    //int atom_number = boost::lexical_cast< int >(_row.at(0));
                    std::string atom_type = _row.at(1);
                    double atom_charge = boost::lexical_cast< double >(_row.at(3));
                    if ( tools::globals::verbose ) cout << "... ... " << atom_id << " " << atom_type << " " << atom_charge << endl;
                    getline(_input_file, _line);
                    boost::trim(_line);
                    boost::algorithm::split(_row, _line, boost::is_any_of("\t "), boost::algorithm::token_compress_on);
                    nfields = _row.size();

                    if (_has_atoms == false) {
                        _orbitals->AddAtom(atom_type, 0, 0, 0, atom_charge);
                    //_orbitals->_has_atoms = true;
                    }else {
                        ctp::QMAtom* pAtom = _orbitals->_atoms.at(atom_id - 1);
                        pAtom->type = atom_type;
                        pAtom->charge = atom_charge;
                    }
{'fraction_non_alphanumeric': 0.14308112393200034, 'fraction_numerical': 0.003294283449308553, 'mean_word_length': 1.769748719199805, 'pattern_counts': {'":': 0, '<': 1043, '<?xml version=': 0, '>': 306, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7563268', 'n_tokens_mistral': 17815, 'n_tokens_neox': 17028, 'n_words': 3607}