/*
* 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 .
*/
/*
* FragmentationAutomationAction.cpp
*
* Created on: May 18, 2012
* Author: heber
*/
// include config.h
#ifdef HAVE_CONFIG_H
#include
#endif
#include
// boost asio needs specific operator new
#include
//#include "CodePatterns/MemDebug.hpp"
//// include headers that implement a archive in simple text format
#include
#include
//
//#include
//#include
//#include
#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
#include
#include
#include
#include
#include
#include
#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
std::pair getResiduals(const ResultClass &_result)
{ return std::pair(0.,0.); }
std::pair getResiduals(const MPQCData &_result)
{ return std::make_pair(_result.accuracy, _result.desired_accuracy); }
std::pair 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
bool checkResults(const typename std::map &_results)
{
bool result = true;
for (typename std::map::const_iterator iter = _results.begin();
iter != _results.end(); ++iter) {
std::pair 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
*/
/* Is not used currently, commented out to disable warning message.
static void calculateNucleiNucleiLongRangeContribution(
const std::map &_shortrangedata,
std::map& _longrangedata
)
{
ASSERT( _shortrangedata.size() == _longrangedata.size(),
"calculateNucleiNucleiLongRangeContribution() - shortrange and longrange data have inequal size.");
std::map::iterator longrangeiter = _longrangedata.begin();
std::map::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 >::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 >::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();
const FragmentJobQueue::edges_per_fragment_t &edges_per_fragment =
FragmentJobQueue::getInstance().getEdgesPerFragment();
size_t Exitflag = 0;
std::map 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) {
ELOG(1, "Could not request all ids from Server.");
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) {
ELOG(1, "Could not obtain all results from Server.");
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(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 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 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 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 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) {
ELOG(1, "VMGFragmentController returned non-zero exit flag before getting results.");
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::iterator destiter = longrangedata.begin();
std::map::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 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, edges_per_fragment, shortrangedata, longrangedata);
} else {
container.clear();
container.addShortRangeResults(keysets, forcekeysets, edges_per_fragment, shortrangedata);
}
#else
container.clear();
container.addShortRangeResults(keysets, forcekeysets, edges_per_fragment, 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 ====================== */