Changeset 1ed5fc for src/bin/mpqc/mpqc.cc

Timestamp:

Aug 1, 2012, 4:19:45 PM (13 years ago)

Author:

Frederik Heber <heber@…>

Children:

d901f7

Parents:

5ce17f

Message:

FIX: Corrected times and added some so far unset values.

• get_cpu_time() and alikes just gives the cputime, not the time from start. We now use get_cpu_times(double*) andn pick the first element.
• electron exchange and coulomb (two-body) is obtained via CLHF.
• correlation is obtained via MBPT2.
• kinetic is obtained in a similar way to hcore, only we require just one half of the core_hamiltonian().

File:

• src/bin/mpqc/mpqc.cc (modified) (4 diffs)

src/bin/mpqc/mpqc.cc

@@ -139 +139 @@
 #include "Jobs/MPQCJob.hpp"
 #include "Jobs/MPQCData.hpp"
+
+#include <chemistry/qc/basis/obint.h>
+#include <chemistry/qc/basis/symmint.h>
 #endif
 
@@ -1480 +1483 @@
      data.energies.total = mole->energy();
      data.energies.nuclear_repulsion = wfn->nuclear_repulsion_energy();
-//     data.energies.electron_repulsion = integral->electron_repulsion();
-//     data.energies.correlation = wfn->corr_energy();
-     // taken from clscf.cc: CLSCF::scf_energy() (but see also Szabo/Ostlund)
      {
+       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 CLHF information available." << endl;
+         data.energies.electron_coulomb = 0.;
+         data.energies.electron_exchange = 0.;
+       }
+     }
+     {
+       MBPT2 *mbpt2 = dynamic_cast<MBPT2*>(wfn.pointer());
+       if (mbpt2 != NULL) {
+         data.energies.correlation = mbpt2->corr_energy();
+         mbpt2 = 0;
+       } else {
+         ExEnv::out0() << "INFO: There is no MBPT2 information available." << endl;
+         data.energies.correlation = 0.;
+       }
+     }
+     {
+       // taken from clscf.cc: CLSCF::scf_energy() (but see also Szabo/Ostlund)
        RefSymmSCMatrix t = wfn->overlap();
        RefSymmSCMatrix cl_dens_ = wfn->ao_density();
@@ -1500 +1525 @@
        cl_dens_ = 0;
      }
-//     data.energies.kinetic = integral->kinetic();
+     {
+       // taken from Wavefunction::core_hamiltonian()
+       RefSymmSCMatrix hao(wfn->basis()->basisdim(), wfn->basis()->matrixkit());
+       hao.assign(0.0);
+       Ref<PetiteList> pl = wfn->integral()->petite_list();
+       Ref<SCElementOp> hc =
+         new OneBodyIntOp(new SymmOneBodyIntIter(wfn->integral()->kinetic(), pl));
+       hao.element_op(hc);
+       hc=0;
+
+       RefSymmSCMatrix h(wfn->so_dimension(), wfn->basis_matrixkit());
+       pl->symmetrize(hao,h);
+
+       // taken from clscf.cc: CLSCF::scf_energy() (but see also Szabo/Ostlund)
+       RefSymmSCMatrix cl_dens_ = wfn->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 = eop->result();
+
+       delete eop;
+       hao = 0;
+       h = 0;
+       cl_dens_ = 0;
+     }
      {
        RefSymmSCMatrix t = wfn->core_hamiltonian();
@@ -1545 +1599 @@
        }
      }
-     // times
-     data.times.walltime = tim->get_wall_time();
-     data.times.cputime = tim->get_cpu_time();
-     data.times.flops = tim->get_flops();
-
-     // eigenvalues (this only works if we have a OneBodyWavefunction, i.e. SCF procedure)
-     SCF *scf = require_dynamic_cast<SCF*>(wfn, "mainFunction");
-     RefDiagSCMatrix evals = scf->eigenvalues();
-
-     for(int i=0;i<wfn->oso_dimension(); ++i) {
-       data.energies.eigenvalues.push_back(evals(i));
-       std::cout << i << "th eigenvalue is " << evals(i) << std::endl;
+     {
+       // times obtain from key "mpqc" which should be the first
+       double *cpu_time = new double[tim->nregion()];
+       double *wall_time = new double[tim->nregion()];
+       double *flops = new double[tim->nregion()];
+       tim->get_cpu_times(cpu_time);
+       tim->get_wall_times(wall_time);
+       tim->get_flops(flops);
+       if (cpu_time != NULL)
+         data.times.cputime = cpu_time[0];
+       if (wall_time != NULL)
+         data.times.walltime = wall_time[0];
+       if (flops != NULL)
+         data.times.flops = flops[0];
+       delete cpu_time;
+       delete wall_time;
+       delete flops;
      }
 
-     // we now need to sample the density on the grid
-     // 1. get max and min over all basis function positions
-     assert( wfn->basis()->ncenter() > 0 );
-     SCVector3 min( wfn->basis()->r(0,0), wfn->basis()->r(0,1), wfn->basis()->r(0,2) );
-     SCVector3 max( wfn->basis()->r(0,0), wfn->basis()->r(0,1), wfn->basis()->r(0,2) );
-     for (int nr = 1; nr< wfn->basis()->ncenter(); ++nr) {
-       for (int i=0; i < 3; ++i) {
-         if (wfn->basis()->r(nr,i) < min(i))
-           min(i) = wfn->basis()->r(nr,i);
-         if (wfn->basis()->r(nr,i) > max(i))
-           max(i) = wfn->basis()->r(nr,i);
+     {
+       // eigenvalues (this only works if we have a OneBodyWavefunction, i.e. SCF procedure)
+       SCF *scf = dynamic_cast<SCF*>(wfn.pointer());
+       if (scf != NULL) {
+         RefDiagSCMatrix evals = scf->eigenvalues();
+
+         for(int i=0;i<wfn->oso_dimension(); ++i) {
+           data.energies.eigenvalues.push_back(evals(i));
+           //std::cout << i << "th eigenvalue is " << evals(i) << std::endl;
+         }
        }
      }
-     std::cout << "Min is at " << min << " and max is at " << max << std::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
-     double boundary = 5.;  // boundary extent around compact domain containing basis functions
-     double delta = 1.; // step width in density sampling
-
-     // for the moment we always generate a grid of full size
-     min.x() = data.sampled_grid.begin[0];
-     min.y() = data.sampled_grid.begin[1];
-     min.z() = data.sampled_grid.begin[2];
-     max.x() = min.x() + data.sampled_grid.size;
-     max.y() = min.y() + data.sampled_grid.size;
-     max.z() = min.z() + data.sampled_grid.size;
-     const int gridpoints = pow(2,data.sampled_grid.level);
-     delta = data.sampled_grid.size / (double) gridpoints;
-     std::cout << "Grid starts at " << min
-         << " and ends at " << max
-         << " with a delta of " << delta
-         << " to get " << gridpoints << " gridpoints."
-         << std::endl;
-
-     // 3. sample the atomic density
-     int nbasis = wfn->basis()->nbasis();
-     double *b_val = new double[nbasis];
-     Ref<Integral> intgrl = Integral::get_default_integral();
-     GaussianBasisSet::ValueData vdat(wfn->basis(), integral);
-     SCVector3 r;
-     data.sampled_grid.sampled_grid.clear();
-     data.sampled_grid.sampled_grid.reserve(gridpoints*gridpoints*gridpoints);
-     for (r.x() = min.x() ; r.x() < max.x(); r.x() += delta)
-       for (r.y() = min.y() ; r.y() < max.y(); r.y() += delta)
-         for (r.z() =min.z() ; r.z() < max.z(); r.z() += delta) {
-           wfn->basis()->values(r, &vdat, b_val);
-           double sum = 0.;
-           for (int i=0; i<nbasis; i++)
-             sum += wfn->ao_density()->get_element(i,i)*b_val[i];
-           data.sampled_grid.sampled_grid.push_back(sum);
+     {
+       // we now need to sample the density on the grid
+       // 1. get max and min over all basis function positions
+       assert( wfn->basis()->ncenter() > 0 );
+       SCVector3 min( wfn->basis()->r(0,0), wfn->basis()->r(0,1), wfn->basis()->r(0,2) );
+       SCVector3 max( wfn->basis()->r(0,0), wfn->basis()->r(0,1), wfn->basis()->r(0,2) );
+       for (int nr = 1; nr< wfn->basis()->ncenter(); ++nr) {
+         for (int i=0; i < 3; ++i) {
+           if (wfn->basis()->r(nr,i) < min(i))
+             min(i) = wfn->basis()->r(nr,i);
+           if (wfn->basis()->r(nr,i) > max(i))
+             max(i) = wfn->basis()->r(nr,i);
          }
+       }
+       std::cout << "Min is at " << min << " and max is at " << max << std::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
+       double boundary = 5.;  // boundary extent around compact domain containing basis functions
+       double delta = 1.; // step width in density sampling
+
+       // for the moment we always generate a grid of full size
+       min.x() = data.sampled_grid.begin[0];
+       min.y() = data.sampled_grid.begin[1];
+       min.z() = data.sampled_grid.begin[2];
+       max.x() = min.x() + data.sampled_grid.size;
+       max.y() = min.y() + data.sampled_grid.size;
+       max.z() = min.z() + data.sampled_grid.size;
+       const int gridpoints = pow(2,data.sampled_grid.level);
+       delta = data.sampled_grid.size / (double) gridpoints;
+       std::cout << "Grid starts at " << min
+           << " and ends at " << max
+           << " with a delta of " << delta
+           << " to get " << gridpoints << " gridpoints."
+           << std::endl;
+
+       // 3. sample the atomic density
+       int nbasis = wfn->basis()->nbasis();
+       double *b_val = new double[nbasis];
+       Ref<Integral> intgrl = Integral::get_default_integral();
+       GaussianBasisSet::ValueData vdat(wfn->basis(), integral);
+       SCVector3 r;
+       data.sampled_grid.sampled_grid.clear();
+       data.sampled_grid.sampled_grid.reserve(gridpoints*gridpoints*gridpoints);
+       for (r.x() = min.x() ; r.x() < max.x(); r.x() += delta)
+         for (r.y() = min.y() ; r.y() < max.y(); r.y() += delta)
+           for (r.z() =min.z() ; r.z() < max.z(); r.z() += delta) {
+             wfn->basis()->values(r, &vdat, b_val);
+             double sum = 0.;
+             for (int i=0; i<nbasis; i++)
+               sum += wfn->ao_density()->get_element(i,i)*b_val[i];
+             data.sampled_grid.sampled_grid.push_back(sum);
+           }
+     }
 
 //     // GaussianShell is the actual orbital functions it seems ...