/*
* Project: MoleCuilder
* Description: creates and alters molecular systems
* Copyright (C) 2010-2012 University of Bonn. 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
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "JobMarket/Controller/FragmentController.hpp"
#include "JobMarket/Jobs/FragmentJob.hpp"
#include "Atom/atom.hpp"
#include "Box.hpp"
#include "Fragmentation/EnergyMatrix.hpp"
#include "Fragmentation/ForceMatrix.hpp"
#include "Fragmentation/Fragmentation.hpp"
#include "Fragmentation/SetValues/Histogram.hpp"
#include "Fragmentation/SetValues/IndexedVectors.hpp"
#include "Fragmentation/HydrogenSaturation_enum.hpp"
#include "Fragmentation/KeySet.hpp"
#include "Fragmentation/KeySetsContainer.hpp"
#include "Fragmentation/Summation/printOrthogonalSum.hpp"
#include "Fragmentation/Summation/printSum.hpp"
#include "Graph/DepthFirstSearchAnalysis.hpp"
#include "Jobs/MPQCJob.hpp"
#include "Jobs/MPQCData.hpp"
#include "Jobs/MPQCData_printKeyNames.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#ifdef HAVE_VMG
#include "Jobs/VMGJob.hpp"
#endif
#include "molecule.hpp"
#include "World.hpp"
#include
#include
#include
#include
#include "Actions/FragmentationAction/FragmentationAutomationAction.hpp"
using namespace MoleCuilder;
using namespace boost::assign;
// 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;
}
const int LEVEL = 5;
/** Creates a MPQCCommandJob with argument \a filename.
*
* @param jobs created job is added to this vector
* @param command mpqc command to execute
* @param filename filename being argument to job
* @param nextid id for this job
*/
void parsejob(
std::vector &jobs,
const std::string &command,
const std::string &filename,
const JobId_t nextid)
{
std::ifstream file;
file.open(filename.c_str());
ASSERT( file.good(), "parsejob() - file "+filename+" does not exist.");
std::string output((std::istreambuf_iterator(file)),
std::istreambuf_iterator());
double begin[NDIM] = { 0., 0., 0. };
const RealSpaceMatrix& M = World::getInstance().getDomain().getM();
const double size = M.at(0,0);
ASSERT( M.determinant() == size*size*size,
"parsejob() - current domain matrix "+toString(M)+" is not cubic.");
const int level = LEVEL;
FragmentJob::ptr testJob( new MPQCJob(nextid, output, begin, size, level) );
jobs.push_back(testJob);
file.close();
LOG(1, "INFO: Added MPQCCommandJob from file "+filename+".");
}
/** Helper function to get number of atoms somehow.
*
* Here, we just parse the number of lines in the adjacency file as
* it should correspond to the number of atoms, except when some atoms
* are not bonded, but then fragmentation makes no sense.
*
* @param path path to the adjacency file
*/
size_t getNoAtomsFromAdjacencyFile(const std::string &path)
{
size_t NoAtoms = 0;
// parse in special file to get atom count (from line count)
std::string filename(path);
filename += FRAGMENTPREFIX;
filename += ADJACENCYFILE;
std::ifstream adjacency(filename.c_str());
if (adjacency.fail()) {
LOG(0, endl << "getNoAtomsFromAdjacencyFile() - Unable to open " << filename << ", is the directory correct?");
return false;
}
std::string buffer;
while (getline(adjacency, buffer))
NoAtoms++;
LOG(1, "INFO: There are " << NoAtoms << " atoms.");
return NoAtoms;
}
/** Extracts MPQCData from received vector of FragmentResults.
*
* @param results results to extract MPQCData from
* @param fragmentData on return array filled with extracted MPQCData
*/
void ConvertFragmentResultToMPQCData(
const std::vector &results,
std::vector &fragmentData)
{
// extract results
fragmentData.clear();
fragmentData.reserve(results.size());
LOG(2, "DEBUG: Parsing now through " << results.size() << " results.");
for (std::vector::const_iterator iter = results.begin();
iter != results.end(); ++iter) {
//LOG(1, "RESULT: job #"+toString((*iter)->getId())+": "+toString((*iter)->result));
MPQCData extractedData;
std::stringstream inputstream((*iter)->result);
LOG(2, "DEBUG: First 50 characters FragmentResult's string: "+(*iter)->result.substr(0, 50));
boost::archive::text_iarchive ia(inputstream);
ia >> extractedData;
LOG(1, "INFO: extracted data is " << extractedData << ".");
fragmentData.push_back(extractedData);
}
ASSERT( results.size() == fragmentData.size(),
"ConvertFragmentResultToMPQCData() - the number of extracted data differs from the number of results.");
}
/** Print MPQCData from received results.
*
* @param results results with ids to associate with fragment number
* @param fragmentData MPQCData resulting from the jobs
* @param KeySetFilename filename with keysets to associate forces correctly
* @param NoAtoms total number of atoms
*/
bool printReceivedMPQCResults(
const std::vector &results,
const std::vector &fragmentData,
const std::string &KeySetFilename,
size_t NoAtoms)
{
EnergyMatrix Energy;
EnergyMatrix EnergyFragments;
ForceMatrix Force;
ForceMatrix ForceFragments;
// align fragments
std::map< JobId_t, size_t > MatrixNrLookup;
size_t FragmentCounter = 0;
{
// bring ids in order ...
typedef std::map< JobId_t, FragmentResult::ptr> IdResultMap_t;
IdResultMap_t IdResultMap;
for (std::vector::const_iterator iter = results.begin();
iter != results.end(); ++iter) {
#ifndef NDEBUG
std::pair< IdResultMap_t::iterator, bool> inserter =
#endif
IdResultMap.insert( make_pair((*iter)->getId(), *iter) );
ASSERT( inserter.second,
"ExtractMPQCDataFromResults() - two results have same id "
+toString((*iter)->getId())+".");
}
// ... and fill lookup
for(IdResultMap_t::const_iterator iter = IdResultMap.begin();
iter != IdResultMap.end(); ++iter)
MatrixNrLookup.insert( make_pair(iter->first, FragmentCounter++) );
}
LOG(1, "INFO: There are " << FragmentCounter << " fragments.");
std::vector::const_iterator dataiter = fragmentData.begin();
std::vector::const_iterator resultiter = results.begin();
for (; dataiter != fragmentData.end(); ++dataiter, ++resultiter) {
const MPQCData &extractedData = *dataiter;
// place results into EnergyMatrix ...
{
MatrixContainer::MatrixArray matrix;
matrix.resize(1);
matrix[0].resize(1, extractedData.energies.total);
if (!Energy.AddMatrix(
std::string("MPQCJob ")+toString((*resultiter)->getId()),
matrix,
MatrixNrLookup[(*resultiter)->getId()])) {
ELOG(1, "Adding energy matrix failed.");
return false;
}
}
// ... and ForceMatrix (with two empty columns in front)
{
MatrixContainer::MatrixArray matrix;
const size_t rows = extractedData.forces.size();
matrix.resize(rows);
for (size_t i=0;igetId()),
matrix,
MatrixNrLookup[(*resultiter)->getId()])) {
ELOG(1, "Adding force matrix failed.");
return false;
}
}
}
// add one more matrix (not required for energy)
MatrixContainer::MatrixArray matrix;
matrix.resize(1);
matrix[0].resize(1, 0.);
if (!Energy.AddMatrix(std::string("MPQCJob total"), matrix, FragmentCounter))
return false;
// but for energy because we need to know total number of atoms
matrix.resize(NoAtoms);
for (size_t i = 0; i< NoAtoms; ++i)
matrix[i].resize(2+NDIM, 0.);
if (!Force.AddMatrix(std::string("MPQCJob total"), matrix, FragmentCounter))
return false;
// initialise indices
KeySetsContainer KeySet;
if (!Energy.InitialiseIndices()) return false;
if (!Force.ParseIndices(KeySetFilename.c_str())) return false;
if (!KeySet.ParseKeySets(KeySetFilename.c_str(), Force.RowCounter, Force.MatrixCounter)) return false;
/// prepare for OrthogonalSummation
// gather all present indices in AllIndices
IndexSet::ptr AllIndices(new IndexSet);
for (KeySetsContainer::ArrayOfIntVectors::const_iterator iter = KeySet.KeySets.begin();
iter != KeySet.KeySets.end(); ++iter)
for(KeySetsContainer::IntVector::const_iterator keyiter = (*iter).begin();
keyiter != (*iter).end(); ++keyiter) {
if (*keyiter != -1)
(*AllIndices) += *keyiter;
}
LOG(1, "INFO: AllIndices is " << AllIndices << ".");
// create container with all keysets
IndexSetContainer::ptr container(new IndexSetContainer(AllIndices));
for (KeySetsContainer::ArrayOfIntVectors::const_iterator iter = KeySet.KeySets.begin();
iter != KeySet.KeySets.end(); ++iter) {
IndexSet tempset;
for(KeySetsContainer::IntVector::const_iterator keyiter = (*iter).begin();
keyiter != (*iter).end(); ++keyiter)
if (*keyiter != -1)
tempset += *keyiter;
container->insert(tempset);
}
// create the map of all keysets
SubsetMap::ptr subsetmap(new SubsetMap(*container));
/// convert all MPQCData to MPQCDataMap_t
{
// energy_t
std::vector MPQCData_Energy_fused;
MPQCData_Energy_fused.reserve(fragmentData.size());
for(std::vector::const_iterator dataiter = fragmentData.begin();
dataiter != fragmentData.end(); ++dataiter) {
const MPQCData &extractedData = *dataiter;
LOG(2, "DEBUG: Current extracted Data is " << extractedData << ".");
MPQCDataEnergyMap_t instance;
boost::fusion::at_key(instance) = extractedData.energies.total;
boost::fusion::at_key(instance) = extractedData.energies.nuclear_repulsion;
boost::fusion::at_key(instance) = extractedData.energies.electron_coulomb;
boost::fusion::at_key(instance) = extractedData.energies.electron_exchange;
boost::fusion::at_key(instance) = extractedData.energies.correlation;
boost::fusion::at_key(instance) = extractedData.energies.overlap;
boost::fusion::at_key(instance) = extractedData.energies.kinetic;
boost::fusion::at_key(instance) = extractedData.energies.hcore;
boost::fusion::at_key(instance) = extractedData.energies.eigenvalues;
MPQCData_Energy_fused.push_back(instance);
}
// forces
const IndexSetContainer::Container_t &indices = container->getContainer();
ASSERT( indices.size() == fragmentData.size(),
"FragmentationAutomationAction::performCall() - indices and fragmentData differ in size.");
std::vector MPQCData_Force_fused;
MPQCData_Force_fused.reserve(fragmentData.size());
std::vector::const_iterator dataiter = fragmentData.begin();
IndexSetContainer::Container_t::const_iterator indexiter = indices.begin();
for(;dataiter != fragmentData.end(); ++dataiter, ++indexiter) {
const MPQCData &extractedData = *dataiter;
LOG(2, "DEBUG: Current extracted Data is " << extractedData << ".");
MPQCDataForceMap_t instance;
boost::fusion::at_key(instance) =
IndexedVectors(**indexiter, extractedData.forces);
MPQCData_Force_fused.push_back(instance);
}
// sampled_grid
std::vector MPQCData_Grid_fused;
MPQCData_Grid_fused.reserve(fragmentData.size());
for(std::vector::const_iterator dataiter = fragmentData.begin();
dataiter != fragmentData.end(); ++dataiter) {
const MPQCData &extractedData = *dataiter;
LOG(2, "DEBUG: Current extracted Data is " << extractedData << ".");
MPQCDataGridMap_t instance;
boost::fusion::at_key(instance) = extractedData.sampled_grid;
MPQCData_Grid_fused.push_back(instance);
}
// times
std::vector MPQCData_Time_fused;
MPQCData_Time_fused.reserve(fragmentData.size());
for(std::vector::const_iterator dataiter = fragmentData.begin();
dataiter != fragmentData.end(); ++dataiter) {
const MPQCData &extractedData = *dataiter;
LOG(2, "DEBUG: Current extracted Data is " << extractedData << ".");
MPQCDataTimeMap_t instance;
boost::fusion::at_key(instance) = extractedData.times.walltime;
boost::fusion::at_key(instance) = extractedData.times.cputime;
boost::fusion::at_key(instance) = extractedData.times.flops;
MPQCData_Time_fused.push_back(instance);
}
// create a vector of all job ids
std::vector jobids(results.size(), JobId::IllegalJob);
std::transform(results.begin(), results.end(), jobids.begin(),
boost::bind(&FragmentResult::getId,
boost::bind(&FragmentResult::ptr::operator->, _1)));
// sum up and print energies
boost::mpl::for_each(
printOrthogonalSum(
subsetmap,
MPQCData_Energy_fused,
jobids,
container->getContainer(),
MatrixNrLookup)
);
// sum up and print forces
boost::mpl::for_each(
printOrthogonalSum(
subsetmap,
MPQCData_Force_fused,
jobids,
container->getContainer(),
MatrixNrLookup)
);
// sum up and print grids
boost::mpl::for_each(
printOrthogonalSum(
subsetmap,
MPQCData_Grid_fused,
jobids,
container->getContainer(),
MatrixNrLookup)
);
// sum up and print times
boost::mpl::for_each(
printSum(
subsetmap,
MPQCData_Time_fused,
jobids,
container->getContainer(),
MatrixNrLookup)
);
}
// combine all found data
if (!KeySet.ParseManyBodyTerms()) return false;
if (!EnergyFragments.AllocateMatrix(Energy.Header, Energy.MatrixCounter, Energy.RowCounter, Energy.ColumnCounter)) return false;
if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return false;
if(!Energy.SetLastMatrix(0., 0)) return false;
if(!Force.SetLastMatrix(0., 2)) return false;
for (int BondOrder=0;BondOrder &jobfiles)
{
std::vector jobs;
for (std::vector< boost::filesystem::path >::const_iterator iter = jobfiles.begin();
iter != jobfiles .end(); ++iter) {
const std::string &filename = (*iter).string();
if (boost::filesystem::exists(filename)) {
const JobId_t next_id = controller.getAvailableId();
LOG(1, "INFO: Creating MPQCCommandJob with filename'"
+filename+"', and id "+toString(next_id)+".");
parsejob(jobs, params.executable.get().string(), filename, next_id);
} else {
ELOG(1, "Fragment job "+filename+" does not exist.");
return false;
}
}
controller.addJobs(jobs);
controller.sendJobs(params.host.get(), params.port.get());
return true;
}
#ifdef HAVE_VMG
bool createLongRangeJobs(
FragmentController &controller,
const FragmentationFragmentationAutomationAction::FragmentationFragmentationAutomationParameters ¶ms,
const std::vector &fragmentData)
{
std::vector jobs;
// add one job for each fragment as the short-range correction which we need
// to subtract from the obtained full potential to get the long-range part only
for (std::vector::const_iterator iter = fragmentData.begin();
iter != fragmentData.end(); ++iter) {
const JobId_t next_id = controller.getAvailableId();
LOG(1, "INFO: Creating VMGJob.");
FragmentJob::ptr testJob(
new VMGJob(next_id, iter->sampled_grid) );
jobs.push_back(testJob);
}
// add one more job for the full calculation
// TODO: Here, we actually have to combine all the other sampled_grids
{
const int level = LEVEL;
const int GRID = pow(2, level);
std::vector full_sample(GRID*GRID*GRID, 0.);
double begin[NDIM] = { 0., 0., 0. };
const RealSpaceMatrix& M = World::getInstance().getDomain().getM();
const double size = M.at(0,0);
ASSERT( M.determinant() == size*size*size,
"createLongRangeJobs() - current domain matrix "+toString(M)+" is not cubic.");
const SamplingGrid full_sampled_grid(begin, size, level, full_sample);
const JobId_t next_id = controller.getAvailableId();
FragmentJob::ptr testJob(
new VMGJob(next_id, full_sampled_grid) );
jobs.push_back(testJob);
}
// then send jobs to controller
controller.addJobs(jobs);
controller.sendJobs(params.host.get(), params.port.get());
return true;
}
#endif
void WaitforResults(
boost::asio::io_service &io_service,
FragmentController &controller,
const FragmentationFragmentationAutomationAction::FragmentationFragmentationAutomationParameters ¶ms,
const size_t NoExpectedResults
)
{
size_t NoCalculatedResults = 0;
while (NoCalculatedResults != NoExpectedResults) {
// wait a bit
boost::asio::deadline_timer timer(io_service);
timer.expires_from_now(boost::posix_time::milliseconds(500));
timer.wait();
// then request status
controller.checkResults(params.host.get(), params.port.get());
RunService(io_service, "Checking on results");
const std::pair JobStatus = controller.getJobStatus();
LOG(1, "INFO: #" << JobStatus.first << " are waiting in the queue and #" << JobStatus.second << " jobs are calculated so far.");
NoCalculatedResults = JobStatus.second;
}
}
Action::state_ptr FragmentationFragmentationAutomationAction::performCall() {
boost::asio::io_service io_service;
FragmentController controller(io_service);
// TODO: Have io_service run in second thread and merge with current again eventually
// Phase One: obtain ids
std::vector< boost::filesystem::path > jobfiles = params.jobfiles.get();
requestIds(controller, params, jobfiles.size());
RunService(io_service, "Requesting ids");
// Phase Two: create and add MPQCJobs
if (!createJobsFromFiles(controller, params, jobfiles))
return Action::failure;
RunService(io_service, "Adding MPQCJobs");
// Phase Three: calculate result
WaitforResults(io_service, controller, params, jobfiles.size());
controller.receiveResults(params.host.get(), params.port.get());
RunService(io_service, "Requesting short-range results");
std::vector MPQCresults = controller.getReceivedResults();
std::vector fragmentData;
ConvertFragmentResultToMPQCData(MPQCresults, fragmentData);
// print intermediate short-range results
{
LOG(1, "INFO: Parsing fragment files from " << params.path.get() << ".");
printReceivedMPQCResults(
MPQCresults,
fragmentData,
params.path.get(),
getNoAtomsFromAdjacencyFile(params.path.get()));
}
#ifdef HAVE_VMG
if (params.DoLongrange.get()) {
// Phase Four: obtain more ids
requestIds(controller, params, fragmentData.size()+1);
RunService(io_service, "Requesting ids");
// Phase Five: create VMGJobs
if (!createLongRangeJobs(controller, params, fragmentData))
return Action::failure;
RunService(io_service, "Adding VMGJobs");
// Phase Six: calculate result
WaitforResults(io_service, controller, params, fragmentData.size()+1);
controller.receiveResults(params.host.get(), params.port.get());
RunService(io_service, "Requesting long-range results");
std::vector VMGresults = controller.getReceivedResults();
ASSERT( MPQCresults.size()+1 == VMGresults.size(),
"FragmentationFragmentationAutomationAction::performCall() - number of MPQCresultd and VMGresults don't match.");
// Final phase: print result
{
LOG(1, "INFO: Parsing fragment files from " << params.path.get() << ".");
printReceivedMPQCResults(
MPQCresults,
fragmentData,
params.path.get(),
getNoAtomsFromAdjacencyFile(params.path.get()));
}
}
#endif
size_t Exitflag = controller.getExitflag();
return (Exitflag == 0) ? Action::success : Action::failure;
}
Action::state_ptr FragmentationFragmentationAutomationAction::performUndo(Action::state_ptr _state) {
return Action::success;
}
Action::state_ptr FragmentationFragmentationAutomationAction::performRedo(Action::state_ptr _state){
return Action::success;
}
bool FragmentationFragmentationAutomationAction::canUndo() {
return false;
}
bool FragmentationFragmentationAutomationAction::shouldUndo() {
return false;
}
/** =========== end of function ====================== */