Changeset 4dfef04

Timestamp:

Aug 10, 2012, 9:34:32 PM (13 years ago)

Author:

Frederik Heber <heber@…>

Children:

499cea

Parents:

d901f7

Message:

FIX: We now always convert from angstroem to bohr radii and vice versa, also rescale density to correct electron number.

File:

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

src/bin/mpqc/mpqc.cc

@@ -1638 +1638 @@
      }
      {
-       // fill positions and charges
+       // fill positions and charges (converting from bohr radii to angstroem)
+       const double AtomicLengthToAngstroem = 0.52917721;
        data.positions.reserve(wfn->molecule()->natom());
        data.charges.reserve(wfn->molecule()->natom());
@@ -1645 +1646 @@
          std::vector<double> pos(3, 0.);
          for (int j=0;j<3;++j)
-           pos[j] = wfn->molecule()->r(iatom, j);
+           pos[j] = wfn->molecule()->r(iatom, j)*AtomicLengthToAngstroem;
          data.positions.push_back(pos);
        }
@@ -1675 +1676 @@
 
        // for the moment we always generate a grid of full size
+       // (converting grid dimensions from angstroem to bohr radii)
        const double AtomicLengthToAngstroem = 0.52917721;
        SCVector3 min;
@@ -1685 +1687 @@
        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;
+       delta = (data.sampled_grid.size/AtomicLengthToAngstroem) / (double) gridpoints;
        std::cout << "Grid starts at " << min
            << " and ends at " << max
@@ -1696 +1698 @@
        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) {
+       r = min;
+       const double element_volume_conversion =
+           1./AtomicLengthToAngstroem/AtomicLengthToAngstroem/AtomicLengthToAngstroem;
+       for (size_t x=0; x< gridpoints; ++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) {
-             // only evaluate if we close the basis functions
+         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))
@@ -1706 +1713 @@
                data.sampled_grid.sampled_grid.push_back(0.);
              } else {
-               const double dens_at_r = scf->density(r);
+               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;
@@ -1712 +1719 @@
              }
            }
+           r.z() = min.z();
+         }
+         r.y() = min.y();
        }
-       assert( data.sampled_grid.size() == gridpoints*gridpoints*gridpoints);
+       assert( data.sampled_grid.sampled_grid.size() == gridpoints*gridpoints*gridpoints);
+       // normalization of electron charge to equal electron number
+       {
+         double integral_value = 0.;
+         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);
+         const double normalization = scf->nelectron()/integral_value;
+         std::cout << "Created " << data.sampled_grid.sampled_grid.size() << " grid points "
+             << " with integral value of " << integral_value
+             << " against nelectron of " << scf->nelectron() << "." << std::endl;
+         for (std::vector<double>::iterator diter = data.sampled_grid.sampled_grid.begin();
+             diter != data.sampled_grid.sampled_grid.end(); ++diter)
+             *diter *= normalization;
+       }
      }
      scf = 0;