Changeset 3ad723 for src

Timestamp:

Feb 15, 2013, 4:48:27 PM (13 years ago)

Author:

Frederik Heber <heber@…>

Children:

8b3553

Parents:

cf3275

git-author:

Frederik Heber <heber@…> (01/25/13 15:29:34)

git-committer:

Frederik Heber <heber@…> (02/15/13 16:48:27)

Message:

Rewrote sampled value gathering due to new non-zero window.

File:

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

src/bin/mpqc/mpqc.cc

@@ -1667 +1667 @@
        // 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
+       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
@@ -1675 +1684 @@
        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);
+       SCVector3 delta;
+       SCVector3 samplepoints;
+       SCVector3 length;
+       SCVector3 total;
        // 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
-       delta = (data.sampled_grid.size/AtomicLengthToAngstroem) / (double) gridpoints;
+       const int gridpoints = data.sampled_grid.getGridPointsPerAxis();
+       for (size_t i=0;i<3;++i) {
+         min(i) = bmin(i)/AtomicLengthToAngstroem;
+         max(i) = bmax(i)/AtomicLengthToAngstroem;
+         length(i) = bmax(i) - bmin(i);
+         total(i) = data.sampled_grid.end[i] - data.sampled_grid.begin[i];
+         delta(i) = total(i)/AtomicLengthToAngstroem/(double)gridpoints;
+         samplepoints(i) = floor(gridpoints*(length(i)/total(i))/AtomicLengthToAngstroem)+1;
+       }
        std::cout << "Grid starts at " << min
            << " and ends at " << max
            << " with a delta of " << delta
-           << " to get " << gridpoints << " gridpoints."
+           << " to get " << samplepoints << " samplepoints."
            << std::endl;
+       assert( data.sampled_grid.sampled_grid.size() == samplepoints(0)*samplepoints(1)*samplepoints(2));
 
        // 3. sample the atomic density
-       SCVector3 r;
-       data.sampled_grid.sampled_grid.clear();
-       data.sampled_grid.sampled_grid.reserve(gridpoints*gridpoints*gridpoints);
-       r = min;
        const double element_volume_conversion =
            1./AtomicLengthToAngstroem/AtomicLengthToAngstroem/AtomicLengthToAngstroem;
-       for (size_t x=0; x< gridpoints; ++x, r.x() += delta) {
+       SCVector3 r = min;
+       std::vector<double>::iterator griditer = data.sampled_grid.sampled_grid.begin();
+       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< gridpoints; ++y, r.y() += delta) {
-           for (size_t z=0; z< gridpoints; ++z, r.z() += delta) {
-             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.z()+boundary)))) {
-               data.sampled_grid.sampled_grid.push_back(0.);
-             } else {
-               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;
-               data.sampled_grid.sampled_grid.push_back(dens_at_r);
-             }
+         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)) {
+             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();
@@ -1718 +1728 @@
          r.y() = min.y();
        }
-       assert( data.sampled_grid.sampled_grid.size() == gridpoints*gridpoints*gridpoints);
+       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 *= pow(delta*AtomicLengthToAngstroem,3);
+         integral_value *= volume_element;
          const double normalization =
              ((integral_value == 0) && (scf->nelectron() == 0)) ?
                  1. : scf->nelectron()/integral_value;
-         std::cout << "Created " << data.sampled_grid.sampled_grid.size() << " grid points "
+         std::cout << "Created " << data.sampled_grid.sampled_grid.size() << " grid points"
              << " with integral value of " << integral_value
              << " against nelectron of " << scf->nelectron() << "." << std::endl;