source: src/Actions/FragmentationAction/FragmentationAutomationAction.cpp@ 635645

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since 635645 was 635645, checked in by Frederik Heber <heber@…>, 12 years ago

Translated convertMPQCDatato...() into specialized template functions convertMPQCDataTo().

  • Property mode set to 100644
File size: 28.8 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
5 *
6 *
7 * This file is part of MoleCuilder.
8 *
9 * MoleCuilder is free software: you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation, either version 2 of the License, or
12 * (at your option) any later version.
13 *
14 * MoleCuilder is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
21 */
22
23/*
24 * FragmentationAutomationAction.cpp
25 *
26 * Created on: May 18, 2012
27 * Author: heber
28 */
29
30// include config.h
31#ifdef HAVE_CONFIG_H
32#include <config.h>
33#endif
34
35#include <boost/archive/text_iarchive.hpp>
36// boost asio needs specific operator new
37#include <boost/asio.hpp>
38
39#include "CodePatterns/MemDebug.hpp"
40
41#include <boost/mpl/remove.hpp>
42
43#include "CodePatterns/Assert.hpp"
44#include "CodePatterns/Info.hpp"
45#include "CodePatterns/Log.hpp"
46#include "JobMarket/Controller/FragmentController.hpp"
47#include "JobMarket/Jobs/FragmentJob.hpp"
48
49#include "Atom/atom.hpp"
50#include "Box.hpp"
51#include "Fragmentation/EnergyMatrix.hpp"
52#include "Fragmentation/ForceMatrix.hpp"
53#include "Fragmentation/Fragmentation.hpp"
54#include "Fragmentation/SetValues/Histogram.hpp"
55#include "Fragmentation/SetValues/IndexedVectors.hpp"
56#include "Fragmentation/HydrogenSaturation_enum.hpp"
57#include "Fragmentation/KeySet.hpp"
58#include "Fragmentation/KeySetsContainer.hpp"
59#include "Fragmentation/Summation/AllLevelOrthogonalSummator.hpp"
60#include "Fragmentation/Summation/AllLevelSummator.hpp"
61#include "Fragmentation/Summation/writeTable.hpp"
62#include "Graph/DepthFirstSearchAnalysis.hpp"
63#include "Helpers/defs.hpp"
64#include "Jobs/MPQCJob.hpp"
65#include "Jobs/MPQCData.hpp"
66#include "Jobs/MPQCData_printKeyNames.hpp"
67#include "LinearAlgebra/RealSpaceMatrix.hpp"
68#ifdef HAVE_VMG
69#include "Jobs/VMGJob.hpp"
70#endif
71#include "molecule.hpp"
72#include "World.hpp"
73
74#include <iostream>
75#include <string>
76#include <vector>
77
78#include <boost/mpl/for_each.hpp>
79
80#include "Actions/FragmentationAction/FragmentationAutomationAction.hpp"
81
82using namespace MoleCuilder;
83
84// and construct the stuff
85#include "FragmentationAutomationAction.def"
86#include "Action_impl_pre.hpp"
87/** =========== define the function ====================== */
88
89class controller_AddOn;
90
91// needs to be defined for using the FragmentController
92controller_AddOn *getAddOn()
93{
94 return NULL;
95}
96
97const int LEVEL = 5;
98
99/** Creates a MPQCCommandJob with argument \a filename.
100 *
101 * @param jobs created job is added to this vector
102 * @param command mpqc command to execute
103 * @param filename filename being argument to job
104 * @param nextid id for this job
105 */
106void parsejob(
107 std::vector<FragmentJob::ptr> &jobs,
108 const std::string &command,
109 const std::string &filename,
110 const JobId_t nextid)
111{
112 std::ifstream file;
113 file.open(filename.c_str());
114 ASSERT( file.good(), "parsejob() - file "+filename+" does not exist.");
115 std::string output((std::istreambuf_iterator<char>(file)),
116 std::istreambuf_iterator<char>());
117 double begin[NDIM] = { 0., 0., 0. };
118 const RealSpaceMatrix& M = World::getInstance().getDomain().getM();
119 const double size = M.at(0,0);
120 ASSERT( M.determinant() == size*size*size,
121 "parsejob() - current domain matrix "+toString(M)+" is not cubic.");
122 const int level = LEVEL;
123 FragmentJob::ptr testJob( new MPQCJob(nextid, output, begin, size, level) );
124 jobs.push_back(testJob);
125 file.close();
126 LOG(1, "INFO: Added MPQCCommandJob from file "+filename+".");
127}
128
129/** Helper function to get number of atoms somehow.
130 *
131 * Here, we just parse the number of lines in the adjacency file as
132 * it should correspond to the number of atoms, except when some atoms
133 * are not bonded, but then fragmentation makes no sense.
134 *
135 * @param path path to the adjacency file
136 */
137size_t getNoAtomsFromAdjacencyFile(const std::string &path)
138{
139 size_t NoAtoms = 0;
140
141 // parse in special file to get atom count (from line count)
142 std::string filename(path);
143 filename += FRAGMENTPREFIX;
144 filename += ADJACENCYFILE;
145 std::ifstream adjacency(filename.c_str());
146 if (adjacency.fail()) {
147 LOG(0, endl << "getNoAtomsFromAdjacencyFile() - Unable to open " << filename << ", is the directory correct?");
148 return false;
149 }
150 std::string buffer;
151 while (getline(adjacency, buffer))
152 NoAtoms++;
153 LOG(1, "INFO: There are " << NoAtoms << " atoms.");
154
155 return NoAtoms;
156}
157
158/** Extracts MPQCData from received vector of FragmentResults.
159 *
160 * @param results results to extract MPQCData from
161 * @param fragmentData on return array filled with extracted MPQCData
162 */
163void ConvertFragmentResultToMPQCData(
164 const std::vector<FragmentResult::ptr> &results,
165 std::vector<MPQCData> &fragmentData)
166{
167 // extract results
168 fragmentData.clear();
169 fragmentData.reserve(results.size());
170
171 LOG(2, "DEBUG: Parsing now through " << results.size() << " results.");
172 for (std::vector<FragmentResult::ptr>::const_iterator iter = results.begin();
173 iter != results.end(); ++iter) {
174 //LOG(1, "RESULT: job #"+toString((*iter)->getId())+": "+toString((*iter)->result));
175 MPQCData extractedData;
176 std::stringstream inputstream((*iter)->result);
177 LOG(2, "DEBUG: First 50 characters FragmentResult's string: "+(*iter)->result.substr(0, 50));
178 boost::archive::text_iarchive ia(inputstream);
179 ia >> extractedData;
180 LOG(1, "INFO: extracted data is " << extractedData << ".");
181 fragmentData.push_back(extractedData);
182 }
183
184 ASSERT( results.size() == fragmentData.size(),
185 "ConvertFragmentResultToMPQCData() - the number of extracted data differs from the number of results.");
186}
187
188/** Creates a lookup from FragmentJob::id to the true fragment number.
189 *
190 * @param results result with job ids
191 * @param MatrixNrLookup Lookup up-map, filled on return
192 * @param FragmentCounter total number of fragments on return
193 */
194void createMatrixNrLookup(
195 const std::vector<FragmentResult::ptr> &results,
196 std::map< JobId_t, size_t > &MatrixNrLookup,
197 size_t &FragmentCounter)
198{
199 // align fragments
200 MatrixNrLookup.clear();
201 FragmentCounter = 0;
202 {
203 // bring ids in order ...
204 typedef std::map< JobId_t, FragmentResult::ptr> IdResultMap_t;
205 IdResultMap_t IdResultMap;
206 for (std::vector<FragmentResult::ptr>::const_iterator iter = results.begin();
207 iter != results.end(); ++iter) {
208 #ifndef NDEBUG
209 std::pair< IdResultMap_t::iterator, bool> inserter =
210 #endif
211 IdResultMap.insert( make_pair((*iter)->getId(), *iter) );
212 ASSERT( inserter.second,
213 "ExtractMPQCDataFromResults() - two results have same id "
214 +toString((*iter)->getId())+".");
215 }
216 // ... and fill lookup
217 for(IdResultMap_t::const_iterator iter = IdResultMap.begin();
218 iter != IdResultMap.end(); ++iter)
219 MatrixNrLookup.insert( make_pair(iter->first, FragmentCounter++) );
220 }
221 LOG(1, "INFO: There are " << FragmentCounter << " fragments.");
222}
223
224/** Place results from FragmentResult into EnergyMatrix and ForceMatrix.
225 *
226 * @param results results with ids to associate with fragment number
227 * @param fragmentData MPQCData resulting from the jobs
228 * @param MatrixNrLookup Lookup up-map from job id to fragment number
229 * @param FragmentCounter total number of fragments
230 * @param NoAtoms total number of atoms
231 * @param Energy energy matrix to be filled on return
232 * @param Force force matrix to be filled on return
233 * @return true - everything ok, false - else
234 */
235bool putResultsintoMatrices(
236 const std::vector<FragmentResult::ptr> &results,
237 const std::vector<MPQCData> &fragmentData,
238 std::map< JobId_t, size_t > &MatrixNrLookup,
239 const size_t FragmentCounter,
240 const size_t NoAtoms,
241 EnergyMatrix &Energy,
242 ForceMatrix &Force)
243{
244 ASSERT( results.size() == fragmentData.size(),
245 "printReceivedMPQCResults() - results and fragmentData differ in size.");
246 std::vector<MPQCData>::const_iterator dataiter = fragmentData.begin();
247 std::vector<FragmentResult::ptr>::const_iterator resultiter = results.begin();
248 for (; dataiter != fragmentData.end(); ++dataiter, ++resultiter) {
249 const MPQCData &extractedData = *dataiter;
250 // place results into EnergyMatrix ...
251 {
252 MatrixContainer::MatrixArray matrix;
253 matrix.resize(1);
254 matrix[0].resize(1, extractedData.energies.total);
255 if (!Energy.AddMatrix(
256 std::string("MPQCJob ")+toString((*resultiter)->getId()),
257 matrix,
258 MatrixNrLookup[(*resultiter)->getId()])) {
259 ELOG(1, "Adding energy matrix failed.");
260 return false;
261 }
262 }
263 // ... and ForceMatrix (with two empty columns in front)
264 {
265 MatrixContainer::MatrixArray matrix;
266 const size_t rows = extractedData.forces.size();
267 matrix.resize(rows);
268 for (size_t i=0;i<rows;++i) {
269 const size_t columns = 2+extractedData.forces[i].size();
270 matrix[i].resize(columns, 0.);
271 // for (size_t j=0;j<2;++j)
272 // matrix[i][j] = 0.;
273 for (size_t j=2;j<columns;++j)
274 matrix[i][j] = extractedData.forces[i][j-2];
275 }
276 if (!Force.AddMatrix(
277 std::string("MPQCJob ")+toString((*resultiter)->getId()),
278 matrix,
279 MatrixNrLookup[(*resultiter)->getId()])) {
280 ELOG(1, "Adding force matrix failed.");
281 return false;
282 }
283 }
284 }
285 // add one more matrix (not required for energy)
286 MatrixContainer::MatrixArray matrix;
287 matrix.resize(1);
288 matrix[0].resize(1, 0.);
289 if (!Energy.AddMatrix(std::string("MPQCJob total"), matrix, FragmentCounter))
290 return false;
291 // but for energy because we need to know total number of atoms
292 matrix.resize(NoAtoms);
293 for (size_t i = 0; i< NoAtoms; ++i)
294 matrix[i].resize(2+NDIM, 0.);
295 if (!Force.AddMatrix(std::string("MPQCJob total"), matrix, FragmentCounter))
296 return false;
297
298 return true;
299}
300
301template <typename T>
302void convertMPQCDataTo(
303 const std::vector<MPQCData> &fragmentData,
304 std::vector<T> &MPQCData_fused)
305{
306 MPQCData_fused.clear();
307}
308
309template <>
310void convertMPQCDataTo<MPQCDataEnergyMap_t>(
311 const std::vector<MPQCData> &fragmentData,
312 std::vector<MPQCDataEnergyMap_t> &MPQCData_Energy_fused)
313{
314 // energy_t
315 MPQCData_Energy_fused.clear();
316 MPQCData_Energy_fused.reserve(fragmentData.size());
317 for(std::vector<MPQCData>::const_iterator dataiter = fragmentData.begin();
318 dataiter != fragmentData.end(); ++dataiter) {
319 const MPQCData &extractedData = *dataiter;
320 LOG(2, "DEBUG: Current extracted Data is " << extractedData << ".");
321 MPQCDataEnergyMap_t instance;
322 boost::fusion::at_key<MPQCDataFused::energy_total>(instance) = extractedData.energies.total;
323 boost::fusion::at_key<MPQCDataFused::energy_nuclear_repulsion>(instance) = extractedData.energies.nuclear_repulsion;
324 boost::fusion::at_key<MPQCDataFused::energy_electron_coulomb>(instance) = extractedData.energies.electron_coulomb;
325 boost::fusion::at_key<MPQCDataFused::energy_electron_exchange>(instance) = extractedData.energies.electron_exchange;
326 boost::fusion::at_key<MPQCDataFused::energy_correlation>(instance) = extractedData.energies.correlation;
327 boost::fusion::at_key<MPQCDataFused::energy_overlap>(instance) = extractedData.energies.overlap;
328 boost::fusion::at_key<MPQCDataFused::energy_kinetic>(instance) = extractedData.energies.kinetic;
329 boost::fusion::at_key<MPQCDataFused::energy_hcore>(instance) = extractedData.energies.hcore;
330 boost::fusion::at_key<MPQCDataFused::energy_eigenvalues>(instance) = extractedData.energies.eigenvalues;
331 MPQCData_Energy_fused.push_back(instance);
332 }
333}
334
335
336void convertMPQCDatatoForceMap(
337 const std::vector<MPQCData> &fragmentData,
338 const KeySetsContainer &ForceKeySet,
339 std::vector<MPQCDataForceMap_t> &MPQCData_Force_fused)
340{
341 // forces
342 ASSERT( ForceKeySet.KeySets.size() == fragmentData.size(),
343 "FragmentationAutomationAction::performCall() - indices and fragmentData differ in size.");
344 MPQCData_Force_fused.clear();
345 MPQCData_Force_fused.reserve(fragmentData.size());
346 std::vector<MPQCData>::const_iterator dataiter = fragmentData.begin();
347 KeySetsContainer::ArrayOfIntVectors::const_iterator arrayiter = ForceKeySet.KeySets.begin();
348 for(;dataiter != fragmentData.end(); ++dataiter, ++arrayiter) {
349 const MPQCData &extractedData = *dataiter;
350 LOG(2, "DEBUG: Current extracted Data is " << extractedData << ".");
351 MPQCDataForceMap_t instance;
352 // must convert int to index_t
353 IndexedVectors::indices_t indices(arrayiter->begin(), arrayiter->end());
354 boost::fusion::at_key<MPQCDataFused::forces>(instance) =
355 IndexedVectors(indices, extractedData.forces);
356 MPQCData_Force_fused.push_back(instance);
357 }
358}
359
360template <>
361void convertMPQCDataTo<MPQCDataGridMap_t>(
362 const std::vector<MPQCData> &fragmentData,
363 std::vector<MPQCDataGridMap_t> &MPQCData_Grid_fused)
364{
365 // sampled_grid
366 MPQCData_Grid_fused.clear();
367 MPQCData_Grid_fused.reserve(fragmentData.size());
368 for(std::vector<MPQCData>::const_iterator dataiter = fragmentData.begin();
369 dataiter != fragmentData.end(); ++dataiter) {
370 const MPQCData &extractedData = *dataiter;
371 LOG(2, "DEBUG: Current extracted Data is " << extractedData << ".");
372 MPQCDataGridMap_t instance;
373 boost::fusion::at_key<MPQCDataFused::sampled_grid>(instance) = extractedData.sampled_grid;
374 MPQCData_Grid_fused.push_back(instance);
375 }
376}
377
378template <>
379void convertMPQCDataTo<MPQCDataTimeMap_t>(
380 const std::vector<MPQCData> &fragmentData,
381 std::vector<MPQCDataTimeMap_t> &MPQCData_Time_fused)
382{
383 // times
384 MPQCData_Time_fused.clear();
385 MPQCData_Time_fused.reserve(fragmentData.size());
386 for(std::vector<MPQCData>::const_iterator dataiter = fragmentData.begin();
387 dataiter != fragmentData.end(); ++dataiter) {
388 const MPQCData &extractedData = *dataiter;
389 LOG(2, "DEBUG: Current extracted Data is " << extractedData << ".");
390 MPQCDataTimeMap_t instance;
391 boost::fusion::at_key<MPQCDataFused::times_walltime>(instance) = extractedData.times.walltime;
392 boost::fusion::at_key<MPQCDataFused::times_cputime>(instance) = extractedData.times.cputime;
393 boost::fusion::at_key<MPQCDataFused::times_flops>(instance) = extractedData.times.flops;
394 MPQCData_Time_fused.push_back(instance);
395 }
396}
397
398/** Print MPQCData from received results.
399 *
400 * @param results results with ids to associate with fragment number
401 * @param fragmentData MPQCData resulting from the jobs
402 * @param KeySetFilename filename with keysets to associate forces correctly
403 * @param NoAtoms total number of atoms
404 */
405bool printReceivedMPQCResults(
406 const std::vector<FragmentResult::ptr> &results,
407 const std::vector<MPQCData> &fragmentData,
408 const std::string &KeySetFilename,
409 size_t NoAtoms)
410{
411 // create lookup from job nr to fragment number
412 std::map< JobId_t, size_t > MatrixNrLookup;
413 size_t FragmentCounter = 0;
414 createMatrixNrLookup(results, MatrixNrLookup, FragmentCounter);
415
416 // place results into maps
417 EnergyMatrix Energy;
418 ForceMatrix Force;
419 if (!putResultsintoMatrices(results, fragmentData, MatrixNrLookup, FragmentCounter, NoAtoms, Energy, Force))
420 return false;
421
422 // initialise keysets
423 KeySetsContainer KeySet;
424 KeySetsContainer ForceKeySet;
425 if (!Energy.InitialiseIndices()) return false;
426
427 if (!Force.ParseIndices(KeySetFilename.c_str())) return false;
428
429 {
430 // else needs keysets without hydrogens
431 std::stringstream filename;
432 filename << FRAGMENTPREFIX << KEYSETFILE;
433 if (!KeySet.ParseKeySets(KeySetFilename, filename.str(), FragmentCounter)) return false;
434 }
435
436 {
437 // forces need keysets including hydrogens
438 std::stringstream filename;
439 filename << FRAGMENTPREFIX << FORCESFILE;
440 if (!ForceKeySet.ParseKeySets(KeySetFilename, filename.str(), FragmentCounter)) return false;
441 }
442
443 /// prepare for OrthogonalSummation
444
445 // convert KeySetContainer to IndexSetContainer
446 IndexSetContainer::ptr container(new IndexSetContainer(KeySet));
447 // create the map of all keysets
448 SubsetMap::ptr subsetmap(new SubsetMap(*container));
449
450 // create a vector of all job ids
451 std::vector<JobId_t> jobids(results.size(), JobId::IllegalJob);
452 std::transform(results.begin(), results.end(), jobids.begin(),
453 boost::bind(&FragmentResult::getId,
454 boost::bind(&FragmentResult::ptr::operator->, _1)));
455
456 /// convert all MPQCData to MPQCDataMap_t
457 {
458 ASSERT( ForceKeySet.KeySets.size() == fragmentData.size(),
459 "FragmentationAutomationAction::performCall() - ForceKeySet's KeySets and fragmentData differ in size.");
460
461 std::vector<MPQCDataEnergyMap_t> MPQCData_Energy_fused;
462 convertMPQCDataTo<MPQCDataEnergyMap_t>(fragmentData, MPQCData_Energy_fused);
463 std::vector<MPQCDataForceMap_t> MPQCData_Force_fused;
464 convertMPQCDatatoForceMap(fragmentData, ForceKeySet, MPQCData_Force_fused);
465 std::vector<MPQCDataGridMap_t> MPQCData_Grid_fused;
466 convertMPQCDataTo<MPQCDataGridMap_t>(fragmentData, MPQCData_Grid_fused);
467 std::vector<MPQCDataTimeMap_t> MPQCData_Time_fused;
468 convertMPQCDataTo<MPQCDataTimeMap_t>(fragmentData, MPQCData_Time_fused);
469
470 // create summation instances
471 std::vector<MPQCDataEnergyMap_t> Result_Energy_fused(subsetmap->getMaximumSubsetLevel());
472 AllLevelOrthogonalSummator<MPQCDataEnergyMap_t> energySummer(
473 subsetmap,
474 MPQCData_Energy_fused,
475 jobids,
476 container->getContainer(),
477 MatrixNrLookup,
478 Result_Energy_fused);
479 std::vector<MPQCDataForceMap_t> Result_Force_fused(subsetmap->getMaximumSubsetLevel());
480 AllLevelOrthogonalSummator<MPQCDataForceMap_t> forceSummer(
481 subsetmap,
482 MPQCData_Force_fused,
483 jobids,
484 container->getContainer(),
485 MatrixNrLookup,
486 Result_Force_fused);
487 std::vector<MPQCDataGridMap_t> Result_Grid_fused(subsetmap->getMaximumSubsetLevel());
488 AllLevelOrthogonalSummator<MPQCDataGridMap_t> gridSummer(
489 subsetmap,
490 MPQCData_Grid_fused,
491 jobids,
492 container->getContainer(),
493 MatrixNrLookup,
494 Result_Grid_fused);
495 std::vector<MPQCDataTimeMap_t> Result_Time_fused(subsetmap->getMaximumSubsetLevel());
496 AllLevelSummator<MPQCDataTimeMap_t> timeSummer(
497 subsetmap,
498 MPQCData_Time_fused,
499 jobids,
500 container->getContainer(),
501 MatrixNrLookup,
502 Result_Time_fused);
503
504 // sum up
505 boost::mpl::for_each<MPQCDataEnergyVector_t>(boost::ref(energySummer));
506 boost::mpl::for_each<MPQCDataForceVector_t>(boost::ref(forceSummer));
507 boost::mpl::for_each<MPQCDataGridVector_t>(boost::ref(gridSummer));
508 boost::mpl::for_each<MPQCDataTimeVector_t>(boost::ref(timeSummer));
509
510 // print tables (without eigenvalues, they go extra)
511 const size_t MaxLevel = subsetmap->getMaximumSubsetLevel();
512 typedef boost::mpl::remove<MPQCDataEnergyVector_t, MPQCDataFused::energy_eigenvalues>::type MPQCDataEnergyVector_noeigenvalues_t;
513 const std::string energyresult =
514 writeTable<MPQCDataEnergyMap_t, MPQCDataEnergyVector_noeigenvalues_t >()(
515 Result_Energy_fused, MaxLevel);
516 LOG(0, "Energy table is \n" << energyresult);
517 const std::string eigenvalueresult;
518
519 LOG(0, "Eigenvalue table is \n" << eigenvalueresult);
520 const std::string forceresult =
521 writeTable<MPQCDataForceMap_t, MPQCDataForceVector_t>()(
522 Result_Force_fused, MaxLevel);
523 LOG(0, "Force table is \n" << forceresult);
524 // we don't want to print grid to a table
525 // print times (without flops for now)
526 typedef boost::mpl::remove<MPQCDataTimeVector_t, MPQCDataFused::times_flops>::type MPQCDataTimeVector_noflops_t;
527 const std::string timesresult =
528 writeTable<MPQCDataTimeMap_t, MPQCDataTimeVector_noflops_t >()(
529 Result_Time_fused, MaxLevel);
530 LOG(0, "Times table is \n" << timesresult);
531 }
532
533 // combine all found data
534 if (!KeySet.ParseManyBodyTerms()) return false;
535
536 EnergyMatrix EnergyFragments;
537 ForceMatrix ForceFragments;
538 if (!EnergyFragments.AllocateMatrix(Energy.Header, Energy.MatrixCounter, Energy.RowCounter, Energy.ColumnCounter)) return false;
539 if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return false;
540
541 if(!Energy.SetLastMatrix(0., 0)) return false;
542 if(!Force.SetLastMatrix(0., 2)) return false;
543
544 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
545 // --------- sum up energy --------------------
546 LOG(1, "INFO: Summing energy of order " << BondOrder+1 << " ...");
547 if (!EnergyFragments.SumSubManyBodyTerms(Energy, KeySet, BondOrder)) return false;
548 if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return false;
549
550 // --------- sum up Forces --------------------
551 LOG(1, "INFO: Summing forces of order " << BondOrder+1 << " ...");
552 if (!ForceFragments.SumSubManyBodyTerms(Force, KeySet, BondOrder)) return false;
553 if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder, 1.)) return false;
554 }
555
556 // for debugging print resulting energy and forces
557 LOG(1, "INFO: Resulting energy is " << Energy.Matrix[ FragmentCounter ][0][0]);
558 std::stringstream output;
559 for (int i=0; i< Force.RowCounter[FragmentCounter]; ++i) {
560 for (int j=0; j< Force.ColumnCounter[FragmentCounter]; ++j)
561 output << Force.Matrix[ FragmentCounter ][i][j] << " ";
562 output << "\n";
563 }
564 LOG(1, "INFO: Resulting forces are " << std::endl << output.str());
565
566 return true;
567}
568
569
570void RunService(
571 boost::asio::io_service &io_service,
572 std::string message)
573{
574 message = std::string("io_service: ") + message;
575 io_service.reset();
576 Info info(message.c_str());
577 io_service.run();
578}
579
580void requestIds(
581 FragmentController &controller,
582 const FragmentationFragmentationAutomationAction::FragmentationFragmentationAutomationParameters &params,
583 const size_t numberjobs)
584{
585 controller.requestIds(params.host.get(), params.port.get(), numberjobs);
586}
587
588bool createJobsFromFiles(
589 FragmentController &controller,
590 const FragmentationFragmentationAutomationAction::FragmentationFragmentationAutomationParameters &params,
591 const std::vector< boost::filesystem::path > &jobfiles)
592{
593 std::vector<FragmentJob::ptr> jobs;
594 for (std::vector< boost::filesystem::path >::const_iterator iter = jobfiles.begin();
595 iter != jobfiles .end(); ++iter) {
596 const std::string &filename = (*iter).string();
597 if (boost::filesystem::exists(filename)) {
598 const JobId_t next_id = controller.getAvailableId();
599 LOG(1, "INFO: Creating MPQCCommandJob with filename'"
600 +filename+"', and id "+toString(next_id)+".");
601 parsejob(jobs, params.executable.get().string(), filename, next_id);
602 } else {
603 ELOG(1, "Fragment job "+filename+" does not exist.");
604 return false;
605 }
606 }
607 controller.addJobs(jobs);
608 controller.sendJobs(params.host.get(), params.port.get());
609 return true;
610}
611
612#ifdef HAVE_VMG
613bool createLongRangeJobs(
614 FragmentController &controller,
615 const FragmentationFragmentationAutomationAction::FragmentationFragmentationAutomationParameters &params,
616 const std::vector<MPQCData> &fragmentData)
617{
618 std::vector<FragmentJob::ptr> jobs;
619 // add one job for each fragment as the short-range correction which we need
620 // to subtract from the obtained full potential to get the long-range part only
621 for (std::vector<MPQCData>::const_iterator iter = fragmentData.begin();
622 iter != fragmentData.end(); ++iter) {
623 const JobId_t next_id = controller.getAvailableId();
624 LOG(1, "INFO: Creating VMGJob with " << iter->sampled_grid.sampled_grid.size()
625 << " gridpoints and " << iter->charges.size() << " particle charges.");
626 FragmentJob::ptr testJob(
627 new VMGJob(next_id, iter->sampled_grid, iter->positions, iter->charges) );
628 jobs.push_back(testJob);
629 }
630
631 // add one more job for the full calculation
632 // TODO: Here, we actually have to combine all the other sampled_grids
633 {
634 const int level = LEVEL;
635 const int GRID = pow(2, level);
636 std::vector<double> full_sample(GRID*GRID*GRID, 0.);
637 double begin[NDIM] = { 0., 0., 0. };
638 const RealSpaceMatrix& M = World::getInstance().getDomain().getM();
639 const double size = M.at(0,0);
640 ASSERT( M.determinant() == size*size*size,
641 "createLongRangeJobs() - current domain matrix "+toString(M)+" is not cubic.");
642 const SamplingGrid full_sampled_grid(begin, size, level, full_sample);
643 const std::vector< std::vector<double> > positions;
644 const std::vector<double> charges;
645 const JobId_t next_id = controller.getAvailableId();
646 LOG(1, "INFO: Creating full VMGJob with " << full_sample.size()
647 << " gridpoints and " << charges.size() << " particle charges.");
648 FragmentJob::ptr testJob(
649 new VMGJob(next_id, full_sampled_grid, positions, charges) );
650 jobs.push_back(testJob);
651 }
652
653 // then send jobs to controller
654 controller.addJobs(jobs);
655 controller.sendJobs(params.host.get(), params.port.get());
656 return true;
657}
658#endif
659
660void WaitforResults(
661 boost::asio::io_service &io_service,
662 FragmentController &controller,
663 const FragmentationFragmentationAutomationAction::FragmentationFragmentationAutomationParameters &params,
664 const size_t NoExpectedResults
665 )
666{
667 size_t NoCalculatedResults = 0;
668 while (NoCalculatedResults != NoExpectedResults) {
669 // wait a bit
670 boost::asio::deadline_timer timer(io_service);
671 timer.expires_from_now(boost::posix_time::milliseconds(500));
672 timer.wait();
673 // then request status
674 controller.checkResults(params.host.get(), params.port.get());
675 RunService(io_service, "Checking on results");
676
677 const std::pair<size_t, size_t> JobStatus = controller.getJobStatus();
678 LOG(1, "INFO: #" << JobStatus.first << " are waiting in the queue and #" << JobStatus.second << " jobs are calculated so far.");
679 NoCalculatedResults = JobStatus.second;
680 }
681}
682
683
684Action::state_ptr FragmentationFragmentationAutomationAction::performCall() {
685 boost::asio::io_service io_service;
686 FragmentController controller(io_service);
687
688 // TODO: Have io_service run in second thread and merge with current again eventually
689
690 // Phase One: obtain ids
691 std::vector< boost::filesystem::path > jobfiles = params.jobfiles.get();
692 requestIds(controller, params, jobfiles.size());
693 RunService(io_service, "Requesting ids");
694
695 // Phase Two: create and add MPQCJobs
696 if (!createJobsFromFiles(controller, params, jobfiles))
697 return Action::failure;
698 RunService(io_service, "Adding MPQCJobs");
699
700 // Phase Three: calculate result
701 WaitforResults(io_service, controller, params, jobfiles.size());
702 controller.receiveResults(params.host.get(), params.port.get());
703 RunService(io_service, "Requesting short-range results");
704 std::vector<FragmentResult::ptr> MPQCresults = controller.getReceivedResults();
705 std::vector<MPQCData> fragmentData;
706 ConvertFragmentResultToMPQCData(MPQCresults, fragmentData);
707
708 // print intermediate short-range results
709 {
710 LOG(1, "INFO: Parsing fragment files from " << params.path.get() << ".");
711 printReceivedMPQCResults(
712 MPQCresults,
713 fragmentData,
714 params.path.get(),
715 getNoAtomsFromAdjacencyFile(params.path.get()));
716 }
717
718#ifdef HAVE_VMG
719 if (params.DoLongrange.get()) {
720 // Phase Four: obtain more ids
721 requestIds(controller, params, fragmentData.size()+1);
722 RunService(io_service, "Requesting ids");
723
724 // Phase Five: create VMGJobs
725 if (!createLongRangeJobs(controller, params, fragmentData))
726 return Action::failure;
727 RunService(io_service, "Adding VMGJobs");
728
729 // Phase Six: calculate result
730 WaitforResults(io_service, controller, params, fragmentData.size()+1);
731 controller.receiveResults(params.host.get(), params.port.get());
732 RunService(io_service, "Requesting long-range results");
733 std::vector<FragmentResult::ptr> VMGresults = controller.getReceivedResults();
734 ASSERT( MPQCresults.size()+1 == VMGresults.size(),
735 "FragmentationFragmentationAutomationAction::performCall() - number of MPQCresultd and VMGresults don't match.");
736
737 // Final phase: print result
738 {
739 LOG(1, "INFO: Parsing fragment files from " << params.path.get() << ".");
740 printReceivedMPQCResults(
741 MPQCresults,
742 fragmentData,
743 params.path.get(),
744 getNoAtomsFromAdjacencyFile(params.path.get()));
745 }
746 }
747#endif
748 size_t Exitflag = controller.getExitflag();
749
750 return (Exitflag == 0) ? Action::success : Action::failure;
751}
752
753Action::state_ptr FragmentationFragmentationAutomationAction::performUndo(Action::state_ptr _state) {
754 return Action::success;
755}
756
757Action::state_ptr FragmentationFragmentationAutomationAction::performRedo(Action::state_ptr _state){
758 return Action::success;
759}
760
761bool FragmentationFragmentationAutomationAction::canUndo() {
762 return false;
763}
764
765bool FragmentationFragmentationAutomationAction::shouldUndo() {
766 return false;
767}
768/** =========== end of function ====================== */
Note: See TracBrowser for help on using the repository browser.