source: src/Actions/FragmentationAction/FragmentationAutomationAction.cpp@ 8a0e78

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 Candidate_v1.7.0 Candidate_v1.7.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 8a0e78 was a3a382, checked in by Frederik Heber <heber@…>, 13 years ago

FragmentationAutomationAction creates instance of print(Orthogonal)Sum.

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