//
// mpqc_extract.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.
//
// \note This was extracted from \file mpqc.cc for refactoring into a library.

#ifdef HAVE_CONFIG_H
#include <scconfig.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 "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

#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

#include "mpqc_extract.h"

using namespace std;
using namespace sc;

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;
}

/** 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)
{
  int 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;
}

/** Extractor function that is called after all calculations have been made.
 *
 * \param data result structure to fill
 */
void extractResults(
	 Ref<MolecularEnergy> &mole,
	 void *_data
	 )
{
  MPQCData &data = *static_cast<MPQCData *>(_data);
	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();
	  if (t.dim()->equiv(cl_dens_.dim())) {
	    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();
	  if (h.dim()->equiv(cl_dens_.dim())) {
	    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();
	  if (h.dim()->equiv(cl_dens_.dim())) {
	    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();
	}
	// gradient calculation needs to be activated in the configuration
	// some methods such as open shell MBPT2 do not allow for gradient calc.
//     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;
	  }
	} else {
	  ExEnv::out0() << "INFO: There is no gradient information available." << endl;
	}
	grad = NULL;

	{
	  // eigenvalues (this only works if we have a OneBodyWavefunction, i.e. SCF procedure)
	  // eigenvalues seem to be invalid for unrestricted SCF calculation
	  // (see UnrestrictedSCF::eigenvalues() implementation)
	  UnrestrictedSCF *uscf = dynamic_cast<UnrestrictedSCF*>(wfn.pointer());
	  if ((scf != NULL) && (uscf == 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;
	    }
	  } else {
	      ExEnv::out0() << "INFO: There is no eigenvalue information available." << 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;
}

void extractTimings(
	 Ref<RegionTimer> &tim,
	 void *_data)
{
  MPQCData &data = *static_cast<MPQCData *>(_data);
  // 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;
}