source: src/Actions/FragmentationAction/FragmentationAutomationAction.cpp@ 27afbf

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 27afbf was 27afbf, checked in by Frederik Heber <heber@…>, 12 years ago

FragmentationAutomationAction now prints table of results for energy, force, and times.

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