source: src/Actions/FragmentationAction/FragmentationAutomationAction.cpp@ 82034f1

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 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/>.
 */

/*
 * FragmentationAutomationAction.cpp
 *
 *  Created on: May 18, 2012
 *      Author: heber
 */

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

#include <boost/archive/text_iarchive.hpp>
// boost asio needs specific operator new
#include <boost/asio.hpp>

#include "CodePatterns/MemDebug.hpp"

//// include headers that implement a archive in simple text format
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>

//
//#include <boost/mpl/remove.hpp>
//#include <boost/lambda/lambda.hpp>

//#include <iostream>

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"

#ifdef HAVE_JOBMARKET
#include "JobMarket/Jobs/FragmentJob.hpp"
#else
#include "Jobs/JobMarket/FragmentJob.hpp"
#endif

#include "Fragmentation/Automation/FragmentJobQueue.hpp"
#ifdef HAVE_JOBMARKET
#include "Fragmentation/Automation/MPQCFragmentController.hpp"
#else
#include "Fragmentation/Automation/MPQCCommandFragmentController.hpp"
#endif
#include "Fragmentation/Exporters/ExportGraph_ToJobs.hpp"
#include "Fragmentation/Summation/Containers/FragmentationChargeDensity.hpp"
#include "Fragmentation/Summation/Containers/FragmentationLongRangeResults.hpp"
#include "Fragmentation/Summation/Containers/FragmentationResultContainer.hpp"
#include "Fragmentation/Summation/Containers/FragmentationShortRangeResults.hpp"
#include "Fragmentation/Summation/Containers/MPQCData.hpp"
#include "Fragmentation/KeySetsContainer.hpp"
#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
#include "Fragmentation/Automation/VMGDebugGridFragmentController.hpp"
#include "Fragmentation/Automation/VMGFragmentController.hpp"
#include "Fragmentation/Summation/Containers/GridDownsampler.hpp"
#include "Fragmentation/Summation/Containers/VMGData.hpp"
#include "Fragmentation/Summation/Containers/VMGDataFused.hpp"
#include "Fragmentation/Summation/Containers/VMGDataMap.hpp"
#include "Fragmentation/Summation/Containers/VMGData_printKeyNames.hpp"
#endif
#include "World.hpp"

#include <boost/bind.hpp>
#include <boost/thread.hpp>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>

#include <boost/mpl/for_each.hpp>

#include "Actions/FragmentationAction/FragmentationAutomationAction.hpp"

using namespace MoleCuilder;

// and construct the stuff
#include "FragmentationAutomationAction.def"
#include "Action_impl_pre.hpp"
/** =========== define the function ====================== */

class controller_AddOn;

// needs to be defined for using the FragmentController
controller_AddOn *getAddOn()
{
  return NULL;
}

static void updateSteps(Process &p, const size_t step, const size_t total)
{
  LOG(1, "There are " << step << " steps out of " << total << " done.");
  p.setCurrStep(step);
}

template <class ResultClass>
std::pair<double, double> getResiduals(const ResultClass &_result)
{ return std::pair<double, double>(0.,0.); }

std::pair<double, double> getResiduals(const MPQCData &_result)
{ return std::make_pair(_result.accuracy, _result.desired_accuracy); }

std::pair<double, double> getResiduals(const VMGData &_result)
{
  // sensibility check on absolute residual
  if (_result.residual > 1e-5)
    ELOG(2, "Encountered absolute residual greater than 1e-5: " << _result.residual);
  // take smaller value of the two as just one of them needs to be less than precision
  const double residual = std::min(_result.relative_residual, _result.residual);
  return std::make_pair(residual, _result.precision);
}

template <class ResultClass>
bool checkResults(const typename std::map<JobId_t, ResultClass> &_results)
{
  bool result = true;
  for (typename std::map<JobId_t, ResultClass>::const_iterator iter = _results.begin();
      iter != _results.end(); ++iter) {
    std::pair<double, double> residuals = getResiduals(iter->second);
    if (residuals.second != 0.) {
      if (residuals.first >= residuals.second) {
        ELOG(1, "Fragment Job " << iter->first << " converged to " << residuals.first
            << " instead of " << residuals.second);
        result = false;
      }
    } else {
      LOG(4, "DEBUG: Not checking accuracy because desired precision not given.");
    }
  }
  return result;
}


#ifdef HAVE_JOBMARKET
/** This performs a direct sum calculation for the nuclei-nuclei interaction energy
 * and forces for each fragment contained in \a shortrangedata.
 *
 * Results are (re)placed in \a _longrangedata.
 *
 * \warning This is only valid for fully open boundary conditions.
 *
 * \param _shortrangedata containing information on position and charges per fragment
 * \param _longrangedata containing energy and forces on return
 */
static void calculateNucleiNucleiLongRangeContribution(
    const std::map<JobId_t, MPQCData> &_shortrangedata,
    std::map<JobId_t, VMGData>& _longrangedata
    )
{
  ASSERT( _shortrangedata.size() == _longrangedata.size(),
      "calculateNucleiNucleiLongRangeContribution() - shortrange and longrange data have inequal size.");
  std::map<JobId_t, VMGData>::iterator longrangeiter = _longrangedata.begin();
  std::map<JobId_t, MPQCData>::const_iterator iter = _shortrangedata.begin();
  for (;iter != _shortrangedata.end(); ++iter, ++longrangeiter) {
    const MPQCData &data = iter->second;
    VMGData &longrange_data = longrangeiter->second;
    // set long-range contributions to zero
    longrange_data.nuclei_long = 0.;
    longrange_data.forces.clear();
    longrange_data.forces.resize(data.positions.size(), FragmentForces::force_t(3,0.));
    longrange_data.hasForces = true;
    // go through positions and evaluate sum naively
    ASSERT( data.positions.size() == data.charges.size(),
        "calculateNucleiNucleiLongRangeContribution() - positions and charges differ in size.");
    std::vector< std::vector<double> >::const_iterator positer = data.positions.begin();
    std::vector< double >::const_iterator chargeiter = data.charges.begin();
    FragmentForces::iterator forceiter = longrange_data.forces.begin();
    for (; (positer != data.positions.end()) && (chargeiter != data.charges.end());
        ++positer, ++chargeiter, ++forceiter) {
      ASSERT( positer->size() == NDIM,
          "calculateNucleiNucleiLongRangeContribution() - position "
          +toString(std::distance(data.positions.begin(), positer))+" has not 3 components.");
      const Vector position((*positer)[0], (*positer)[1], (*positer)[2]);

      std::vector< std::vector<double> >::const_iterator otherpositer = data.positions.begin();
      std::vector< double >::const_iterator otherchargeiter = data.charges.begin();
      for (; (otherpositer != data.positions.end()) && (otherchargeiter != data.charges.end());
          ++otherpositer, ++otherchargeiter) {
        if ((positer == otherpositer) || (chargeiter == otherchargeiter))
          continue;
        ASSERT( otherpositer->size() == NDIM,
            "calculateNucleiNucleiLongRangeContribution() - other position "
          +toString(std::distance(data.positions.begin(), otherpositer))+" has not 3 components.");
        const Vector otherposition((*otherpositer)[0], (*otherpositer)[1], (*otherpositer)[2]);
        const Vector distance = position - otherposition;
        const double invsqrdist = 1./distance.Norm();
        const double factor = (*otherchargeiter)*invsqrdist*invsqrdist*invsqrdist;
        (*forceiter)[0] = factor*distance[0];
        (*forceiter)[1] = factor*distance[1];
        (*forceiter)[2] = factor*distance[2];
        longrange_data.nuclei_long += 0.5*(*chargeiter)*(*otherchargeiter)*invsqrdist;
      }
    }
  }
}
#endif

ActionState::ptr FragmentationFragmentationAutomationAction::performCall() {
  boost::asio::io_service io_service;

  // TODO: Have io_service run in second thread and merge with current again eventually

  FragmentationResultContainer &container =
      FragmentationResultContainer::getInstance();
  const KeySetsContainer& keysets = FragmentJobQueue::getInstance().getKeySets();
  const KeySetsContainer& forcekeysets = FragmentJobQueue::getInstance().getFullKeySets();

  size_t Exitflag = 0;
  std::map<JobId_t, MPQCData> shortrangedata;
#ifdef HAVE_JOBMARKET
  {
    const size_t NumberJobs = FragmentJobQueue::getInstance().size();
    MPQCFragmentController mpqccontroller(io_service);
    mpqccontroller.setHost(params.host.get());
    mpqccontroller.setPort(params.port.get());
    // Phase One: obtain ids
    mpqccontroller.requestIds(NumberJobs);
    if (mpqccontroller.getExitflag() != 0)
      return Action::failure;

    // Phase Two: add MPQCJobs and send
    const bool AddJobsStatus = mpqccontroller.addJobsFromQueue(
        params.DoLongrange.get() ? MPQCData::DoSampleDensity : MPQCData::DontSampleDensity,
        params.DoValenceOnly.get() ? MPQCData::DoSampleValenceOnly : MPQCData::DontSampleValenceOnly
        );
    if (AddJobsStatus)
      LOG(1, "INFO: Added jobs from FragmentJobsQueue.");
    else {
      ELOG(1, "Adding jobs failed.");
      return Action::failure;
    }
    mpqccontroller.run();

    // Phase Three: calculate result
    setMaxSteps(NumberJobs);
    mpqccontroller.setUpdateHandler(
        boost::bind(&updateSteps, boost::ref(*this), _1, _2)
        );
    start();
    boost::thread wait_thread(
        boost::bind(&MPQCFragmentController::waitforResults, boost::ref(mpqccontroller), boost::cref(NumberJobs))
    );
    wait_thread.join();
    stop();
    if (mpqccontroller.getExitflag() != 0)
      return Action::failure;

    mpqccontroller.getResults(shortrangedata);
    if (!checkResults(shortrangedata)) {
      ELOG(1, "At least one of the fragments' energy could not be properly minimized.");
      return Action::failure;
    }
    if (mpqccontroller.getExitflag() != 0) {
      ELOG(1, "MPQCCommandFragmentController returned non-zero exit flag.");
      return Action::failure;
    }

    Exitflag += mpqccontroller.getExitflag();
  }
#else
  {
    const size_t NumberJobs = FragmentJobQueue::getInstance().size();
    MPQCCommandFragmentController mpqccontroller;
    // Phase One: obtain ids: not needed, we have infinite pool

    // Phase Two: add MPQCJobs and send
    const size_t NoJobs = mpqccontroller.addJobsFromQueue(
        params.DoLongrange.get() ? MPQCData::DoSampleDensity : MPQCData::DontSampleDensity,
        params.DoValenceOnly.get() ? MPQCData::DoSampleValenceOnly : MPQCData::DontSampleValenceOnly,
        params.executable.get().string()
        );
    if (NoJobs != NumberJobs) {
      ELOG(1, "Not all jobs were added succesfully.");
      return Action::failure;
    }
    LOG(1, "INFO: Added " << NoJobs << " from FragmentJobsQueue.");

    // prepare process
    setMaxSteps(NumberJobs);
    mpqccontroller.setUpdateHandler(
        boost::bind(&updateSteps, boost::ref(*this), _1, _2)
        );
    start();
    mpqccontroller.run();
    stop();
    if (mpqccontroller.getExitflag() != 0) {
      ELOG(1, "MPQCCommandFragmentController returned non-zero exit flag before getting results.");
      return Action::failure;
    }

    // get back the results and place them in shortrangedata
    mpqccontroller.getResults(shortrangedata);
    ASSERT( shortrangedata.size() == NumberJobs,
        "FragmentationFragmentationAutomationAction::performCall() - number of converted results "
        +toString(shortrangedata.size())+" and number of jobs "+toString(NumberJobs)+ " differ.");
    if (!checkResults(shortrangedata)) {
      ELOG(1, "At least one of the fragments' energy could not be properly minimized.");
      return Action::failure;
    }
    if (mpqccontroller.getExitflag() != 0) {
      ELOG(1, "MPQCCommandFragmentController returned non-zero exit flag after getting results.");
      return Action::failure;
    }

    Exitflag += mpqccontroller.getExitflag();
  }
#endif

#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
  if (params.DoLongrange.get()) {
  if ( const_cast<const World &>(World::getInstance()).getAllAtoms().size() == 0) {
    STATUS("Please load the full molecule into the world before starting this action.");
    return Action::failure;
  }

  // obtain combined charge density
  std::map<JobId_t, MPQCData> shortrangedata_downsampled;
  SamplingGridProperties domain(ExportGraph_ToJobs::getDomainGrid(params.level.get()));
  FragmentationChargeDensity summedChargeDensity(
      shortrangedata);
  {
    SamplingGrid zero_globalgrid(domain);
    summedChargeDensity(
        shortrangedata,
        FragmentJobQueue::getInstance().getKeySets(),
        zero_globalgrid);
  }
  const std::vector<SamplingGrid> full_sample = summedChargeDensity.getFullSampledGrid();
  LOG(1, "INFO: There are " << shortrangedata.size() << " short-range and "
      << full_sample.size() << " level-wise long-range jobs.");

  // check boundary conditions
  const BoundaryConditions::Conditions_t &conditions =
      World::getInstance().getDomain().getConditions();
  const bool OpenBoundaryConditions =
      !((conditions[0] == BoundaryConditions::Wrap) &&
      (conditions[1] == BoundaryConditions::Wrap) &&
      (conditions[2] == BoundaryConditions::Wrap));
  LOG(1, std::string("INFO: Using ") 
      << (OpenBoundaryConditions ? "open" : "periodic")
      << " boundary conditions.");

  // Phase Four: obtain more ids
  std::map<JobId_t, VMGData> longrangedata;
  {
    VMGFragmentController vmgcontroller(io_service);
    vmgcontroller.setHost(params.host.get());
    vmgcontroller.setPort(params.port.get());
    const size_t NoJobs = shortrangedata.size()+full_sample.size();
    vmgcontroller.requestIds(2*NoJobs);
    if (vmgcontroller.getExitflag() != 0) {
      ELOG(1, "VMGFragmentController returned non-zero exit flag on requesting long-range ids.");
      return Action::failure;
    }

    // Phase Five a: create VMGJobs for electronic charge distribution
    const size_t near_field_cells = params.near_field_cells.get();
    const size_t interpolation_degree = params.interpolation_degree.get();
    if (!vmgcontroller.createLongRangeJobs(
        shortrangedata,
        full_sample,
        near_field_cells,
        interpolation_degree,
        VMGFragmentController::DontSampleParticles,
        VMGFragmentController::DoTreatGrid,
        params.DoValenceOnly.get() ? MPQCData::DoSampleValenceOnly : MPQCData::DontSampleValenceOnly,
        params.DoPrintDebug.get(),
        OpenBoundaryConditions,
        params.DoSmearCharges.get(),
        false)) {
      STATUS("Could not create long-range jobs for electronic charge distribution.");
      return Action::failure;
    }

    // Phase Six a: calculate result
    vmgcontroller.waitforResults(NoJobs);
    if (vmgcontroller.getExitflag() != 0) {
      ELOG(1, "VMGFragmentController returned non-zero exit flag before getting results.");
      return Action::failure;
    }
    vmgcontroller.getResults(longrangedata);
    ASSERT( NoJobs == longrangedata.size(),
        "FragmentationFragmentationAutomationAction::performCall() - number of MPQCresults+"
        +toString(full_sample.size())+"="+toString(NoJobs)
        +" and first VMGresults "+toString(longrangedata.size())+" don't match.");
    if (!checkResults(longrangedata)) {
      ELOG(1, "At least one of the fragments' electronic long-range potential could not be properly minimized.");
      return Action::failure;
    }
    Exitflag += vmgcontroller.getExitflag();

    {
      std::map<JobId_t, VMGData> longrangedata_both;
      // Phase Five b: create VMGJobs for nuclei charge distributions
      const size_t near_field_cells = params.near_field_cells.get();
      const size_t interpolation_degree = params.interpolation_degree.get();
      if (!vmgcontroller.createLongRangeJobs(
          shortrangedata,
          full_sample,
          near_field_cells,
          interpolation_degree,
          VMGFragmentController::DoSampleParticles,
          VMGFragmentController::DoTreatGrid,
          params.DoValenceOnly.get() ? MPQCData::DoSampleValenceOnly : MPQCData::DontSampleValenceOnly,
          params.DoPrintDebug.get(),
          OpenBoundaryConditions,
          params.DoSmearCharges.get(),
          params.UseImplicitCharges.get())) {
        STATUS("Could not create long-range jobs for nuclei charge distribution.");
        return Action::failure;
      }

      // Phase Six b: calculate result
      vmgcontroller.waitforResults(NoJobs);
      if (vmgcontroller.getExitflag() != 0)
        return Action::failure;
      vmgcontroller.getResults(longrangedata_both);
      ASSERT( NoJobs == longrangedata_both.size(),
          "FragmentationFragmentationAutomationAction::performCall() - number of MPQCresults+"
          +toString(full_sample.size())+"="+toString(NoJobs)
          +" and second VMGresults "+toString(longrangedata_both.size())+" don't match.");
      if (!checkResults(longrangedata_both)) {
        ELOG(1, "At least one of the fragments' nuclei long-range potential could not be properly minimized.");
        return Action::failure;
      }
      if (vmgcontroller.getExitflag() != 0) {
        ELOG(1, "VMGFragmentController returned non-zero exit flag after getting results.");
        return Action::failure;
      }
      Exitflag += vmgcontroller.getExitflag();

      // now replace nuclei-nuclei contribution with values from direct summation (is exact for
      // open boundary conditions
      //
      // NOTE: We need to calculate both_sampled_potential for fitting partial charges,
      // nowhere else is this quantity needed.
      //if (OpenBoundaryConditions)
      //  calculateNucleiNucleiLongRangeContribution(shortrangedata, longrangedata_both);

      // go through either data and replace nuclei_long with contribution from both
      ASSERT( longrangedata.size() == longrangedata_both.size(),
          "FragmentationFragmentationAutomationAction::performCall() - longrange results have different sizes.");
      std::map<JobId_t, VMGData>::iterator destiter = longrangedata.begin();
      std::map<JobId_t, VMGData>::const_iterator srciter = longrangedata_both.begin();
      for (;destiter != longrangedata.end(); ++srciter, ++destiter) {
        destiter->second.both_sampled_potential = srciter->second.sampled_potential;
        destiter->second.nuclei_long = srciter->second.nuclei_long;
        destiter->second.forces = srciter->second.forces;
        destiter->second.hasForces = srciter->second.hasForces;
      }
    }
  }

  if (params.DoPrintDebug.get()) {
    std::map<JobId_t, std::string> debugData;
    {
      if (!full_sample.empty()) {
        // create debug jobs for each level to print the summed-up potential to vtk files
        VMGDebugGridFragmentController debugcontroller(io_service);
        debugcontroller.setHost(params.host.get());
        debugcontroller.setPort(params.port.get());
        debugcontroller.requestIds(full_sample.size());
        if (!debugcontroller.createDebugJobs(full_sample, OpenBoundaryConditions)) {
          STATUS("Could not create debug jobs.");
          return Action::failure;
        }
        debugcontroller.waitforResults(full_sample.size());
        debugcontroller.getResults(debugData);
        Exitflag += debugcontroller.getExitflag();
      }
    }
  }
  container.clear();
  container.addFullResults(keysets, forcekeysets, shortrangedata, longrangedata);
  } else {
    container.clear();
    container.addShortRangeResults(keysets, forcekeysets, shortrangedata);
  }
#else
  container.clear();
  container.addShortRangeResults(keysets, forcekeysets, shortrangedata);
#endif

  // now clear all present jobs as we are done
  FragmentJobQueue::getInstance().clear();

  if (Exitflag != 0)
    STATUS("Controller has returned failure.");
  return (Exitflag == 0) ? Action::success : Action::failure;
}

ActionState::ptr FragmentationFragmentationAutomationAction::performUndo(ActionState::ptr _state) {
  return Action::success;
}

ActionState::ptr FragmentationFragmentationAutomationAction::performRedo(ActionState::ptr _state){
  return Action::success;
}

bool FragmentationFragmentationAutomationAction::canUndo() {
  return false;
}

bool FragmentationFragmentationAutomationAction::shouldUndo() {
  return false;
}
/** =========== end of function ====================== */