//
// mpqc.cc
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Edward Seidl <seidl@janed.com>
// Maintainer: LPS
//
// This file is part of MPQC.
//
// MPQC 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; either version 2, or (at your option)
// any later version.
//
// MPQC 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.
//
// You should have received a copy of the GNU General Public License
// along with the MPQC; see the file COPYING.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//

// This is needed to make GNU extensions available, such as
// feenableexcept and fedisableexcept.
#ifndef _GNU_SOURCE
#  define _GNU_SOURCE
#endif

#ifdef HAVE_CONFIG_H
#include <scconfig.h>
#endif

#ifdef HAVE_TIME_H
#include <time.h>
#endif

#include <scdirlist.h>

#include <new>
#include <stdexcept>
#include <string.h>
#include <unistd.h>
#include <sys/stat.h>
#include <fstream>

#include <boost/bind.hpp>
#include <boost/function.hpp>

#include <scconfig.h>
#ifdef HAVE_SSTREAM
#  include <sstream>
#else
#  include <strstream.h>
#endif

#ifdef HAVE_SYS_RESOURCE_H
#  include <sys/resource.h>
#endif
#ifdef HAVE_SYS_TIME_H
#  include <sys/time.h>
#endif

#include <util/options/GetLongOpt.h>
#include <util/class/scexception.h>
#include <util/misc/newstring.h>
#include <util/keyval/keyval.h>
#include <util/state/state_bin.h>
#include <util/group/message.h>
#include <util/group/memory.h>
#include <util/group/mstate.h>
#include <util/group/thread.h>
#include <util/group/pregtime.h>
#include <util/misc/bug.h>
#include <util/misc/formio.h>
#include <util/misc/exenv.h>
#ifdef HAVE_CHEMISTRY_CCA
  #include <util/misc/ccaenv.h>
#endif
#include <util/render/render.h>

#include <math/optimize/opt.h>

#include <chemistry/molecule/coor.h>
#include <chemistry/molecule/energy.h>
#include <chemistry/molecule/molfreq.h>
#include <chemistry/molecule/fdhess.h>
#include <chemistry/molecule/formula.h>
#include <chemistry/qc/wfn/wfn.h>

// Force linkages:
#include <util/group/linkage.h>
#include <chemistry/qc/wfn/linkage.h>
#include <chemistry/qc/scf/linkage.h>
#include <chemistry/qc/dft/linkage.h>
#include <chemistry/qc/mbpt/linkage.h>
#ifdef HAVE_SC_SRC_LIB_CHEMISTRY_QC_MBPTR12
#  include <chemistry/qc/mbptr12/linkage.h>
#endif
#ifdef HAVE_SC_SRC_LIB_CHEMISTRY_QC_CINTS
#  include <chemistry/qc/cints/linkage.h>
#endif
//#include <chemistry/qc/psi/linkage.h>
#include <util/state/linkage.h>
#ifdef HAVE_SC_SRC_LIB_CHEMISTRY_QC_CC
#  include <chemistry/qc/cc/linkage.h>
#endif
#ifdef HAVE_SC_SRC_LIB_CHEMISTRY_QC_PSI
#  include <chemistry/qc/psi/linkage.h>
#endif
#ifdef HAVE_SC_SRC_LIB_CHEMISTRY_QC_INTCCA
#  include <chemistry/qc/intcca/linkage.h>
#endif

#ifdef HAVE_MPI
#define MPICH_SKIP_MPICXX
#include <mpi.h>
#include <util/group/messmpi.h>
#endif

#ifdef HAVE_JOBMARKET
// include headers that implement a archive in simple text format
// otherwise BOOST_CLASS_EXPORT_IMPLEMENT has no effect
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <boost/bind.hpp>
#include <boost/tokenizer.hpp>

#include "JobMarket/Results/FragmentResult.hpp"
#include "JobMarket/poolworker_main.hpp"

#include "chemistry/qc/scf/scfops.h"

#ifdef HAVE_MPQCDATA
#include "Jobs/MPQCJob.hpp"
#include "Fragmentation/Summation/Containers/MPQCData.hpp"

#include <chemistry/qc/basis/obint.h>
#include <chemistry/qc/basis/symmint.h>
#endif

#include <algorithm>
#include <stdlib.h>
#endif

using namespace std;
using namespace sc;

#include "mpqcin.h"

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

const KeyValValueboolean truevalue(1), falsevalue(0);
const char *devnull = "/dev/null";


static void
trash_stack_b(int &i, char *&ichar)
{
  char stack;
  ichar = &stack;
  ichar -= 10;
  for (i=0; i<1000; i++) {
    *ichar-- = 0xfe;
  }
}

static void
trash_stack()
{
  int i;
  char *ichar;
  trash_stack_b(i,ichar);
}

static void
clean_up(void)
{
  MemoryGrp::set_default_memorygrp(0);
  ThreadGrp::set_default_threadgrp(0);
  SCMatrixKit::set_default_matrixkit(0);
  Integral::set_default_integral(0);
  RegionTimer::set_default_regiontimer(0);
}

static void
final_clean_up(void)
{
  MessageGrp::set_default_messagegrp(0);
}

#include <signal.h>

#ifdef HAVE_FENV_H
#  include <fenv.h>
#endif

static void
print_unseen(const Ref<ParsedKeyVal> &parsedkv,
             const char *input)
{
  if (parsedkv->have_unseen()) {
    ExEnv::out0() << endl;
    ExEnv::out0() << indent
         << "The following keywords in \"" << input << "\" were ignored:"
         << endl;
    ExEnv::out0() << incindent;
    parsedkv->print_unseen(ExEnv::out0());
    ExEnv::out0() << decindent;
  }
}

double EvaluateDensity(
    SCVector3 &r,
    Ref<Integral> &intgrl,
    GaussianBasisSet::ValueData &vdat,
    Ref<Wavefunction> &wfn);

/** Places all known options into \a options and parses them from argc,argv.
 *
 * \param options options structure
 * \param argc argument count
 * \param argv argument array
 * \return return value by GetLongOpt::parse() function
 */
int ParseOptions(
    GetLongOpt &options,
    int argc,
    char **argv)
{
  options.usage("[options] [filename]");
  options.enroll("f", GetLongOpt::MandatoryValue,
                 "the name of an object format input file", 0);
  options.enroll("o", GetLongOpt::MandatoryValue,
                 "the name of the output file", 0);
  options.enroll("n", GetLongOpt::MandatoryValue,
                 "listen for incoming object format input files", 0);
  options.enroll("messagegrp", GetLongOpt::MandatoryValue,
                 "which message group to use", 0);
  options.enroll("threadgrp", GetLongOpt::MandatoryValue,
                 "which thread group to use", 0);
  options.enroll("memorygrp", GetLongOpt::MandatoryValue,
                 "which memory group to use", 0);
  options.enroll("integral", GetLongOpt::MandatoryValue,
                 "which integral evaluator to use", 0);
  options.enroll("l", GetLongOpt::MandatoryValue, "basis set limit", "0");
  options.enroll("W", GetLongOpt::MandatoryValue,
                 "set the working directory", ".");
  options.enroll("c", GetLongOpt::NoValue, "check input then exit", 0);
  options.enroll("v", GetLongOpt::NoValue, "print the version number", 0);
  options.enroll("w", GetLongOpt::NoValue, "print the warranty", 0);
  options.enroll("L", GetLongOpt::NoValue, "print the license", 0);
  options.enroll("k", GetLongOpt::NoValue, "print key/value assignments", 0);
  options.enroll("i", GetLongOpt::NoValue, "convert simple to OO input", 0);
  options.enroll("d", GetLongOpt::NoValue, "debug", 0);
  options.enroll("h", GetLongOpt::NoValue, "print this message", 0);
  options.enroll("cca-path", GetLongOpt::OptionalValue,
                 "cca component path", "");
  options.enroll("cca-load", GetLongOpt::OptionalValue,
                 "cca components to load", "");

  int optind = options.parse(argc, argv);

  return optind;
}

/** Checks for each known option and acts accordingly.
 *
 * \param options option structure
 * \param *output name of outputfile on return
 * \param *outstream open output stream on return
 */
void ComputeOptions(
    GetLongOpt &options,
    const char *&output,
    ostream *&outstream)
{
  output = options.retrieve("o");
  outstream = 0;
  if (output != 0) {
    outstream = new ofstream(output);
    ExEnv::set_out(outstream);
  }

  if (options.retrieve("h")) {
    ExEnv::out0()
         << indent << "MPQC version " << SC_VERSION << endl
         << indent << "compiled for " << TARGET_ARCH << endl
         << SCFormIO::copyright << endl;
    options.usage(ExEnv::out0());
    exit(0);
  }

  if (options.retrieve("v")) {
    ExEnv::out0()
         << indent << "MPQC version " << SC_VERSION << endl
         << indent << "compiled for " << TARGET_ARCH << endl
         << SCFormIO::copyright;
    exit(0);
  }

  if (options.retrieve("w")) {
    ExEnv::out0()
         << indent << "MPQC version " << SC_VERSION << endl
         << indent << "compiled for " << TARGET_ARCH << endl
         << SCFormIO::copyright << endl
         << SCFormIO::warranty;
    exit(0);
  }

  if (options.retrieve("L")) {
    ExEnv::out0()
         << indent << "MPQC version " << SC_VERSION << endl
         << indent << "compiled for " << TARGET_ARCH << endl
         << SCFormIO::copyright << endl
         << SCFormIO::license;
    exit(0);
  }

  if (options.retrieve("d"))
    SCFormIO::set_debug(1);

  // set the working dir
  if (strcmp(options.retrieve("W"),"."))
    int retval = chdir(options.retrieve("W"));

  // check that n and f/o are not given at the same time
  if ((options.retrieve("n")) && ((options.retrieve("f")) || (options.retrieve("o")))) {
    throw invalid_argument("-n must not be given with -f or -o");
  }
}

/** Temporary structure for storing information from command-line
 *
 * This structure has been introduced to gather the various calls to GetLongOpts
 * at one (initial) place and to abstract it from the source of command-lines.
 * This temporary object can be set by other means, too. I.e. we become
 * independent of usage in command-line programs.
 */
struct OptionValues {
  const char *keyvalue; // option "k"
  const char *debug; // option ""
  int limit; // option "l"
  const char *check; // option "c"
  const char *simple_input; // option "i"
  string executablename;

#ifdef HAVE_CHEMISTRY_CCA
  string cca_load;  // option "cca-path"
  string cc_path; // option "cca-load"
#endif
};

/** Parse remainder options not treated by ComputeOptions() into temporary storage.
 *
 * \param options option structure to obtain values from
 * \param values remaining option values which are processed later and now
 *        stored in a temporary structure
 */
void parseRemainderOptions(
    GetLongOpt &options,
    struct OptionValues &values,
    int argc,
    char **argv)
{
  values.keyvalue = options.retrieve("k");
  values.debug = options.retrieve("d");
  values.limit = atoi(options.retrieve("l"));
  values.check = options.retrieve("c");
  values.simple_input = options.retrieve("i");
  values.executablename = argv[0];

#ifdef HAVE_CHEMISTRY_CCA
  values.cca_load = options.retrieve("cca-load");
  values.cca_path = options.retrieve("cca-path");
#endif
}

/** Sets object and generic input file names.
 *
 * \param object_input filename of object-oriented input
 * \param generic_input filename of generic input
 * \param options (command-line)option structure
 * \param argc argument count
 * \param argv argument array
 */
void getInputFileNames(
    const char *&object_input,
    const char *&generic_input,
    GetLongOpt &options,
    int optind,
    int argc,
    char **argv)
{
  // initialize keyval input
  object_input = options.retrieve("f");
  generic_input = 0;
  if (argc - optind == 0) {
    generic_input = 0;
  }
  else if (argc - optind == 1) {
    generic_input = argv[optind];
  }
  else if (!options.retrieve("n")) {
    options.usage();
    throw invalid_argument("extra arguments given");
  }

  if (object_input == 0 && generic_input == 0) {
    generic_input = "mpqc.in";
    }
  else if (object_input && !options.retrieve("n") && generic_input) {
    options.usage();
    throw invalid_argument("only one of -f and a file argument can be given");
    }
}

/** Gets the MPI Message group.
 *
 * \param grp reference to obtained group
 * \param argc argument count
 * \param argv argument array
 */
void getMessageGroup(
    Ref<MessageGrp> &grp,
    int argc,
    char **argv)
{
#if defined(HAVE_MPI) && defined(ALWAYS_USE_MPI)
  grp = new MPIMessageGrp(&argc, &argv);
#endif
  if (grp.null()) grp = MessageGrp::initial_messagegrp(argc, argv);
  if (grp.nonnull())
    MessageGrp::set_default_messagegrp(grp);
  else
    grp = MessageGrp::get_default_messagegrp();
}

/** Sets the base name of output files.
 *
 * \param input input file name
 * \param output output file name
 */
void setOutputBaseName(const char *input, const char *output)
{
  const char *basename_source;
  if (output) basename_source = output;
  else        basename_source = input;
  const char *baseprefix = ::strrchr(basename_source, '.');
  int nfilebase = 1;
  if (baseprefix == NULL) {
    std::cerr << "setOutputBaseName() - ERROR: basename_source "
        << basename_source << " contains no dot (.)." << std::endl;
    nfilebase = ::strlen(basename_source);
  } else
    nfilebase = (int) (baseprefix - basename_source);
  char *basename = new char[nfilebase + 1];
  strncpy(basename, basename_source, nfilebase);
  basename[nfilebase] = '\0';
  SCFormIO::set_default_basename(basename);
  delete[] basename;
}

/** Prints current key values.
 *
 * \param keyval key value structure
 * \param opt optimization structure
 * \param molname name of molecule
 * \param restartfile name of restartfile
 */
void printOptions(
    Ref<KeyVal> &keyval,
    Ref<Optimize> &opt,
    const char *molname,
    const char *restartfile)
{
  int restart = keyval->booleanvalue("restart",truevalue);

  int checkpoint = keyval->booleanvalue("checkpoint",truevalue);

  int savestate = keyval->booleanvalue("savestate",truevalue);

  int do_energy = keyval->booleanvalue("do_energy",truevalue);

  int do_grad = keyval->booleanvalue("do_gradient",falsevalue);

  int do_opt = keyval->booleanvalue("optimize",truevalue);

  int do_pdb = keyval->booleanvalue("write_pdb",falsevalue);

  int print_mole = keyval->booleanvalue("print_mole",truevalue);

  int print_timings = keyval->booleanvalue("print_timings",truevalue);

  // sanity checks for the benefit of reasonable looking output
  if (opt.null()) do_opt=0;

  ExEnv::out0() << endl << indent
       << "MPQC options:" << endl << incindent
       << indent << "matrixkit     = <"
       << SCMatrixKit::default_matrixkit()->class_name() << ">" << endl
       << indent << "filename      = " << molname << endl
       << indent << "restart_file  = " << restartfile << endl
       << indent << "restart       = " << (restart ? "yes" : "no") << endl
       << indent << "checkpoint    = " << (checkpoint ? "yes" : "no") << endl
       << indent << "savestate     = " << (savestate ? "yes" : "no") << endl
       << indent << "do_energy     = " << (do_energy ? "yes" : "no") << endl
       << indent << "do_gradient   = " << (do_grad ? "yes" : "no") << endl
       << indent << "optimize      = " << (do_opt ? "yes" : "no") << endl
       << indent << "write_pdb     = " << (do_pdb ? "yes" : "no") << endl
       << indent << "print_mole    = " << (print_mole ? "yes" : "no") << endl
       << indent << "print_timings = " << (print_timings ? "yes" : "no")
       << endl << decindent;

}

/** Saves the current state to checkpoint file.
 *
 * \param keyval key value structure
 * \param opt optimization structure
 * \param grp message group
 * \param mole MolecularEnergy object
 * \param molname name of molecule
 * \param ckptfile name of check point file
 */
void saveState(
    char *wfn_file,
    int savestate,
    Ref<Optimize> &opt,
    Ref<MessageGrp> &grp,
    Ref<MolecularEnergy> &mole,
    char *&molname,
    char *&ckptfile)
{
  // function stuff
  if (savestate) {
    if (opt.nonnull()) {
      if (grp->me() == 0) {
        ckptfile = new char[strlen(molname)+6];
        sprintf(ckptfile,"%s.ckpt",molname);
      }
      else {
        ckptfile = new char[strlen(devnull)+1];
        strcpy(ckptfile, devnull);
      }

      StateOutBin so(ckptfile);
      SavableState::save_state(opt.pointer(),so);
      so.close();

      delete[] ckptfile;
    }

    if (mole.nonnull()) {
      if (grp->me() == 0) {
        if (wfn_file == 0) {
          wfn_file = new char[strlen(molname)+6];
          sprintf(wfn_file,"%s.wfn",molname);
        }
      }
      else {
        delete[] wfn_file;
        wfn_file = new char[strlen(devnull)+1];
        strcpy(wfn_file, devnull);
      }

      StateOutBin so(wfn_file);
      SavableState::save_state(mole.pointer(),so);
      so.close();

    }
  }
  delete[] wfn_file;
}

/** Sets up indentation and output modes.
 *
 * \param grp message group
 */
void setupSCFormIO(
    Ref<MessageGrp> &grp
    )
{
  SCFormIO::setindent(ExEnv::outn(), 2);
  SCFormIO::setindent(ExEnv::errn(), 2);
  SCFormIO::setindent(cout, 2);
  SCFormIO::setindent(cerr, 2);

  SCFormIO::set_printnode(0);
  if (grp->n() > 1)
    SCFormIO::init_mp(grp->me());
}

/** Initialises the timer.
 *
 * \param grp message group
 * \param keyval key value structure
 * \param tim timing structure
 */
void initTimings(
    Ref<MessageGrp> &grp,
    Ref<KeyVal> &keyval,
    Ref<RegionTimer> &tim
    )
{
  grp->sync(); // make sure nodes are sync'ed before starting timings
  if (keyval->exists("timer")) tim << keyval->describedclassvalue("timer");
  else                         tim = new ParallelRegionTimer(grp,"mpqc",1,1);
  RegionTimer::set_default_regiontimer(tim);

  if (tim.nonnull()) tim->enter("input");
}

/** Prints the header of the output.
 *
 * \param tim timing structure
 */
void makeAnnouncement(
    Ref<RegionTimer> &tim
    )
{
  const char title1[] = "MPQC: Massively Parallel Quantum Chemistry";
  int ntitle1 = sizeof(title1);
  const char title2[] = "Version " SC_VERSION;
  int ntitle2 = sizeof(title2);
  ExEnv::out0() << endl;
  ExEnv::out0() << indent;
  for (int i=0; i<(80-ntitle1)/2; i++) ExEnv::out0() << ' ';
  ExEnv::out0() << title1 << endl;
  ExEnv::out0() << indent;
  for (int i=0; i<(80-ntitle2)/2; i++) ExEnv::out0() << ' ';
  ExEnv::out0() << title2 << endl << endl;

  const char *tstr = 0;
#if defined(HAVE_TIME) && defined(HAVE_CTIME)
  time_t t;
  time(&t);
  tstr = ctime(&t);
#endif
  if (!tstr) {
    tstr = "UNKNOWN";
  }

  ExEnv::out0()
       << indent << scprintf("Machine:    %s", TARGET_ARCH) << endl
       << indent << scprintf("User:       %s@%s",
                             ExEnv::username(), ExEnv::hostname()) << endl
       << indent << scprintf("Start Time: %s", tstr) << endl;
}

/** Parse the input configuration from char array into keyvalue container.
 *
 * \param parsedkv key value container to foll
 * \param values temporary options value structure
 * \param in_char_array char array with input file
 * \param use_simple_input whether the format in \a in_char_array is simple (1)
 *        or object-oriented (0)
 */
void parseIntoKeyValue(
    Ref<ParsedKeyVal> &parsedkv,
    struct OptionValues &values,
    char *&in_char_array,
    int use_simple_input)
{
  if (use_simple_input) {
    MPQCIn mpqcin;
    char *simple_input_text = mpqcin.parse_string(in_char_array);
    if (values.simple_input) {
      ExEnv::out0() << "Generated object-oriented input file:" << endl
                   << simple_input_text
                   << endl;
      exit(0);
    }
    parsedkv = new ParsedKeyVal();
    parsedkv->parse_string(simple_input_text);
    delete[] simple_input_text;
  } else {
    parsedkv = new ParsedKeyVal();
    parsedkv->parse_string(in_char_array);
  }
}

/** Parse the input file into the key value container.
 *
 * \param grp message group
 * \param parsedkev keyvalue container on return
 * \param values (command-line) options structure
 * \param input input file name
 * \param generic_input filename of generic input
 * \param in_char_array char array with input file's contents on return
 * \param use_simple_input whether the file contents is in simple format (1)
 *        or object-oriented (0)
 */
void parseInputfile(
    Ref<MessageGrp> &grp,
    Ref<ParsedKeyVal> &parsedkv,
    struct OptionValues &values,
    const char *&input,
    const char *&generic_input,
    char *&in_char_array,
    int &use_simple_input
    )
{
  // read the input file on only node 0
  if (grp->me() == 0) {
    ifstream is(input);
#ifdef HAVE_SSTREAM
    ostringstream ostrs;
    is >> ostrs.rdbuf();
    int n = 1 + strlen(ostrs.str().c_str());
    in_char_array = strcpy(new char[n],ostrs.str().c_str());
#else
    ostrstream ostrs;
    is >> ostrs.rdbuf();
    ostrs << ends;
    in_char_array = ostrs.str();
    int n = ostrs.pcount();
#endif
    grp->bcast(n);
    grp->bcast(in_char_array, n);
    }
  else {
    int n;
    grp->bcast(n);
    in_char_array = new char[n];
    grp->bcast(in_char_array, n);
    }

  if (generic_input && grp->me() == 0) {
    MPQCIn mpqcin;
    use_simple_input = mpqcin.check_string(in_char_array);
  }
  else {
    use_simple_input = 0;
  }
  grp->bcast(use_simple_input);
}

/** Get the thread group.
 *
 * \param keyval keyvalue container
 * \param thread thread group on return
 * \param argc argument count
 * \param argv argument array
 */
void getThreadGroup(
    Ref<KeyVal> &keyval,
    Ref<ThreadGrp> &thread,
    int argc,
    char **argv)
{
  //first try the commandline and environment
  thread = ThreadGrp::initial_threadgrp(argc, argv);

  // if we still don't have a group, try reading the thread group
  // from the input
  if (thread.null()) {
    thread << keyval->describedclassvalue("thread");
  }

  if (thread.nonnull())
    ThreadGrp::set_default_threadgrp(thread);
  else
    thread = ThreadGrp::get_default_threadgrp();
}

/** Get the memory group.
 *
 * \param keyval keyvalue container
 * \param memory memory group on return
 * \param argc argument count
 * \param argv argument array
 */
void getMemoryGroup(
    Ref<KeyVal> &keyval,
    Ref<MemoryGrp> &memory,
    int argc,
    char **argv)
{
  // first try the commandline and environment
  memory = MemoryGrp::initial_memorygrp(argc, argv);

  // if we still don't have a group, try reading the memory group
  // from the input
  if (memory.null()) {
    memory << keyval->describedclassvalue("memory");
  }

  if (memory.nonnull())
    MemoryGrp::set_default_memorygrp(memory);
  else
    memory = MemoryGrp::get_default_memorygrp();
}

/** Prepares CCA component if available.
 *
 * \param keyval keyvalue container
 * \param values parsed (command-line) options
 */
void prepareCCA(
    Ref<KeyVal> &keyval,
    struct OptionValues &values
    )
{
#ifdef HAVE_CHEMISTRY_CCA
  // initialize cca framework
  KeyValValuestring emptystring("");
  bool do_cca = keyval->booleanvalue("do_cca",falsevalue);

  string cca_path(values.cca_path);
  string cca_load(values.cca_load);
  if(cca_path.size()==0)
    cca_path = keyval->stringvalue("cca_path",emptystring);
  if(cca_load.size()==0)
    cca_load = keyval->stringvalue("cca_load",emptystring);

  if( !do_cca && (cca_load.size() > 0 || cca_path.size() > 0) )
    do_cca = true;

  if(cca_path.size()==0) {
    #ifdef CCA_PATH
      cca_path = CCA_PATH;
    #endif
  }
  if(cca_load.size()==0) {
    cca_load += "MPQC.IntegralEvaluatorFactory";
  }

  if( cca_load.size() > 0 && cca_path.size() > 0 && do_cca ) {
    string cca_args = "--path " + cca_path + " --load " + cca_load;
    ExEnv::out0() << endl << indent << "Initializing CCA framework with args: "
                  << endl << indent << cca_args << endl;
    CCAEnv::init( cca_args );
  }
#endif
}

/** Setup debugger.
 *
 * \param keyval keyvalue container
 * \param grp message group
 * \param debugger debugger structure
 * \param options parsed command line options
 */
void setupDebugger(
    Ref<KeyVal> &keyval,
    Ref<MessageGrp> &grp,
    Ref<Debugger> &debugger,
    struct OptionValues &values)
{
  debugger << keyval->describedclassvalue("debug");
  if (debugger.nonnull()) {
    Debugger::set_default_debugger(debugger);
    debugger->set_exec(values.executablename.c_str());
    debugger->set_prefix(grp->me());
    if (values.debug)
      debugger->debug("Starting debugger because -d given on command line.");
  }
}

/** Get integral factory.
 *
 * \param keyval keyvalue container
 * \param integral integral group on return
 * \param argc argument count
 * \param argv argument array
 */
void getIntegralFactory(
    Ref<KeyVal> &keyval,
    Ref<Integral> &integral,
    int argc,
    char **argv)
{
  // first try commandline and environment
  integral = Integral::initial_integral(argc, argv);

  // if we still don't have a integral, try reading the integral
  // from the input
  if (integral.null()) {
    integral << keyval->describedclassvalue("integrals");
  }

  if (integral.nonnull())
    Integral::set_default_integral(integral);
  else
    integral = Integral::get_default_integral();

}

void performRestart(
    Ref<KeyVal> &keyval,
    Ref<MessageGrp> &grp,
    Ref<Optimize> &opt,
    Ref<MolecularEnergy> &mole,
    char *&restartfile
    )
{
  int restart = keyval->booleanvalue("restart",truevalue);
  struct stat sb;
  int statresult, statsize;
  if (restart) {
    if (grp->me() == 0) {
      statresult = stat(restartfile,&sb);
      statsize = (statresult==0) ? sb.st_size : 0;
    }
    grp->bcast(statresult);
    grp->bcast(statsize);
  }
  if (restart && statresult==0 && statsize) {
    BcastStateInBin si(grp,restartfile);
    if (keyval->exists("override")) {
      si.set_override(new PrefixKeyVal(keyval,"override"));
    }
    char *suf = strrchr(restartfile,'.');
    if (!strcmp(suf,".wfn")) {
      mole << SavableState::key_restore_state(si,"mole");
      ExEnv::out0() << endl
                   << indent << "Restored <" << mole->class_name()
                   << "> from " << restartfile << endl;

      opt << keyval->describedclassvalue("opt");
      if (opt.nonnull())
        opt->set_function(mole.pointer());
    }
    else {
      opt << SavableState::key_restore_state(si,"opt");
      if (opt.nonnull()) {
        mole << opt->function();
        ExEnv::out0() << endl << indent
             << "Restored <Optimize> from " << restartfile << endl;
      }
    }
  } else {
    mole << keyval->describedclassvalue("mole");
    opt << keyval->describedclassvalue("opt");
  }
}

char *setMolecularCheckpointFile(
    Ref<KeyVal> &keyval,
    Ref<MessageGrp> &grp,
    Ref<MolecularEnergy> &mole,
    char *mole_ckpt_file
    )
{
  int checkpoint = keyval->booleanvalue("checkpoint",truevalue);
  int checkpoint_freq = keyval->intvalue("checkpoint_freq",KeyValValueint(1));
  if (mole.nonnull()) {
    MolecularFormula mf(mole->molecule());
    ExEnv::out0() << endl << indent
         << "Molecular formula " << mf.formula() << endl;
    if (checkpoint) {
      mole->set_checkpoint();
      if (grp->me() == 0) mole->set_checkpoint_file(mole_ckpt_file);
      else mole->set_checkpoint_file(devnull);
      mole->set_checkpoint_freq(checkpoint_freq);
    }
  }
}

/** Checks whether limit on command-line exceeds the basis functions.
 *
 * \param mole molecular energy object
 * \param values temporarily storage for (command-line) options
 * \return 0 - not exceeded, 1 - exceeded
 */
int checkBasisSetLimit(
    Ref<MolecularEnergy> &mole,
    struct OptionValues &values
    )
{
  int check = (values.check != (const char *)0);
  int limit = values.limit;
  if (limit) {
    Ref<Wavefunction> wfn; wfn << mole;
    if (wfn.nonnull() && wfn->ao_dimension()->n() > limit) {
      ExEnv::out0() << endl << indent
           << "The limit of " << limit << " basis functions has been exceeded."
           << endl;
      check = 1;
    }
  }
  return check;
}

/** Performs the energy optimization.
 *
 * \param opt optimization object
 * \param mole molecular energy object
 * \return 0 - not read for frequency calculation, 1 - ready
 */
int performEnergyOptimization(
    Ref<Optimize> &opt,
    Ref<MolecularEnergy> &mole
    )
{
  int ready_for_freq = 0;
  int result = opt->optimize();
  if (result) {
    ExEnv::out0() << indent
         << "The optimization has converged." << endl << endl;
    ExEnv::out0() << indent
         << scprintf("Value of the MolecularEnergy: %15.10f",
                     mole->energy())
         << endl << endl;
    ready_for_freq = 1;
  } else {
    ExEnv::out0() << indent
         << "The optimization has NOT converged." << endl << endl;
    ready_for_freq = 0;
  }
  return ready_for_freq;
}

/** Performs gradient calculation.
 *
 * \param mole molecular energy object
 */
void performGradientCalculation(
    Ref<MolecularEnergy> &mole
    )
{
  mole->do_gradient(1);
  ExEnv::out0() << endl << indent
       << scprintf("Value of the MolecularEnergy: %15.10f",
                   mole->energy())
       << endl;
  if (mole->value_result().actual_accuracy()
      > mole->value_result().desired_accuracy()) {
    ExEnv::out0() << indent
         << "WARNING: desired accuracy not achieved in energy" << endl;
  }
  ExEnv::out0() << endl;
  // Use result_noupdate since the energy might not have converged
  // to the desired accuracy in which case grabbing the result will
  // start up the calculation again.  However the gradient might
  // not have been computed (if we are restarting and the gradient
  // isn't in the save file for example).
  RefSCVector grad;
  if (mole->gradient_result().computed()) {
    grad = mole->gradient_result().result_noupdate();
  }
  else {
    grad = mole->gradient();
  }
  if (grad.nonnull()) {
    grad.print("Gradient of the MolecularEnergy:");
    if (mole->gradient_result().actual_accuracy()
        > mole->gradient_result().desired_accuracy()) {
      ExEnv::out0() << indent
           << "WARNING: desired accuracy not achieved in gradient" << endl;
    }
  }
}

/** Performs frequency calculation.
 *
 * \param mole molecular energy object
 * \param molhess molecular hessian object
 * \param molfreq molecular frequency object
 */
void performFrequencyCalculation(
    Ref<MolecularEnergy> &mole,
    Ref<MolecularHessian> &molhess,
    Ref<MolecularFrequencies> &molfreq

    )
{
  RefSymmSCMatrix xhessian;
  if (molhess.nonnull()) {
    // if "hess" input was given, use it to compute the hessian
    xhessian = molhess->cartesian_hessian();
  }
  else if (mole->hessian_implemented()) {
    // if mole can compute the hessian, use that hessian
    xhessian = mole->get_cartesian_hessian();
  }
  else if (mole->gradient_implemented()) {
    // if mole can compute gradients, use gradients at finite
    // displacements to compute the hessian
    molhess = new FinDispMolecularHessian(mole);
    xhessian = molhess->cartesian_hessian();
  }
  else {
    ExEnv::out0() << "mpqc: WARNING: Frequencies cannot be computed" << endl;
  }

  if (xhessian.nonnull()) {
    char *hessfile = SCFormIO::fileext_to_filename(".hess");
    MolecularHessian::write_cartesian_hessian(hessfile,
                                              mole->molecule(), xhessian);
    delete[] hessfile;

    molfreq->compute_frequencies(xhessian);
    // DEGENERACY IS NOT CORRECT FOR NON-SINGLET CASES:
    molfreq->thermochemistry(1);
  }
}

/** Renders some objects.
 *
 * \param renderer renderer object
 * \param keyval keyvalue container
 * \param tim timing object
 * \param grp message group
 */
void renderObjects(
    Ref<Render> &renderer,
    Ref<KeyVal> &keyval,
    Ref<RegionTimer> &tim,
    Ref<MessageGrp> &grp
    )
{
  Ref<RenderedObject> rendered;
  rendered << keyval->describedclassvalue("rendered");
  Ref<AnimatedObject> animated;
  animated << keyval->describedclassvalue("rendered");
  if (rendered.nonnull()) {
    if (tim.nonnull()) tim->enter("render");
    if (grp->me() == 0) renderer->render(rendered);
    if (tim.nonnull()) tim->exit("render");
  }
  else if (animated.nonnull()) {
    if (tim.nonnull()) tim->enter("render");
    if (grp->me() == 0) renderer->animate(animated);
    if (tim.nonnull()) tim->exit("render");
  }
  else {
    if (tim.nonnull()) tim->enter("render");
    int n = keyval->count("rendered");
    for (int i=0; i<n; i++) {
      rendered << keyval->describedclassvalue("rendered",i);
      animated << keyval->describedclassvalue("rendered",i);
      if (rendered.nonnull()) {
        // make sure the object has a name so we don't overwrite its file
        if (rendered->name() == 0) {
          char ic[64];
          sprintf(ic,"%02d",i);
          rendered->set_name(ic);
        }
        if (grp->me() == 0) renderer->render(rendered);
      }
      else if (animated.nonnull()) {
        // make sure the object has a name so we don't overwrite its file
        if (animated->name() == 0) {
          char ic[64];
          sprintf(ic,"%02d",i);
          animated->set_name(ic);
        }
        if (grp->me() == 0) renderer->animate(animated);
      }
    }
    if (tim.nonnull()) tim->exit("render");
  }
}

/** Save the molecule to PDB file.
 *
 * \param do_pdb whether to save as pdb (1) or not (0)
 * \param grp message group
 * \param mole molecular energy object
 * \param molname name of output file
 */
void saveToPdb(
    int do_pdb,
    Ref<MessageGrp> &grp,
    Ref<MolecularEnergy> &mole,
    const char *molname
    )
{
  if (do_pdb && grp->me() == 0) {
    char *ckptfile = new char[strlen(molname)+5];
    sprintf(ckptfile, "%s.pdb", molname);
    ofstream pdbfile(ckptfile);
    mole->molecule()->print_pdb(pdbfile);
    delete[] ckptfile;
  }
}

static int getCoreElectrons(const int z)
{
  int n=0;
  if (z > 2) n += 2;
  if (z > 10) n += 8;
  if (z > 18) n += 8;
  if (z > 30) n += 10;
  if (z > 36) n += 8;
  if (z > 48) n += 10;
  if (z > 54) n += 8;
  return n;
}

void init()
{
  //trash_stack();

  int i;
  atexit(clean_up);

#ifdef HAVE_FEENABLEEXCEPT
  // this uses a glibc extension to trap on individual exceptions
# ifdef FE_DIVBYZERO
  feenableexcept(FE_DIVBYZERO);
# endif
# ifdef FE_INVALID
  feenableexcept(FE_INVALID);
# endif
# ifdef FE_OVERFLOW
  feenableexcept(FE_OVERFLOW);
# endif
#endif

#ifdef HAVE_FEDISABLEEXCEPT
  // this uses a glibc extension to not trap on individual exceptions
# ifdef FE_UNDERFLOW
  fedisableexcept(FE_UNDERFLOW);
# endif
# ifdef FE_INEXACT
  fedisableexcept(FE_INEXACT);
# endif
#endif

#if defined(HAVE_SETRLIMIT)
  struct rlimit rlim;
  rlim.rlim_cur = 0;
  rlim.rlim_max = 0;
  setrlimit(RLIMIT_CORE,&rlim);
#endif
}

/** Finds the region index to a given timer region name.
 *
 * @param nregion number of regions
 * @param region_names array with name of each region
 * @param name name of desired region
 * @return index of desired region in array
 */
int findTimerRegion(const int &nregion, const char **&region_names, const char *name)
{
  size_t region=0;
  for (;region<nregion;++region) {
    //std::cout << "Comparing " << region_names[region] << " and " << name << "." << std::endl;
    if (strcmp(region_names[region], name) == 0)
      break;
  }
  if (region == nregion)
    region = 0;
  return region;
}

/** Performs the main work to calculate the ground state energies, gradients, etc.
 *
 * @param grp message group
 * @param values temporary value storage for parsed command-line
 * @param input input file name
 * @param generic_input generic input file name
 * @param in_char_array either NULL or contains char array with read input
 * @param argc argument count
 * @param argv argument array
 */
void mainFunction(
    Ref<MessageGrp> grp,
    struct OptionValues &values,
    const char *&output,
    const char *&input,
    const char *&generic_input,
    char *&in_char_array,
    int argc,
    char **argv
#ifdef HAVE_MPQCDATA
    , MPQCData &data
#endif
    )
{
  // get the basename for output files
  setOutputBaseName(input, output);

  // parse input into keyvalue container
  Ref<ParsedKeyVal> parsedkv;
  int use_simple_input = 0; // default is object-oriented if in_char_array is given
  if (!in_char_array) // obtain from file
    parseInputfile(grp, parsedkv, values, input, generic_input, in_char_array, use_simple_input);
  parseIntoKeyValue(parsedkv, values, in_char_array, use_simple_input);
  delete[] in_char_array;

  // prefix parsed values wit "mpqc"
  if (values.keyvalue) parsedkv->verbose(1);
  Ref<KeyVal> keyval = new PrefixKeyVal(parsedkv.pointer(),"mpqc");

  // set up output classes
  setupSCFormIO(grp);

  // initialize timing for mpqc
  Ref<RegionTimer> tim;
  initTimings(grp, keyval, tim);
  
  // announce ourselves
  makeAnnouncement(tim);

  // get the thread group.
  Ref<ThreadGrp> thread;
  getThreadGroup(keyval, thread, argc, argv);

  // get the memory group.
  Ref<MemoryGrp> memory;
  getMemoryGroup(keyval, memory, argc, argv);

  ExEnv::out0() << indent
       << "Using " << grp->class_name()
       << " for message passing (number of nodes = " << grp->n() << ")." << endl
       << indent
       << "Using " << thread->class_name()
       << " for threading (number of threads = " << thread->nthread() << ")." << endl
       << indent
       << "Using " << memory->class_name()
       << " for distributed shared memory." << endl
       << indent
       << "Total number of processors = " << grp->n() * thread->nthread() << endl;

  // prepare CCA if available
  prepareCCA(keyval, values);

  // now set up the debugger
  Ref<Debugger> debugger;
  setupDebugger(keyval, grp, debugger, values);

  // now check to see what matrix kit to use
  if (keyval->exists("matrixkit"))
    SCMatrixKit::set_default_matrixkit(
      dynamic_cast<SCMatrixKit*>(
        keyval->describedclassvalue("matrixkit").pointer()));

  // get the integral factory.
  Ref<Integral> integral;
  getIntegralFactory(keyval, integral, argc, argv);
  ExEnv::out0() << endl << indent
       << "Using " << integral->class_name()
       << " by default for molecular integrals evaluation" << endl << endl;
  
  // create some filenames for molecule, checkpoint, basename of output
  const char *basename = SCFormIO::default_basename();
  KeyValValueString molnamedef(basename);
  char * molname = keyval->pcharvalue("filename", molnamedef);
  if (strcmp(molname, basename))
    SCFormIO::set_default_basename(molname);

  char * ckptfile = new char[strlen(molname)+6];
  sprintf(ckptfile,"%s.ckpt",molname);

  KeyValValueString restartfiledef(ckptfile);
  char * restartfile = keyval->pcharvalue("restart_file", restartfiledef);
  
  char * wfn_file = keyval->pcharvalue("wfn_file");
  if (wfn_file == 0) {
    wfn_file = new char[strlen(molname)+6];
    sprintf(wfn_file,"%s.wfn",molname);
  }
  char *mole_ckpt_file = new char[strlen(wfn_file)+1];
  sprintf(mole_ckpt_file,"%s",wfn_file);

  int savestate = keyval->booleanvalue("savestate",truevalue);

  // setup molecular energy and optimization instances
  Ref<MolecularEnergy> mole;
  Ref<Optimize> opt;

  // read in restart file if we do restart
  performRestart(keyval, grp, opt, mole, restartfile);

  // setup molecule checkpoint file
  setMolecularCheckpointFile(keyval, grp, mole, mole_ckpt_file);
  delete[] mole_ckpt_file;

  int checkpoint = keyval->booleanvalue("checkpoint",truevalue);
  if (checkpoint && opt.nonnull()) {
    opt->set_checkpoint();
    if (grp->me() == 0) opt->set_checkpoint_file(ckptfile);
    else opt->set_checkpoint_file(devnull);
  }

  // see if frequencies are wanted
  Ref<MolecularHessian> molhess;
  molhess << keyval->describedclassvalue("hess");
  Ref<MolecularFrequencies> molfreq;
  molfreq << keyval->describedclassvalue("freq");

  // check basis set limit
  const int check = checkBasisSetLimit(mole, values);
  if (check) {
    ExEnv::out0() << endl << indent
         << "Exiting since the check option is on." << endl;
    exit(0);
  }
  
  // from now on we time the calculations
  if (tim.nonnull()) tim->change("calc");

  int do_energy = keyval->booleanvalue("do_energy",truevalue);

  int do_grad = keyval->booleanvalue("do_gradient",falsevalue);

  int do_opt = keyval->booleanvalue("optimize",truevalue);

  int do_pdb = keyval->booleanvalue("write_pdb",falsevalue);

  int print_mole = keyval->booleanvalue("print_mole",truevalue);

  int print_timings = keyval->booleanvalue("print_timings",truevalue);

  // print all current options (keyvalues)
  printOptions(keyval, opt, molname, restartfile);

  // see if any pictures are desired
  Ref<Render> renderer;
  renderer << keyval->describedclassvalue("renderer");

  // If we have a renderer, then we will read in some more info
  // below.  Otherwise we can get rid of the keyval's, to eliminate
  // superfluous references to objects that we might otherwise be
  // able to delete.  We cannot read in the remaining rendering
  // objects now, since some of their KeyVal CTOR's are heavyweight,
  // requiring optimized geometries, etc.
  if (renderer.null()) {
    if (parsedkv.nonnull()) print_unseen(parsedkv, input);
    keyval = 0;
    parsedkv = 0;
  }

  delete[] restartfile;
  delete[] ckptfile;

  int ready_for_freq = 1;
  if (mole.nonnull()) {
    if (((do_opt && opt.nonnull()) || do_grad)
        && !mole->gradient_implemented()) {
      ExEnv::out0() << indent
           << "WARNING: optimization or gradient requested but the given"
           << endl
           << "         MolecularEnergy object cannot do gradients."
           << endl;
    }

    if (do_opt && opt.nonnull() && mole->gradient_implemented()) {

      ready_for_freq = performEnergyOptimization(opt, mole);

    } else if (do_grad && mole->gradient_implemented()) {

      performGradientCalculation(mole);

    } else if (do_energy && mole->value_implemented()) {
      ExEnv::out0() << endl << indent
           << scprintf("Value of the MolecularEnergy: %15.10f",
                       mole->energy())
           << endl << endl;
    }
  }

  // stop timing of calculations
  if (tim.nonnull()) tim->exit("calc");

  // save this before doing the frequency stuff since that obsoletes the
  saveState(wfn_file, savestate, opt, grp, mole, molname, ckptfile);
  
  // Frequency calculation.
  if (ready_for_freq && molfreq.nonnull()) {
    performFrequencyCalculation(mole, molhess, molfreq);
  }

  if (renderer.nonnull()) {
    renderObjects(renderer, keyval, tim, grp);

    Ref<MolFreqAnimate> molfreqanim;
    molfreqanim << keyval->describedclassvalue("animate_modes");
    if (ready_for_freq && molfreq.nonnull()
        && molfreqanim.nonnull()) {
      if (tim.nonnull()) tim->enter("render");
      molfreq->animate(renderer, molfreqanim);
      if (tim.nonnull()) tim->exit("render");
    }
  }

  if (mole.nonnull()) {
    if (print_mole)
      mole->print(ExEnv::out0());

    saveToPdb(do_pdb, grp, mole, molname);

  }
  else {
    ExEnv::out0() << "mpqc: The molecular energy object is null" << endl
         << " make sure \"mole\" specifies a MolecularEnergy derivative"
         << endl;
  }
  if (parsedkv.nonnull()) print_unseen(parsedkv, input);

  // here, we may gather the results
  // we start to fill the MPQC_Data object
  if (tim.nonnull()) tim->enter("gather");
  {
     Ref<Wavefunction> wfn;
     wfn << mole;
//     ExEnv::out0() << "The number of atomic orbitals: " << wfn->ao_dimension()->n() << endl;
//     ExEnv::out0() << "The AO density matrix is ";
//     wfn->ao_density()->print(ExEnv::out0());
//     ExEnv::out0() << "The natural density matrix is ";
//     wfn->natural_density()->print(ExEnv::out0());
//     ExEnv::out0() << "The Gaussian basis is " << wfn->basis()->name() << endl;
//     ExEnv::out0() << "The Gaussians sit at the following centers: " << endl;
//     for (int nr = 0; nr< wfn->basis()->ncenter(); ++nr) {
//       ExEnv::out0() << nr << " basis function has its center at ";
//       for (int i=0; i < 3; ++i)
//           ExEnv::out0() << wfn->basis()->r(nr,i) << "\t";
//       ExEnv::out0() << endl;
//     }
     // store accuracies
     data.accuracy = mole->value_result().actual_accuracy();
     data.desired_accuracy = mole->value_result().desired_accuracy();
     // print the energy
     data.energies.total = wfn->energy();
     data.energies.nuclear_repulsion = wfn->nuclear_repulsion_energy();
     {
       CLHF *clhf = dynamic_cast<CLHF*>(wfn.pointer());
       if (clhf != NULL) {
         double ex, ec;
         clhf->two_body_energy(ec, ex);
         data.energies.electron_coulomb = ec;
         data.energies.electron_exchange = ex;
         clhf = NULL;
       } else {
         ExEnv::out0() << "INFO: There is no direct CLHF information available." << endl;
         data.energies.electron_coulomb = 0.;
         data.energies.electron_exchange = 0.;
       }
     }
     SCF *scf = NULL;
     {
       MBPT2 *mbpt2 = dynamic_cast<MBPT2*>(wfn.pointer());
       if (mbpt2 != NULL) {
         data.energies.correlation = mbpt2->corr_energy();
         scf = mbpt2->ref().pointer();
         CLHF *clhf = dynamic_cast<CLHF*>(scf);
         if (clhf != NULL) {
           double ex, ec;
           clhf->two_body_energy(ec, ex);
           data.energies.electron_coulomb = ec;
           data.energies.electron_exchange = ex;
           clhf = NULL;
         } else {
           ExEnv::out0() << "INFO: There is no reference CLHF information available either." << endl;
           data.energies.electron_coulomb = 0.;
           data.energies.electron_exchange = 0.;
         }
         mbpt2 = 0;
       } else {
         ExEnv::out0() << "INFO: There is no MBPT2 information available." << endl;
         data.energies.correlation = 0.;
         scf = dynamic_cast<SCF*>(wfn.pointer());
         if (scf == NULL)
           abort();
       }
     }
     {
       // taken from clscf.cc: CLSCF::scf_energy() (but see also Szabo/Ostlund)

       RefSymmSCMatrix t = scf->overlap();
       RefSymmSCMatrix cl_dens_ = scf->ao_density();

       SCFEnergy *eop = new SCFEnergy;
       eop->reference();
       Ref<SCElementOp2> op = eop;
       t.element_op(op,cl_dens_);
       op=0;
       eop->dereference();

       data.energies.overlap = eop->result();

       delete eop;
       t = 0;
       cl_dens_ = 0;
     }
     {
       // taken from Wavefunction::core_hamiltonian()
       RefSymmSCMatrix hao(scf->basis()->basisdim(), scf->basis()->matrixkit());
       hao.assign(0.0);
       Ref<PetiteList> pl = scf->integral()->petite_list();
       Ref<SCElementOp> hc =
         new OneBodyIntOp(new SymmOneBodyIntIter(scf->integral()->kinetic(), pl));
       hao.element_op(hc);
       hc=0;

       RefSymmSCMatrix h(scf->so_dimension(), scf->basis_matrixkit());
       pl->symmetrize(hao,h);

       // taken from clscf.cc: CLSCF::scf_energy() (but see also Szabo/Ostlund)
       RefSymmSCMatrix cl_dens_ = scf->ao_density();

       SCFEnergy *eop = new SCFEnergy;
       eop->reference();
       Ref<SCElementOp2> op = eop;
       h.element_op(op,cl_dens_);
       op=0;
       eop->dereference();

       data.energies.kinetic = 2.*eop->result();

       delete eop;
       hao = 0;
       h = 0;
       cl_dens_ = 0;
     }
     {
       // set to potential energy between nuclei and electron charge distribution
       RefSymmSCMatrix hao(scf->basis()->basisdim(), scf->basis()->matrixkit());
       hao.assign(0.0);
       Ref<PetiteList> pl = scf->integral()->petite_list();
       Ref<SCElementOp> hc =
         new OneBodyIntOp(new SymmOneBodyIntIter(scf->integral()->nuclear(), pl));
       hao.element_op(hc);
       hc=0;

       RefSymmSCMatrix h(scf->so_dimension(), scf->basis_matrixkit());
       pl->symmetrize(hao,h);

       // taken from clscf.cc: CLSCF::scf_energy() (but see also Szabo/Ostlund)
       RefSymmSCMatrix cl_dens_ = scf->ao_density();

       SCFEnergy *eop = new SCFEnergy;
       eop->reference();
       Ref<SCElementOp2> op = eop;
       h.element_op(op,cl_dens_);
       op=0;
       eop->dereference();

       data.energies.hcore = 2.*eop->result();

       delete eop;
       hao = 0;
       h = 0;
       cl_dens_ = 0;
     }
     ExEnv::out0() << "total is " << data.energies.total << endl;
     ExEnv::out0() << "nuclear_repulsion is " << data.energies.nuclear_repulsion << endl;
     ExEnv::out0() << "electron_coulomb is " << data.energies.electron_coulomb << endl;
     ExEnv::out0() << "electron_exchange is " << data.energies.electron_exchange << endl;
     ExEnv::out0() << "correlation is " << data.energies.correlation << endl;
     ExEnv::out0() << "overlap is " << data.energies.overlap << endl;
     ExEnv::out0() << "kinetic is " << data.energies.kinetic << endl;
     ExEnv::out0() << "hcore is " << data.energies.hcore << endl;
     ExEnv::out0() << "sum is " <<
         data.energies.nuclear_repulsion
         + data.energies.electron_coulomb
         + data.energies.electron_exchange
         + data.energies.correlation
         + data.energies.kinetic
         + data.energies.hcore
         << endl;

     ExEnv::out0() << endl << indent
          << scprintf("Value of the MolecularEnergy: %15.10f",
                      mole->energy())
          << endl;
     // print the gradient
     RefSCVector grad;
     if (mole->gradient_result().computed()) {
       grad = mole->gradient_result().result_noupdate();
     }
     else {
       grad = mole->gradient();
     }
     if (grad.nonnull()) {
       data.forces.resize(grad.dim()/3);
       for (int j=0;j<grad.dim()/3; ++j) {
         data.forces[j].resize(3, 0.);
       }
       ExEnv::out0() << "Gradient of the MolecularEnergy:" << std::endl;
       for (int j=0;j<grad.dim()/3; ++j) {
         ExEnv::out0() << "\t";
         for (int i=0; i< 3; ++i) {
           data.forces[j][i] = grad[3*j+i];
           ExEnv::out0() << grad[3*j+i] << "\t";
         }
         ExEnv::out0() << endl;
       }
     }
     grad = NULL;

     {
       // eigenvalues (this only works if we have a OneBodyWavefunction, i.e. SCF procedure)
//       SCF *scf = dynamic_cast<SCF*>(wfn.pointer());
//       if (scf != NULL) {
//         const double scfernergy = scf->energy();
         RefDiagSCMatrix evals = scf->eigenvalues();

         ExEnv::out0() << "Eigenvalues:" << endl;
         for(int i=0;i<wfn->oso_dimension(); ++i) {
           data.energies.eigenvalues.push_back(evals(i));
           ExEnv::out0() << i << "th eigenvalue is " << evals(i) << endl;
         }
//       }
     }
     // we do sample the density only on request
       {
         // fill positions and charges (NO LONGER converting from bohr radii to angstroem)
         const double AtomicLengthToAngstroem = 1.;//0.52917721;
         data.positions.reserve(wfn->molecule()->natom());
         data.atomicnumbers.reserve(wfn->molecule()->natom());
         data.charges.reserve(wfn->molecule()->natom());
         for (int iatom=0;iatom < wfn->molecule()->natom(); ++iatom) {
           data.atomicnumbers.push_back(wfn->molecule()->Z(iatom));
           double charge = wfn->molecule()->Z(iatom);
           if (data.DoValenceOnly == MPQCData::DoSampleValenceOnly)
             charge -= getCoreElectrons((int)charge);
           data.charges.push_back(charge);
           std::vector<double> pos(3, 0.);
           for (int j=0;j<3;++j)
             pos[j] = wfn->molecule()->r(iatom, j)*AtomicLengthToAngstroem;
           data.positions.push_back(pos);
         }
         ExEnv::out0() << "We have "
             << data.positions.size() << " positions and "
             << data.charges.size() << " charges." << endl;
       }
       if (data.DoLongrange) {
       if (data.sampled_grid.level != 0)
       {
         // we now need to sample the density on the grid
         // 1. get max and min over all basis function positions
         assert( scf->basis()->ncenter() > 0 );
         SCVector3 bmin( scf->basis()->r(0,0), scf->basis()->r(0,1), scf->basis()->r(0,2) );
         SCVector3 bmax( scf->basis()->r(0,0), scf->basis()->r(0,1), scf->basis()->r(0,2) );
         for (int nr = 1; nr< scf->basis()->ncenter(); ++nr) {
           for (int i=0; i < 3; ++i) {
             if (scf->basis()->r(nr,i) < bmin(i))
               bmin(i) = scf->basis()->r(nr,i);
             if (scf->basis()->r(nr,i) > bmax(i))
               bmax(i) = scf->basis()->r(nr,i);
           }
         }
         ExEnv::out0() << "Basis min is at " << bmin << " and max is at " << bmax << endl;

         // 2. choose an appropriately large grid
         // we have to pay attention to capture the right amount of the exponential decay
         // and also to have a power of two size of the grid at best
         SCVector3 boundaryV(5.);  // boundary extent around compact domain containing basis functions
         bmin -= boundaryV;
         bmax += boundaryV;
         for (size_t i=0;i<3;++i) {
           if (bmin(i) < data.sampled_grid.begin[i])
             bmin(i) = data.sampled_grid.begin[i];
           if (bmax(i) > data.sampled_grid.end[i])
             bmax(i) = data.sampled_grid.end[i];
         }
         // set the non-zero window of the sampled_grid
         data.sampled_grid.setWindow(bmin.data(), bmax.data());

         // for the moment we always generate a grid of full size
         // (NO LONGER converting grid dimensions from angstroem to bohr radii)
         const double AtomicLengthToAngstroem = 1.;//0.52917721;
         SCVector3 min;
         SCVector3 max;
         SCVector3 delta;
         size_t samplepoints[3];
         // due to periodic boundary conditions, we don't need gridpoints-1 here
         // TODO: in case of open boundary conditions, we need data on the right
         // hand side boundary as well
         const int gridpoints = data.sampled_grid.getGridPointsPerAxis();
         for (size_t i=0;i<3;++i) {
           min(i) = data.sampled_grid.begin_window[i]/AtomicLengthToAngstroem;
           max(i) = data.sampled_grid.end_window[i]/AtomicLengthToAngstroem;
           delta(i) = data.sampled_grid.getDeltaPerAxis(i)/AtomicLengthToAngstroem;
           samplepoints[i] = data.sampled_grid.getWindowGridPointsPerAxis(i);
         }
         ExEnv::out0() << "Grid starts at " << min
             << " and ends at " << max
             << " with a delta of " << delta
             << " to get "
             << samplepoints[0] << ","
             << samplepoints[1] << ","
             << samplepoints[2] << " samplepoints."
             << endl;
         assert( data.sampled_grid.sampled_grid.size() == samplepoints[0]*samplepoints[1]*samplepoints[2]);

         // 3. sample the atomic density
         const double element_volume_conversion =
             1./AtomicLengthToAngstroem/AtomicLengthToAngstroem/AtomicLengthToAngstroem;
         SCVector3 r = min;

         std::set<int> valence_indices;
         RefDiagSCMatrix evals = scf->eigenvalues();
         if (data.DoValenceOnly == MPQCData::DoSampleValenceOnly) {
           // find valence orbitals
//           std::cout << "All Eigenvalues:" << std::endl;
//           for(int i=0;i<wfn->oso_dimension(); ++i)
//             std::cout << i << "th eigenvalue is " << evals(i) << std::endl;
           int n_electrons = scf->nelectron();
           int n_core_electrons =  wfn->molecule()->n_core_electrons();
           std::set<double> evals_sorted;
           {
             int i=0;
             double first_positive_ev = std::numeric_limits<double>::max();
             for(i=0;i<wfn->oso_dimension(); ++i) {
               if (evals(i) < 0.)
                 evals_sorted.insert(evals(i));
               else
                 first_positive_ev = std::min(first_positive_ev, (double)evals(i));
             }
             // add the first positive for the distance
             evals_sorted.insert(first_positive_ev);
           }
           std::set<double> evals_distances;
           std::set<double>::const_iterator advancer = evals_sorted.begin();
           std::set<double>::const_iterator iter = advancer++;
           for(;advancer != evals_sorted.end(); ++advancer,++iter)
             evals_distances.insert((*advancer)-(*iter));
           const double largest_distance = *(evals_distances.rbegin());
           ExEnv::out0()  << "Largest distance between EV is " << largest_distance << std::endl;
           advancer = evals_sorted.begin();
           iter = advancer++;
           for(;advancer != evals_sorted.begin(); ++advancer,++iter)
             if (fabs(fabs((*advancer)-(*iter)) - largest_distance) < 1e-10)
               break;
           assert( advancer != evals_sorted.begin() );
           const double last_core_ev = (*iter);
           ExEnv::out0()  << "Last core EV might be " << last_core_ev << std::endl;
           ExEnv::out0()  << "First valence index is " << n_core_electrons/2 << std::endl;
           for(int i=n_core_electrons/2;i<wfn->oso_dimension(); ++i)
             if (evals(i) > last_core_ev)
               valence_indices.insert(i);
//           {
//             int i=0;
//             std::cout << "Valence eigenvalues:" << std::endl;
//             for (std::set<int>::const_iterator iter = valence_indices.begin();
//                 iter != valence_indices.end(); ++iter)
//               std::cout << i++ << "th eigenvalue is " << (*iter) << std::endl;
//           }
         } else {
           // just insert all indices
           for(int i=0;i<wfn->oso_dimension(); ++i)
             valence_indices.insert(i);
         }

         // testing alternative routine from SCF::so_density()
         RefSCMatrix oso_vector = scf->oso_eigenvectors();
         RefSCMatrix vector = scf->so_to_orthog_so().t() * oso_vector;
         oso_vector = 0;
         RefSymmSCMatrix occ(scf->oso_dimension(), scf->basis_matrixkit());
         occ.assign(0.0);
         for (std::set<int>::const_iterator iter = valence_indices.begin();
             iter != valence_indices.end(); ++iter) {
           const int i = *iter;
           occ(i,i) = scf->occupation(i);
           ExEnv::out0() << "# " << i << " has ev of " << evals(i) << ", occupied by " << scf->occupation(i) << std::endl;
         }
         RefSymmSCMatrix d2(scf->so_dimension(), scf->basis_matrixkit());
         d2.assign(0.0);
         d2.accumulate_transform(vector, occ);

         // taken from scf::density()
         RefSCMatrix nos
             = scf->integral()->petite_list()->evecs_to_AO_basis(scf->natural_orbitals());
         RefDiagSCMatrix nd = scf->natural_density();
         GaussianBasisSet::ValueData *valdat
           = new GaussianBasisSet::ValueData(scf->basis(), scf->integral());
         std::vector<double>::iterator griditer = data.sampled_grid.sampled_grid.begin();
         const int nbasis = scf->basis()->nbasis();
         double *bs_values = new double[nbasis];

         // TODO: need to take care when we have periodic boundary conditions.
         for (size_t x = 0; x < samplepoints[0]; ++x, r.x() += delta(0)) {
           std::cout << "Sampling now for x=" << r.x() << std::endl;
           for (size_t y = 0; y < samplepoints[1]; ++y, r.y() += delta(1)) {
             for (size_t z = 0; z < samplepoints[2]; ++z, r.z() += delta(2)) {
               scf->basis()->values(r,valdat,bs_values);

               // loop over natural orbitals adding contributions to elec_density
               double elec_density=0.0;
               for (int i=0; i<nbasis; ++i) {
                   double tmp = 0.0;
                   for (int j=0; j<nbasis; ++j) {
                       tmp += d2(j,i)*bs_values[j]*bs_values[i];
                     }
                   elec_density += tmp;
                 }
               const double dens_at_r = elec_density * element_volume_conversion;
  //             const double dens_at_r = scf->density(r) * element_volume_conversion;

  //             if (fabs(dens_at_r) > 1e-4)
  //               std::cout << "Electron density at " << r << " is " << dens_at_r << std::endl;
               if (griditer != data.sampled_grid.sampled_grid.end())
                 *griditer++ = dens_at_r;
               else
                 std::cerr << "PAST RANGE!" << std::endl;
             }
             r.z() = min.z();
           }
           r.y() = min.y();
         }
         delete[] bs_values;
         delete valdat;
         assert( griditer == data.sampled_grid.sampled_grid.end());
         // normalization of electron charge to equal electron number
         {
           double integral_value = 0.;
           const double volume_element = pow(AtomicLengthToAngstroem,3)*delta(0)*delta(1)*delta(2);
           for (std::vector<double>::const_iterator diter = data.sampled_grid.sampled_grid.begin();
               diter != data.sampled_grid.sampled_grid.end(); ++diter)
               integral_value += *diter;
           integral_value *= volume_element;
           int n_electrons = scf->nelectron();
           if (data.DoValenceOnly == MPQCData::DoSampleValenceOnly)
             n_electrons -= wfn->molecule()->n_core_electrons();
           const double normalization =
               ((integral_value == 0) || (n_electrons == 0)) ?
                   1. : n_electrons/integral_value;
           std::cout << "Created " << data.sampled_grid.sampled_grid.size() << " grid points"
               << " with integral value of " << integral_value
               << " against " << ((data.DoValenceOnly == MPQCData::DoSampleValenceOnly) ? "n_valence_electrons" : "n_electrons")
               << " of " << n_electrons << "." << std::endl;
           // with normalization we also get the charge right : -1.
           for (std::vector<double>::iterator diter = data.sampled_grid.sampled_grid.begin();
               diter != data.sampled_grid.sampled_grid.end(); ++diter)
               *diter *= -1.*normalization;
         }
       }
     }
     scf = 0;
  }
  if (tim.nonnull()) tim->exit("gather");

  if (print_timings)
    if (tim.nonnull()) tim->print(ExEnv::out0());


  {
    // times obtain from key "mpqc" which should be the first
    const int nregion = tim->nregion();
    //std::cout << "There are " << nregion << " timed regions." << std::endl;
    const char **region_names = new const char*[nregion];
    tim->get_region_names(region_names);
    // find "gather"
    size_t gather_region = findTimerRegion(nregion, region_names, "gather");
    size_t mpqc_region = findTimerRegion(nregion, region_names, "mpqc");
    delete[] region_names;

    // get timings
    double *cpu_time = new double[nregion];
    double *wall_time = new double[nregion];
    double *flops = new double[nregion];
    tim->get_cpu_times(cpu_time);
    tim->get_wall_times(wall_time);
    tim->get_flops(flops);
    if (cpu_time != NULL) {
      data.times.total_cputime = cpu_time[mpqc_region];
      data.times.gather_cputime = cpu_time[gather_region];
    }
    if (wall_time != NULL) {
      data.times.total_walltime = wall_time[mpqc_region];
      data.times.gather_walltime = wall_time[gather_region];
    }
    if (flops != NULL) {
      data.times.total_flops = flops[mpqc_region];
      data.times.gather_flops = flops[gather_region];
    }
    delete[] cpu_time;
    delete[] wall_time;
    delete[] flops;
  }

  delete[] molname;
  SCFormIO::set_default_basename(0);

  renderer = 0;
  molfreq = 0;
  molhess = 0;
  opt = 0;
  mole = 0;
  integral = 0;
  debugger = 0;
  thread = 0;
  tim = 0;
  keyval = 0;
  parsedkv = 0;
  memory = 0;
  clean_up();

#if defined(HAVE_TIME) && defined(HAVE_CTIME)
  time_t t;
  time(&t);
  const char *tstr = ctime(&t);
#endif
  if (!tstr) {
    tstr = "UNKNOWN";
  }
  ExEnv::out0() << endl
               << indent << scprintf("End Time: %s", tstr) << endl;
}

// static values object
OptionValues values;

#ifdef HAVE_JOBMARKET
FragmentResult::ptr MPQCJob::Work()
{
  char mpqc[] = "mpqc" ;
  char **argv = new char*[1];
  argv[0] = &mpqc[0];
  int argc = 1;
//  init();
//
//  ExEnv::init(argc, argv);
//
//  // parse commandline options
//  GetLongOpt options;
//  int optind = ParseOptions(options, argc, argv);
//  const char *output = 0;
//  ostream *outstream = 0;
//  ComputeOptions(options, output, outstream);
//  OptionValues values;
//  parseRemainderOptions(options, values, argc, argv);
//
//  // get the basename for output files
//  char filename_template[] = "mpqc_temp_XXXXXX";
//  output = mktemp(filename_template);
//  setOutputBaseName(NULL, output);

  // now comes the actual work
  int nfilebase = (int) inputfile.length();
  char *in_char_array = new char[nfilebase + 1];
  strncpy(in_char_array, inputfile.c_str(), nfilebase);
  in_char_array[nfilebase] = '\0';
  const char *input = 0;
  const char *generic_input = 0;
  Ref<MessageGrp> grp = MessageGrp::get_default_messagegrp();
  // create unique, temporary name and check whether it exists
  const char *output = NULL;
  char *tempfilename = NULL;
  {
    std::ifstream test;
    do {
      if (output != NULL) // free buffer from last round
        delete output;
      char filename_template[] = "mpqc_temp_XXXXXX\0";
      char filename_suffix[] = ".in\0";
      tempfilename = (char *) malloc ( (strlen(filename_template)+strlen(filename_suffix)+2)*(sizeof(char)));
      strncpy(tempfilename, mktemp(filename_template), strlen(filename_template));
      tempfilename[strlen(filename_template)] = '\0';
      strncat(tempfilename, filename_suffix, strlen(filename_suffix));
      output = tempfilename;
      //free(tempfilename); // don't free! output takes over pointer!
      test.open(output);
    } while (test.good()); // test whether file does not(!) exist
    test.close();
  }
  // put info how to sample the density into MPQCData
  MPQCData data(grid);
  data.DoLongrange = DoLongrange; // set whether we sample the density
  data.DoValenceOnly = DoValenceOnly; // set whether we sample just valence electron and nuclei densities
// now call work horse
  try {
    mainFunction(grp, values, output, input, generic_input, in_char_array, argc, argv, data);
  }
  catch (SCException &e) {
    cout << argv[0] << ": ERROR: SC EXCEPTION RAISED:" << endl
         << e.what()
         << endl;
    clean_up();
  }

  //delete[] in_char_array; // is deleted in mainFunction()
  if (output != 0) {
    free(tempfilename);
  }
  delete[] argv;
  grp = NULL;

  // place into returnstream
  std::stringstream returnstream;
  boost::archive::text_oarchive oa(returnstream);
  oa << data;

  FragmentResult::ptr s( new FragmentResult(getId(), returnstream.str()) );
  if (s->exitflag != 0)
    cerr << "Job #" << s->getId() << " failed to reach desired accuracy." << endl;

  return s;
}
#endif

// we need to explicitly instantiate the serialization functions as
// its is only serialized through its base class FragmentJob
BOOST_CLASS_EXPORT_IMPLEMENT(MPQCJob)

int
try_main(int argc, char *argv[])
{
  init();

  ExEnv::init(argc, argv);

  // parse commandline options
  GetLongOpt options;
  int optind = ParseOptions(options, argc, argv);
  const char *output = 0;
  ostream *outstream = 0;
  ComputeOptions(options, output, outstream);
  parseRemainderOptions(options, values, argc, argv);

  // get input file names, either object-oriented or generic
  const char *object_input = 0;
  const char *generic_input = 0;
  getInputFileNames(object_input, generic_input, options, optind, argc, argv);
  const char *input = 0;
  if (object_input) input = object_input;
  if (generic_input) input = generic_input;

  // get the message group.  first try the commandline and environment
  Ref<MessageGrp> grp;
  getMessageGroup(grp, argc, argv);

  // check if we got option "-n"
  int exitflag = 0;
#ifdef HAVE_JOBMARKET
  if (options.retrieve("n")) {
    /// create new argc, argv and call by splitting with tokenizer
    std::string networkstring(options.retrieve("n"));
    typedef boost::tokenizer<boost::char_separator<char> > tokenizer;
    boost::char_separator<char> sep(" ");
    tokenizer tok(networkstring, sep);
    // simple count because tokenizer has now size()
    int argc_network = 0;
    for(tokenizer::iterator beg=tok.begin(); beg!=tok.end();++beg)
      ++argc_network;
    // argv[0] is program name
    char **argv_network = new char*[argc_network+1];
    argv_network[0] = new char[5];
    strcpy(argv_network[0], "mpqc");
    // then we place each token as a new argument
    argc_network = 1;
    for(tokenizer::iterator beg=tok.begin(); beg!=tok.end();++beg){
      const size_t strlength = (*beg).length();
      const char *strarray = beg->c_str();
      //std::cout << "Token " << argc_network << " is " << strarray << ", length " << strlength << endl;
      argv_network[argc_network] = new char[strlength+1];
      strcpy(argv_network[argc_network], strarray);
      argv_network[argc_network][strlength] = '\0';
      for (size_t index = 0; index < strlength; ++index)
        if (argv_network[argc_network][index] == '+')
          argv_network[argc_network][index] = '-';
      ++argc_network;
    }
    /// and start listening for MPQCJobs
    exitflag = poolworker_main(argc_network, argv_network);
    /// remove the artifical [argv,argv] again
    for (int i=0;i<argc_network;++i)
      delete[] argv_network[i];
    delete[] argv_network;
  } else
#endif
  {
    // if not, work on the command line input
    char *in_char_array = 0;
#ifdef HAVE_MPQCDATA
    MPQCData data;
    mainFunction(grp, values, output, input, generic_input, in_char_array, argc, argv, data);
#else
    mainFunction(grp, values, output, input, generic_input, in_char_array, argc, argv);
#endif
  }

  if (output != 0) {
    ExEnv::set_out(&cout);
    delete outstream;
  }

  grp = 0;
  final_clean_up();

  return exitflag;
}


double EvaluateDensity(SCVector3 &r, Ref<Integral> &intgrl, GaussianBasisSet::ValueData &vdat, Ref<Wavefunction> &wfn)
{
  ExEnv::out0() << "We get the following values at " << r << "." << endl;
  int nbasis = wfn->basis()->nbasis();
  double *b_val = new double[nbasis];
  wfn->basis()->values(r, &vdat, b_val);
  double sum=0.;
  for (int i=0; i<nbasis; i++)
    sum += b_val[i];
  delete[] b_val;
  return sum;
}

int
main(int argc, char *argv[])
{
  size_t exitflag = 0;
  try {
    exitflag = try_main(argc, argv);
  }
  catch (SCException &e) {
    cout << argv[0] << ": ERROR: SC EXCEPTION RAISED:" << endl
         << e.what()
         << endl;
    clean_up();
    throw;
  }
  catch (bad_alloc &e) {
    cout << argv[0] << ": ERROR: MEMORY ALLOCATION FAILED:" << endl
         << e.what()
         << endl;
    clean_up();
    throw;
  }
  catch (exception &e) {
    cout << argv[0] << ": ERROR: EXCEPTION RAISED:" << endl
         << e.what()
         << endl;
    clean_up();
    throw;
  }
  catch (...) {
    cout << argv[0] << ": ERROR: UNKNOWN EXCEPTION RAISED" << endl;
    clean_up();
    throw;
  }
  return exitflag;
}

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

// Local Variables:
// mode: c++
// c-file-style: "ETS"
// End: