Changeset debc65 for src

Timestamp:

Oct 5, 2012, 10:58:05 AM (13 years ago)

Author:

Frederik Heber <heber@…>

Children:

b88be2

Parents:

cd2e34

Message:

FIX: Major bug in sampling the electronic density as window around basis functions was wrong.

• z part did check bmax.y() and not bmax.z().
• added check whether both no electrons and integral_value is zero. This bug came up because integral_value came up as zero, causing arithmetic exception. We rather keep it this way: If integral_value is zero, this is always strange and should be inspected. It is only acceptable when there truely are no electrons present.

File:

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

src/bin/mpqc/mpqc.cc

@@ -1696 +1696 @@
        max.z() = min.z() + data.sampled_grid.size/AtomicLengthToAngstroem;
        const int gridpoints = pow(2,data.sampled_grid.level);
+       // 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;
        std::cout << "Grid starts at " << min
@@ -1714 +1717 @@
          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.y()+boundary))) {
+             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 {
@@ -1737 +1740 @@
              integral_value += *diter;
          integral_value *= pow(delta*AtomicLengthToAngstroem,3);
-         const double normalization = scf->nelectron()/integral_value;
+         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 "
              << " with integral value of " << integral_value