Changeset d901f7

Timestamp:

Aug 9, 2012, 8:20:56 AM (13 years ago)

Author:

Frederik Heber <heber@…>

Children:

4dfef04

Parents:

1ed5fc

Message:

FIX: MBPT extractions were not working correctly.

• we need to obtain ref() from MBPT to continue with the remainder of the extraction.

File:

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

src/bin/mpqc/mpqc.cc

@@ -1497 +1497 @@
        }
      }
+     SCF *scf = NULL;
      {
        MBPT2 *mbpt2 = dynamic_cast<MBPT2*>(wfn.pointer());
        if (mbpt2 != NULL) {
          data.energies.correlation = mbpt2->corr_energy();
+         scf = mbpt2->ref().pointer();
          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 = wfn->overlap();
-       RefSymmSCMatrix cl_dens_ = wfn->ao_density();
+
+       RefSymmSCMatrix t = scf->overlap();
+       RefSymmSCMatrix cl_dens_ = scf->ao_density();
 
        SCFEnergy *eop = new SCFEnergy;
@@ -1527 +1533 @@
      {
        // taken from Wavefunction::core_hamiltonian()
-       RefSymmSCMatrix hao(wfn->basis()->basisdim(), wfn->basis()->matrixkit());
+       RefSymmSCMatrix hao(scf->basis()->basisdim(), scf->basis()->matrixkit());
        hao.assign(0.0);
-       Ref<PetiteList> pl = wfn->integral()->petite_list();
+       Ref<PetiteList> pl = scf->integral()->petite_list();
        Ref<SCElementOp> hc =
-         new OneBodyIntOp(new SymmOneBodyIntIter(wfn->integral()->kinetic(), pl));
+         new OneBodyIntOp(new SymmOneBodyIntIter(scf->integral()->kinetic(), pl));
        hao.element_op(hc);
        hc=0;
 
-       RefSymmSCMatrix h(wfn->so_dimension(), wfn->basis_matrixkit());
+       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_ = wfn->ao_density();
+       RefSymmSCMatrix cl_dens_ = scf->ao_density();
 
        SCFEnergy *eop = new SCFEnergy;
@@ -1556 +1562 @@
      }
      {
-       RefSymmSCMatrix t = wfn->core_hamiltonian();
-       RefSymmSCMatrix cl_dens_ = wfn->ao_density();
+       RefSymmSCMatrix t = scf->core_hamiltonian();
+       RefSymmSCMatrix cl_dens_ = scf->ao_density();
 
        SCFEnergy *eop = new SCFEnergy;
@@ -1599 +1605 @@
        }
      }
+     grad = NULL;
      {
        // times obtain from key "mpqc" which should be the first
@@ -1613 +1620 @@
        if (flops != NULL)
          data.times.flops = flops[0];
-       delete cpu_time;
-       delete wall_time;
-       delete flops;
+       delete[] cpu_time;
+       delete[] wall_time;
+       delete[] flops;
      }
 
@@ -1631 +1638 @@
      }
      {
+       // fill positions and charges
+       data.positions.reserve(wfn->molecule()->natom());
+       data.charges.reserve(wfn->molecule()->natom());
+       for (int iatom=0;iatom < wfn->molecule()->natom(); ++iatom) {
+         data.charges.push_back(wfn->molecule()->Z(iatom));
+         std::vector<double> pos(3, 0.);
+         for (int j=0;j<3;++j)
+           pos[j] = wfn->molecule()->r(iatom, j);
+         data.positions.push_back(pos);
+       }
+       std::cout << "We have "
+           << data.positions.size() << " positions and "
+           << data.charges.size() << " charges." << std::endl;
+     }
+     {
        // 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) {
+       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 (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);
+           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);
          }
        }
-       std::cout << "Min is at " << min << " and max is at " << max << std::endl;
+       std::cout << "Basis min is at " << bmin << " and max is at " << bmax << std::endl;
 
        // 2. choose an appropriately large grid
@@ -1653 +1675 @@
 
        // 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 double AtomicLengthToAngstroem = 0.52917721;
+       SCVector3 min;
+       SCVector3 max;
+       min.x() = data.sampled_grid.begin[0]/AtomicLengthToAngstroem;
+       min.y() = data.sampled_grid.begin[1]/AtomicLengthToAngstroem;
+       min.z() = data.sampled_grid.begin[2]/AtomicLengthToAngstroem;
+       max.x() = min.x() + data.sampled_grid.size/AtomicLengthToAngstroem;
+       max.y() = min.y() + data.sampled_grid.size/AtomicLengthToAngstroem;
+       max.z() = min.z() + data.sampled_grid.size/AtomicLengthToAngstroem;
        const int gridpoints = pow(2,data.sampled_grid.level);
        delta = data.sampled_grid.size / (double) gridpoints;
@@ -1668 +1693 @@
 
        // 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.x() = min.x() ; r.x() < max.x(); r.x() += delta) {
+         std::cout << "Sampling now for x=" << r.x() << std::endl;
          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);
+             // only evaluate if we close the basis functions
+             if (((r.x() < bmin.x()-boundary) || (r.x() > bmax.x()+boundary))
+               || ((r.y() < bmin.y()-boundary) || (r.y() > bmax.y()+boundary))
+               || ((r.z() < bmin.z()-boundary) || (r.z() > bmax.y()+boundary))) {
+               data.sampled_grid.sampled_grid.push_back(0.);
+             } else {
+               const double dens_at_r = scf->density(r);
+//               if (fabs(dens_at_r) > 1e-4)
+//                 std::cout << "Electron density at " << r << " is " << dens_at_r << std::endl;
+               data.sampled_grid.sampled_grid.push_back(dens_at_r);
+             }
            }
+       }
+       assert( data.sampled_grid.size() == gridpoints*gridpoints*gridpoints);
      }
+     scf = 0;
 
 //     // GaussianShell is the actual orbital functions it seems ...