source: src/Fragmentation/Exporters/ExportGraph_ToJobs.cpp@ f8dae2

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2011 University of Bonn. All rights reserved.
 * Copyright (C)  2013 Frederik Heber. All rights reserved.
 *
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * ExportGraph_ToJobs.cpp
 *
 *  Created on: 08.03.2012
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

// boost asio required before MemDebug due to placement new
#include <boost/asio.hpp>

#include "CodePatterns/MemDebug.hpp"

#include "Fragmentation/Exporters/ExportGraph_ToJobs.hpp"

#include <algorithm>
#include <cmath>
#include <limits>

#include "CodePatterns/Log.hpp"

#include "Box.hpp"
#include "Fragmentation/KeySet.hpp"
#include "Fragmentation/Automation/FragmentJobQueue.hpp"
#include "Helpers/defs.hpp"
#ifdef HAVE_JOBMARKET
#include "Jobs/MPQCJob.hpp"
#else
#include "Jobs/MPQCCommandJob.hpp"
#endif
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "Parser/FormatParserStorage.hpp"
#include "World.hpp"

const double ExportGraph_ToJobs::log_two(log(2.));

ExportGraph_ToJobs::ExportGraph_ToJobs(
    const Graph &_graph,
    const enum HydrogenTreatment _treatment,
    const enum HydrogenSaturation _saturation,
    const SaturatedFragment::GlobalSaturationPositions_t &_globalsaturationpositions) :
    ExportGraph(_graph, _treatment, _saturation,_globalsaturationpositions),
    level(5),
    max_meshwidth(0.),
    minimum_empty_boundary(5./AtomicLengthToAngstroem)
{}

ExportGraph_ToJobs::~ExportGraph_ToJobs()
{}

SamplingGridProperties ExportGraph_ToJobs::getDomainGrid(
    const int _level)
{
  double domain_begin[NDIM] = { 0., 0., 0. };
  RealSpaceMatrix M = World::getInstance().getDomain().getM();
  M *= 1./AtomicLengthToAngstroem;  // scale to atomic length units
  const double size = std::max( std::max(M.at(0,0), M.at(1,1)), M.at(2,2));
  double domain_end[NDIM] = { size, size, size };
  SamplingGridProperties grid(domain_begin, domain_end, _level);
  return grid;
}

SamplingGridProperties ExportGraph_ToJobs::getGridExtentsForGivenBoundingBox(
    const std::pair<Vector, Vector> &_minmax,
    const SamplingGridProperties &_domain,
    const double & _minimum_empty_boundary
    )
{
  // return value
  SamplingGridProperties fragment_window;

  /// determine center of fragment
  const Vector center = .5*(_minmax.first + _minmax.second);
  LOG(2, "DEBUG: center of fragment is at " << center);

  /// associate center to its containing grid cell (defined by boundary points)
  Vector lower_center;
  Vector higher_center;
  bool equal_components[NDIM] = { false, false, false };
  for (size_t i=0;i<NDIM;++i) {
    lower_center[i] = _domain.getNearestLowerGridPoint(center[i], i);
    higher_center[i] = _domain.getNearestHigherGridPoint(center[i], i);
    equal_components[i] =
        fabs(lower_center[i] - higher_center[i]) < std::numeric_limits<double>::epsilon()*1e4;
  }
  LOG(2, "DEBUG: lower_center is " << lower_center << ", higher_center is " << higher_center);

  // fashion min max into cubic box extents of at least extent plus empty
  // boundary and at most three times the extent
  const Vector extent = _minmax.second - _minmax.first;
  double greatest_extent = extent[extent.GreatestComponent()];
  LOG(2, "DEBUG: extent of fragment is " << extent << " with greatest extent " << greatest_extent);
  if (greatest_extent > _minimum_empty_boundary)
    greatest_extent *= 3.;
  else
    greatest_extent += 2.* _minimum_empty_boundary;
  const Vector total(
      _domain.getTotalLengthPerAxis(0),
      _domain.getTotalLengthPerAxis(1),
      _domain.getTotalLengthPerAxis(2)
  );
  const double greatest_total = total[total.GreatestComponent()];
  greatest_extent = std::min(greatest_extent, greatest_total);
  LOG(2, "DEBUG: Using boundary of " << greatest_extent);

  /// increase levels around this center to find the matching window
  const double delta[NDIM] = {
      _domain.getDeltaPerAxis(0),
      _domain.getDeltaPerAxis(1),
      _domain.getDeltaPerAxis(2)
  };
  LOG(2, "DEBUG: delta is " << Vector(delta));
  const double greatest_delta = std::max(delta[0], std::max(delta[1], delta[2]));

  fragment_window.level = (int)ceil(log(greatest_extent/greatest_delta+1)/log_two);
  const size_t half_fragment_gridpoints = 1 << (fragment_window.level-1);
  const size_t domain_gridpoints = _domain.getGridPointsPerAxis();
  int begin_index[NDIM];
  int end_index[NDIM];
  int begin_steps[NDIM] = { 0, 0, 0 };
  int end_steps[NDIM] = { 0, 0, 0 };
  for (size_t i=0;i<NDIM;++i) {
    const unsigned int lower_index = round(lower_center[i]/delta[i]);
    const unsigned int higher_index = round(higher_center[i]/delta[i]);
    LOG(3, "DEBUG: center is within indices [" << lower_index << ", " << higher_index << "]");
    begin_index[i] = lower_index - (half_fragment_gridpoints-1);
    end_index[i] = higher_index + (half_fragment_gridpoints+(unsigned int)(equal_components[i]));
    LOG(3, "DEBUG: Intermediate begin end indices for axis #" << i << " are ["
        << begin_index[i] << ", " << end_index[i] << "]");
    if (begin_index[i] < 0)
      begin_steps[i] = -begin_index[i];
    if (end_index[i] > (int)domain_gridpoints)
      end_steps[i] = end_index[i] - domain_gridpoints;
    if ((begin_steps[i]+end_steps[i]+end_index[i]-begin_index[i] > (int)domain_gridpoints)
        || ((begin_steps[i] != 0) && (end_steps[i] != 0))) {
      begin_index[i] = 0;
      end_index[i] = domain_gridpoints;
      begin_steps[i] = 0;
      end_steps[i] = 0;
    }
    LOG(3, "DEBUG: Resulting final indices [" << begin_index[i] << ", "
        << end_index[i] << "] and step indices for axis #" << i << " are ["
        << begin_steps[i] << ", " << end_steps[i] << "]");
  }
  for (size_t i=0;i<NDIM;++i) {
    fragment_window.begin[i] = (begin_index[i]+begin_steps[i]-end_steps[i]) * delta[i];
    fragment_window.end[i] = (end_index[i]-end_steps[i]+begin_steps[i]) * delta[i];
  }
    LOG(2, "DEBUG: fragment begin is " << Vector(fragment_window.begin)
        << ", fragment end is " << Vector(fragment_window.end));

    /// check whether window is large enough, if not yet continue
#ifndef NDEBUG
    for (size_t i=0;i<NDIM;++i) {
      const double window_length = fragment_window.end[i] - fragment_window.begin[i];
      ASSERT (((greatest_total != greatest_extent) && (greatest_extent+delta[i] <= window_length))
          || ((greatest_total == greatest_extent) && (greatest_extent <= window_length)),
          "ExportGraph_ToJobs::getGridExtentsForGivenBoundingBox() - level "
          +toString(fragment_window.level)+" is insufficient to place fragment inside box "
          +toString(_minmax.first)+" <-> "+toString(_minmax.second));
    }
#endif
#ifndef NDEBUG
  for (size_t i=0;i<NDIM;++i) {
    ASSERT( fragment_window.begin[i] >= _domain.begin[i],
        "ExportGraph_ToJobs::getGridExtentsForGivenBoundingBox() - fragment begin "
        +toString(fragment_window.begin[i])+" is smaller than that of domain "
        +toString(_domain.begin[i]));
    ASSERT( fragment_window.end[i] <= _domain.end[i],
        "ExportGraph_ToJobs::getGridExtentsForGivenBoundingBox() - fragment end "
        +toString(fragment_window.end[i])+" is smaller than that of domain "
        +toString(_domain.end[i]));
  }
#endif
  if (DoLog(3))
    for (size_t i=0;i<NDIM;++i)
      LOG(3, "DEBUG: We have " << (fragment_window.end[i] - fragment_window.begin[i])/delta[i]
          << " gridpoints on axis " << i);

  return fragment_window;
}

void ExportGraph_ToJobs::setAcceptableFragmentLevel(
    SamplingGridProperties &_grid,
    const double &_max_meshwidth
    )
{
  double max_domain_meshwidth = 0.;
  for (size_t i=0;i<NDIM;++i)
    max_domain_meshwidth = std::max(max_domain_meshwidth, _grid.getDeltaPerAxis(i));
  // find power of 2 that's just greater than that ratio of given over desired mesh width
  const int desired_level = ceil(log(max_domain_meshwidth / _max_meshwidth)/log_two);
  _grid.level +=  desired_level;
  if (_grid.level < 2) {
    ELOG(1, "ExportGraph_ToJobs::setAcceptableFragmentLevel() - fragment grid level of 1 or less? "
        << "max-meshwidth " << _max_meshwidth << " chosen too large w.r.t delta of "
        << max_domain_meshwidth << "? Setting to minimum of 2.");
    _grid.level = 2;
  } else {
    LOG(4, "DEBUG: Fragment requires " << desired_level
        << " additional grid levels to reach at least " << _max_meshwidth
        << " from " << max_domain_meshwidth);
  }
}

bool ExportGraph_ToJobs::operator()()
{
  std::vector<FragmentJob::ptr> jobs;
  KeySetsContainer KeySets;
  KeySetsContainer FullKeySets;
  jobs.reserve(TotalGraph.size());
  LOG(1, "INFO: Creating " << TotalGraph.size() << " possible bond fragmentation jobs.");

  // checking that max-meshwidth is not smaller than grid spacing of global level
  const SamplingGridProperties grid(getDomainGrid(level));

  // gather info about the domain
  const double global_grid_delta =
      std::max(grid.getDeltaPerAxis(0),
          std::max(grid.getDeltaPerAxis(1),
              grid.getDeltaPerAxis(2)));
  if (global_grid_delta < max_meshwidth) {
    ELOG(2, "Desired max meswidth " << max_meshwidth
        << " is larger than grid spacing of global grid " << global_grid_delta << ".\n"
        << "We cannot use coarser grids in fragment level, equating mesh width to global delta.");
    max_meshwidth = global_grid_delta;
  }

  // set parser to use for writing job configurations
  const ParserTypes jobtype =
      FormatParserStorage::getInstance().getTypeFromName("mpqc");

  // go through all fragments, output to stream and create job therefrom
  ExportGraph::SaturatedFragment_ptr CurrentFragment = getNextFragment();
  for (; (CurrentFragment != NULL) && (CurrentFragment->getKeySet() != ExportGraph::EmptySet);
      CurrentFragment = getNextFragment()) {
    const KeySet &set = CurrentFragment->getKeySet();
    LOG(2, "INFO: Creating bond fragment job for set " << set << ".");

    SamplingGridProperties fragment_window(grid);
    if (max_meshwidth != 0.) {
      // gather extent of the fragment in atomic length(!)
      std::pair<Vector, Vector> minmax = CurrentFragment->getRoughBoundingBox();
      minmax.first *= 1./AtomicLengthToAngstroem;
      minmax.second *= 1./AtomicLengthToAngstroem;

      /** we need to find the begin and end points of the smaller grid
       * otherwise we would calculate the whole domain with an incrased grid level with just
       * a small window (which is only used for sampling and storing; vmg uses full grid).
       * Hence, we need to make the grid smaller and such that it is still in a power of two
       * and coinciding with the grid points of the global grid.
       */
      fragment_window = getGridExtentsForGivenBoundingBox(minmax, grid, minimum_empty_boundary);
      LOG(4, "DEBUG: Fragment " << CurrentFragment->getKeySet() << " has window from  "
          << Vector(fragment_window.begin) << " to " << Vector(fragment_window.end)
          << " with total level portion of " << fragment_window.level );

      // next we need to check how much we need to increase from the grid level for
      // the total domain to achieve at least maximum mesh width
      setAcceptableFragmentLevel(fragment_window, max_meshwidth/AtomicLengthToAngstroem);
    }

    // store config in stream
    {
      std::stringstream output;
      // save to stream
      CurrentFragment->OutputConfig(output, jobtype);
      // create job and insert
      FragmentJob::ptr testJob(
#ifdef HAVE_JOBMARKET
          new MPQCJob(
              JobId::IllegalJob,
              output.str(),
              fragment_window.begin, fragment_window.end, fragment_window.level)
#else
          new MPQCCommandJob(output.str(), JobId::IllegalJob)
#endif
      );
      testJob->setPriority(CurrentFragment->getKeySet().size());
      jobs.push_back(testJob);

      // order is the same as the number of non-hydrogen atoms
      const KeySet &keyset = CurrentFragment->getKeySet();
      const size_t order = keyset.size();
      const KeySet &fullmolecule = CurrentFragment->getFullMolecule();
      const KeySet &saturationhydrogens = CurrentFragment->getSaturationHydrogens();
      KeySetsContainer::IntVector indices(keyset.begin(), keyset.end());
      KeySetsContainer::IntVector forceindices(fullmolecule.begin(), fullmolecule.end());
      {
        // replace all saturated hydrogen indices by "-1"
        for (KeySetsContainer::IntVector::iterator iter = forceindices.begin();
            iter != forceindices.end();
            ++iter)
          if (saturationhydrogens.find(*iter) != saturationhydrogens.end())
            *iter = -1;
      }
      KeySets.insert(indices, order);
      FullKeySets.insert(forceindices, order);
    }
    // store force index reference file
    // explicitly release fragment
    CurrentFragment.reset();
  }
  if (CurrentFragment == NULL) {
    ELOG(1, "Some error while obtaining the next fragment occured.");
    return false;
  }

  // push final jobs
  FragmentJobQueue::getInstance().addJobs(jobs, KeySets, FullKeySets);

  return true;
}