Changeset 889c27

Timestamp:

Mar 11, 2013, 11:33:01 PM (13 years ago)

Author:

Frederik Heber <heber@…>

Children:

d13755

Parents:

8d3df9

git-author:

Frederik Heber <heber@…> (03/06/13 21:33:55)

git-committer:

Frederik Heber <heber@…> (03/11/13 23:33:01)

Message:

We sample electron density if MPQCData::DoLongrange says so.

File:

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

src/bin/mpqc/mpqc.cc

@@ -1631 +1631 @@
 //       }
      }
-     {
-       // fill positions and charges (NO LONGER converting from bohr radii to angstroem)
-       const double AtomicLengthToAngstroem = 1.;//0.52917721;
-       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)*AtomicLengthToAngstroem;
-         data.positions.push_back(pos);
+     // we do sample the density only on request
+     if (data.DoLongrange) {
+       {
+         // fill positions and charges (NO LONGER converting from bohr radii to angstroem)
+         const double AtomicLengthToAngstroem = 1.;//0.52917721;
+         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)*AtomicLengthToAngstroem;
+           data.positions.push_back(pos);
+         }
+         std::cout << "We have "
+             << data.positions.size() << " positions and "
+             << data.charges.size() << " charges." << std::endl;
        }
-       std::cout << "We have "
-           << data.positions.size() << " positions and "
-           << data.charges.size() << " charges." << std::endl;
-     }
-     if (data.sampled_grid.level != 0)
-     {
-       // we now need to sample the density on the grid
-       // 1. get max and min over all basis function positions
-       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 (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);
+       if (data.sampled_grid.level != 0)
+       {
+         // we now need to sample the density on the grid
+         // 1. get max and min over all basis function positions
+         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 (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 << "Basis min is at " << bmin << " and max is at " << bmax << 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
-       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
-       // (NO LONGER converting grid dimensions from angstroem to bohr radii)
-       const double AtomicLengthToAngstroem = 1.;//0.52917721;
-       SCVector3 min;
-       SCVector3 max;
-       SCVector3 delta;
-       size_t samplepoints[3];
-       // 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
-       const int gridpoints = data.sampled_grid.getGridPointsPerAxis();
-       for (size_t i=0;i<3;++i) {
-         min(i) = data.sampled_grid.begin_window[i]/AtomicLengthToAngstroem;
-         max(i) = data.sampled_grid.end_window[i]/AtomicLengthToAngstroem;
-         delta(i) = data.sampled_grid.getDeltaPerAxis(i)/AtomicLengthToAngstroem;
-         samplepoints[i] = data.sampled_grid.getWindowGridPointsPerAxis(i);
-       }
-       std::cout << "Grid starts at " << min
-           << " and ends at " << max
-           << " with a delta of " << delta
-           << " to get "
-           << samplepoints[0] << ","
-           << samplepoints[1] << ","
-           << samplepoints[2] << " samplepoints."
-           << std::endl;
-       assert( data.sampled_grid.sampled_grid.size() == samplepoints[0]*samplepoints[1]*samplepoints[2]);
-
-       // 3. sample the atomic density
-       const double element_volume_conversion =
-           1./AtomicLengthToAngstroem/AtomicLengthToAngstroem/AtomicLengthToAngstroem;
-       SCVector3 r = min;
-
-       // taken from scf::density()
-       RefSCMatrix nos
-           = scf->integral()->petite_list()->evecs_to_AO_basis(scf->natural_orbitals());
-       RefDiagSCMatrix nd = scf->natural_density();
-       GaussianBasisSet::ValueData *valdat
-         = new GaussianBasisSet::ValueData(scf->basis(), scf->integral());
-       std::vector<double>::iterator griditer = data.sampled_grid.sampled_grid.begin();
-       const int nbasis = scf->basis()->nbasis();
-       double *bs_values = new double[nbasis];
-
-       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 < samplepoints[1]; ++y, r.y() += delta(1)) {
-           for (size_t z = 0; z < samplepoints[2]; ++z, r.z() += delta(2)) {
-             scf->basis()->values(r,valdat,bs_values);
-
-             // loop over natural orbitals adding contributions to elec_density
-             double elec_density=0.0;
-             for (int i=0; i<nbasis; ++i) {
-                 double tmp = 0.0;
-                 for (int j=0; j<nbasis; ++j) {
-                     tmp += nos.get_element(j,i)*bs_values[j];
-                   }
-                 elec_density += nd.get_element(i)*tmp*tmp;
-               }
-             const double dens_at_r = elec_density * element_volume_conversion;
-//             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;
+         std::cout << "Basis min is at " << bmin << " and max is at " << bmax << 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
+         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
+         // (NO LONGER converting grid dimensions from angstroem to bohr radii)
+         const double AtomicLengthToAngstroem = 1.;//0.52917721;
+         SCVector3 min;
+         SCVector3 max;
+         SCVector3 delta;
+         size_t samplepoints[3];
+         // 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
+         const int gridpoints = data.sampled_grid.getGridPointsPerAxis();
+         for (size_t i=0;i<3;++i) {
+           min(i) = data.sampled_grid.begin_window[i]/AtomicLengthToAngstroem;
+           max(i) = data.sampled_grid.end_window[i]/AtomicLengthToAngstroem;
+           delta(i) = data.sampled_grid.getDeltaPerAxis(i)/AtomicLengthToAngstroem;
+           samplepoints[i] = data.sampled_grid.getWindowGridPointsPerAxis(i);
+         }
+         std::cout << "Grid starts at " << min
+             << " and ends at " << max
+             << " with a delta of " << delta
+             << " to get "
+             << samplepoints[0] << ","
+             << samplepoints[1] << ","
+             << samplepoints[2] << " samplepoints."
+             << std::endl;
+         assert( data.sampled_grid.sampled_grid.size() == samplepoints[0]*samplepoints[1]*samplepoints[2]);
+
+         // 3. sample the atomic density
+         const double element_volume_conversion =
+             1./AtomicLengthToAngstroem/AtomicLengthToAngstroem/AtomicLengthToAngstroem;
+         SCVector3 r = min;
+
+         // taken from scf::density()
+         RefSCMatrix nos
+             = scf->integral()->petite_list()->evecs_to_AO_basis(scf->natural_orbitals());
+         RefDiagSCMatrix nd = scf->natural_density();
+         GaussianBasisSet::ValueData *valdat
+           = new GaussianBasisSet::ValueData(scf->basis(), scf->integral());
+         std::vector<double>::iterator griditer = data.sampled_grid.sampled_grid.begin();
+         const int nbasis = scf->basis()->nbasis();
+         double *bs_values = new double[nbasis];
+
+         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 < samplepoints[1]; ++y, r.y() += delta(1)) {
+             for (size_t z = 0; z < samplepoints[2]; ++z, r.z() += delta(2)) {
+               scf->basis()->values(r,valdat,bs_values);
+
+               // loop over natural orbitals adding contributions to elec_density
+               double elec_density=0.0;
+               for (int i=0; i<nbasis; ++i) {
+                   double tmp = 0.0;
+                   for (int j=0; j<nbasis; ++j) {
+                       tmp += nos.get_element(j,i)*bs_values[j];
+                     }
+                   elec_density += nd.get_element(i)*tmp*tmp;
+                 }
+               const double dens_at_r = elec_density * element_volume_conversion;
+  //             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();
            }
-           r.z() = min.z();
+           r.y() = min.y();
          }
-         r.y() = min.y();
-       }
-       delete[] bs_values;
-       delete valdat;
-       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 *= 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"
-             << " 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;
+         delete[] bs_values;
+         delete valdat;
+         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 *= 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"
+               << " 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;
+         }
        }
      }
@@ -1902 +1905 @@
   // put info how to sample the density into MPQCData
   MPQCData data(grid);
+  data.DoLongrange = DoLongrange; // set whether we sample the density
 // now call work horse
   mainFunction(grp, values, output, input, generic_input, in_char_array, argc, argv, data);