source: src/Jobs/mpqc_extract.cc@ 1ca493a

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