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

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

Outsourced parsing of fragment files into ParseFragmentJobsAction from FragmentationAutomationAction.

  • Property mode set to 100644
File size: 24.1 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 headers that implement a archive in simple text format
42#include <boost/archive/text_oarchive.hpp>
43#include <boost/archive/text_iarchive.hpp>
44
45
46#include <boost/mpl/remove.hpp>
47#include <boost/lambda/lambda.hpp>
48
49#include <iostream>
50
51#include "CodePatterns/Assert.hpp"
52#include "CodePatterns/Info.hpp"
53#include "CodePatterns/Log.hpp"
54#include "JobMarket/Jobs/FragmentJob.hpp"
55
56#include "Fragmentation/Automation/createMatrixNrLookup.hpp"
57#include "Fragmentation/Automation/extractJobIds.hpp"
58#include "Fragmentation/Automation/FragmentationChargeDensity.hpp"
59#include "Fragmentation/Automation/FragmentJobQueue.hpp"
60#include "Fragmentation/Automation/FragmentationResults.hpp"
61#include "Fragmentation/Automation/MPQCFragmentController.hpp"
62#include "Fragmentation/Automation/VMGDebugGridFragmentController.hpp"
63#include "Fragmentation/Automation/VMGFragmentController.hpp"
64#include "Fragmentation/EnergyMatrix.hpp"
65#include "Fragmentation/ForceMatrix.hpp"
66#include "Fragmentation/Fragmentation.hpp"
67#include "Fragmentation/Homology/HomologyContainer.hpp"
68#include "Fragmentation/Homology/HomologyGraph.hpp"
69#include "Fragmentation/KeySet.hpp"
70#include "Fragmentation/KeySetsContainer.hpp"
71#include "Fragmentation/SetValues/Fragment.hpp"
72#include "Fragmentation/SetValues/Histogram.hpp"
73#include "Fragmentation/SetValues/IndexedVectors.hpp"
74#include "Fragmentation/Summation/IndexSetContainer.hpp"
75#include "Fragmentation/Summation/writeTable.hpp"
76#include "Graph/DepthFirstSearchAnalysis.hpp"
77#include "Jobs/MPQCJob.hpp"
78#include "Jobs/MPQCData.hpp"
79#include "Jobs/MPQCData_printKeyNames.hpp"
80#ifdef HAVE_VMG
81#include "Jobs/VMGDebugGridJob.hpp"
82#include "Jobs/VMGJob.hpp"
83#include "Jobs/VMGData.hpp"
84#include "Jobs/VMGDataFused.hpp"
85#include "Jobs/VMGDataMap.hpp"
86#include "Jobs/VMGData_printKeyNames.hpp"
87#endif
88#include "World.hpp"
89
90#include <fstream>
91#include <iostream>
92#include <string>
93#include <vector>
94
95#include <boost/mpl/for_each.hpp>
96
97#include "Actions/FragmentationAction/FragmentationAutomationAction.hpp"
98
99using namespace MoleCuilder;
100
101// and construct the stuff
102#include "FragmentationAutomationAction.def"
103#include "Action_impl_pre.hpp"
104/** =========== define the function ====================== */
105
106class controller_AddOn;
107
108// needs to be defined for using the FragmentController
109controller_AddOn *getAddOn()
110{
111 return NULL;
112}
113
114/** Helper function to get number of atoms somehow.
115 *
116 * Here, we just parse the number of lines in the adjacency file as
117 * it should correspond to the number of atoms, except when some atoms
118 * are not bonded, but then fragmentation makes no sense.
119 *
120 * @param path path to the adjacency file
121 */
122size_t getNoAtomsFromAdjacencyFile(const std::string &path)
123{
124 size_t NoAtoms = 0;
125
126 // parse in special file to get atom count (from line count)
127 std::string filename(path);
128 filename += FRAGMENTPREFIX;
129 filename += ADJACENCYFILE;
130 std::ifstream adjacency(filename.c_str());
131 if (adjacency.fail()) {
132 LOG(0, endl << "getNoAtomsFromAdjacencyFile() - Unable to open " << filename << ", is the directory correct?");
133 return false;
134 }
135 std::string buffer;
136 while (getline(adjacency, buffer))
137 NoAtoms++;
138 LOG(1, "INFO: There are " << NoAtoms << " atoms.");
139
140 return NoAtoms;
141}
142
143
144
145/** Place results from FragmentResult into EnergyMatrix and ForceMatrix.
146 *
147 * @param fragmentData MPQCData resulting from the jobs
148 * @param MatrixNrLookup Lookup up-map from job id to fragment number
149 * @param FragmentCounter total number of fragments
150 * @param NoAtoms total number of atoms
151 * @param Energy energy matrix to be filled on return
152 * @param Force force matrix to be filled on return
153 * @return true - everything ok, false - else
154 */
155bool putResultsintoMatrices(
156 const std::map<JobId_t, MPQCData> &fragmentData,
157 const std::map< JobId_t, size_t > &MatrixNrLookup,
158 const size_t FragmentCounter,
159 const size_t NoAtoms,
160 EnergyMatrix &Energy,
161 ForceMatrix &Force)
162{
163 for (std::map<JobId_t, MPQCData>::const_iterator dataiter = fragmentData.begin();
164 dataiter != fragmentData.end(); ++dataiter) {
165 const MPQCData &extractedData = dataiter->second;
166 const JobId_t &jobid = dataiter->first;
167 std::map< JobId_t, size_t >::const_iterator nriter = MatrixNrLookup.find(jobid);
168 ASSERT( nriter != MatrixNrLookup.end(),
169 "putResultsintoMatrices() - MatrixNrLookup does not contain id "
170 +toString(jobid)+".");
171 // place results into EnergyMatrix ...
172 {
173 MatrixContainer::MatrixArray matrix;
174 matrix.resize(1);
175 matrix[0].resize(1, extractedData.energies.total);
176 if (!Energy.AddMatrix(
177 std::string("MPQCJob ")+toString(jobid),
178 matrix,
179 nriter->second)) {
180 ELOG(1, "Adding energy matrix failed.");
181 return false;
182 }
183 }
184 // ... and ForceMatrix (with two empty columns in front)
185 {
186 MatrixContainer::MatrixArray matrix;
187 const size_t rows = extractedData.forces.size();
188 matrix.resize(rows);
189 for (size_t i=0;i<rows;++i) {
190 const size_t columns = 2+extractedData.forces[i].size();
191 matrix[i].resize(columns, 0.);
192 // for (size_t j=0;j<2;++j)
193 // matrix[i][j] = 0.;
194 for (size_t j=2;j<columns;++j)
195 matrix[i][j] = extractedData.forces[i][j-2];
196 }
197 if (!Force.AddMatrix(
198 std::string("MPQCJob ")+toString(jobid),
199 matrix,
200 nriter->second)) {
201 ELOG(1, "Adding force matrix failed.");
202 return false;
203 }
204 }
205 }
206 // add one more matrix (not required for energy)
207 MatrixContainer::MatrixArray matrix;
208 matrix.resize(1);
209 matrix[0].resize(1, 0.);
210 if (!Energy.AddMatrix(std::string("MPQCJob total"), matrix, FragmentCounter))
211 return false;
212 // but for energy because we need to know total number of atoms
213 matrix.resize(NoAtoms);
214 for (size_t i = 0; i< NoAtoms; ++i)
215 matrix[i].resize(2+NDIM, 0.);
216 if (!Force.AddMatrix(std::string("MPQCJob total"), matrix, FragmentCounter))
217 return false;
218
219 return true;
220}
221
222/** Print MPQCData from received results.
223 *
224 * @param fragmentData MPQCData resulting from the jobs, associated to job id
225 * @param MatrixNrLookup map with enumerated jobids
226 * @param KeySet KeySets of all (non-hydrogen) atoms
227 * @param ForceKeySet KeySets of all atoms except those added by saturation
228 * @param NoAtoms total number of atoms
229 */
230bool printReceivedMPQCResults(
231 const std::map<JobId_t, MPQCData> &fragmentData,
232 const std::map< JobId_t, size_t > MatrixNrLookup,
233 KeySetsContainer KeySet,
234 const KeySetsContainer &ForceKeySet,
235 size_t NoAtoms)
236{
237 size_t FragmentCounter = MatrixNrLookup.size();
238
239 // place results into maps
240 EnergyMatrix Energy;
241 ForceMatrix Force;
242 if (!putResultsintoMatrices(fragmentData, MatrixNrLookup, FragmentCounter, NoAtoms, Energy, Force))
243 return false;
244
245 if (!Energy.InitialiseIndices()) return false;
246
247 if (!Force.ParseIndices(ForceKeySet)) return false;
248
249 // combine all found data
250 if (!KeySet.ParseManyBodyTerms()) return false;
251
252 EnergyMatrix EnergyFragments;
253 ForceMatrix ForceFragments;
254 if (!EnergyFragments.AllocateMatrix(Energy.Header, Energy.MatrixCounter, Energy.RowCounter, Energy.ColumnCounter)) return false;
255 if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return false;
256
257 if(!Energy.SetLastMatrix(0., 0)) return false;
258 if(!Force.SetLastMatrix(0., 2)) return false;
259
260 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
261 // --------- sum up energy --------------------
262 LOG(1, "INFO: Summing energy of order " << BondOrder+1 << " ...");
263 if (!EnergyFragments.SumSubManyBodyTerms(Energy, KeySet, BondOrder)) return false;
264 if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return false;
265
266 // --------- sum up Forces --------------------
267 LOG(1, "INFO: Summing forces of order " << BondOrder+1 << " ...");
268 if (!ForceFragments.SumSubManyBodyTerms(Force, KeySet, BondOrder)) return false;
269 if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder, 1.)) return false;
270 }
271
272 // for debugging print resulting energy and forces
273 LOG(1, "INFO: Resulting energy is " << Energy.Matrix[ FragmentCounter ][0][0]);
274 std::stringstream output;
275 for (int i=0; i< Force.RowCounter[FragmentCounter]; ++i) {
276 for (int j=0; j< Force.ColumnCounter[FragmentCounter]; ++j)
277 output << Force.Matrix[ FragmentCounter ][i][j] << " ";
278 output << "\n";
279 }
280 LOG(1, "INFO: Resulting forces are " << std::endl << output.str());
281
282 return true;
283}
284
285/** Print MPQCData from received results.
286 *
287 * @param fragmentData MPQCData resulting from the jobs, associated to job id
288 * @param KeySetFilename filename with keysets to associate forces correctly
289 * @param NoAtoms total number of atoms
290 */
291bool printReceivedMPQCResults(
292 const std::map<JobId_t, MPQCData> &fragmentData,
293 const std::string &KeySetFilename,
294 size_t NoAtoms)
295{
296 // create a vector of all job ids
297 std::vector<JobId_t> jobids;
298 std::transform(fragmentData.begin(),fragmentData.end(),
299 std::back_inserter(jobids),
300 boost::bind( &std::map<JobId_t,MPQCData>::value_type::first, boost::lambda::_1 )
301 );
302
303 // create lookup from job nr to fragment number
304 size_t FragmentCounter = 0;
305 const std::map< JobId_t, size_t > MatrixNrLookup =
306 createMatrixNrLookup(jobids, FragmentCounter);
307
308 // initialise keysets
309 KeySetsContainer KeySet;
310 parseKeySetFile(KeySet, KeySetFilename, FragmentCounter, NonHydrogenKeySets);
311 KeySetsContainer ForceKeySet;
312 parseKeySetFile(ForceKeySet, KeySetFilename, FragmentCounter, HydrogenKeySets);
313
314 return printReceivedMPQCResults(
315 fragmentData,
316 MatrixNrLookup,
317 KeySet,
318 ForceKeySet,
319 NoAtoms);
320}
321
322/** Print MPQCData from received results.
323 *
324 * @param fragmentData MPQCData resulting from the jobs, associated to job id
325 * @param KeySet KeySets of all (non-hydrogen) atoms
326 * @param ForceKeySet KeySets of all atoms except those added by saturation
327 * @param NoAtoms total number of atoms
328 */
329bool printReceivedMPQCResults(
330 const std::map<JobId_t, MPQCData> &fragmentData,
331 const KeySetsContainer &KeySet,
332 const KeySetsContainer &ForceKeySet,
333 size_t NoAtoms)
334{
335 // create a vector of all job ids
336 std::vector<JobId_t> jobids;
337 std::transform(fragmentData.begin(),fragmentData.end(),
338 std::back_inserter(jobids),
339 boost::bind( &std::map<JobId_t,MPQCData>::value_type::first, boost::lambda::_1 )
340 );
341
342 // create lookup from job nr to fragment number
343 size_t FragmentCounter = 0;
344 const std::map< JobId_t, size_t > MatrixNrLookup =
345 createMatrixNrLookup(jobids, FragmentCounter);
346
347 return printReceivedMPQCResults(
348 fragmentData,
349 MatrixNrLookup,
350 KeySet,
351 ForceKeySet,
352 NoAtoms);
353}
354
355void writeToFile(const std::string &filename, const std::string contents)
356{
357 std::ofstream tablefile(filename.c_str());
358 tablefile << contents;
359 tablefile.close();
360}
361
362/** Print MPQCData from received results.
363 *
364 * @param results summed up results container
365 */
366void printReceivedFullResults(
367 const FragmentationResults &results)
368{
369 // print tables (without eigenvalues, they go extra)
370 {
371 typedef boost::mpl::remove<MPQCDataEnergyVector_t, MPQCDataFused::energy_eigenvalues>::type
372 MPQCDataEnergyVector_noeigenvalues_t;
373 const std::string energyresult =
374 writeTable<MPQCDataEnergyMap_t, MPQCDataEnergyVector_noeigenvalues_t >()(
375 results.Result_Energy_fused, results.getMaxLevel());
376 LOG(0, "Energy table is \n" << energyresult);
377 std::string filename;
378 filename += FRAGMENTPREFIX + std::string("_Energy.dat");
379 writeToFile(filename, energyresult);
380 }
381
382 if (results.Result_LongRange_fused.size() >= (results.getMaxLevel()-2))
383 {
384 const std::string gridresult =
385 writeTable<VMGDataMap_t, VMGDataVector_t >()(
386 results.Result_LongRange_fused, results.getMaxLevel(), 2);
387 LOG(0, "VMG table is \n" << gridresult);
388 std::string filename;
389 filename += FRAGMENTPREFIX + std::string("_VMGEnergy.dat");
390 writeToFile(filename, gridresult);
391 }
392
393 if (results.Result_LongRangeIntegrated_fused.size() >= (results.getMaxLevel()-2))
394 {
395 const std::string gridresult =
396 writeTable<VMGDataLongRangeMap_t, VMGDataLongRangeVector_t >()(
397 results.Result_LongRangeIntegrated_fused, results.getMaxLevel(), 2);
398 LOG(0, "LongRange table is \n" << gridresult);
399 std::string filename;
400 filename += FRAGMENTPREFIX + std::string("_LongRangeEnergy.dat");
401 writeToFile(filename, gridresult);
402 }
403
404 {
405 const std::string eigenvalueresult;
406 LOG(0, "Eigenvalue table is \n" << eigenvalueresult);
407 std::string filename;
408 filename += FRAGMENTPREFIX + std::string("_Eigenvalues.dat");
409 writeToFile(filename, eigenvalueresult);
410 }
411
412 {
413 const std::string forceresult =
414 writeTable<MPQCDataForceMap_t, MPQCDataForceVector_t>()(
415 results.Result_Force_fused, results.getMaxLevel());
416 LOG(0, "Force table is \n" << forceresult);
417 std::string filename;
418 filename += FRAGMENTPREFIX + std::string("_Forces.dat");
419 writeToFile(filename, forceresult);
420 }
421 // we don't want to print grid to a table
422 {
423 // print times (without flops for now)
424 typedef boost::mpl::remove<
425 boost::mpl::remove<MPQCDataTimeVector_t, MPQCDataFused::times_total_flops>::type,
426 MPQCDataFused::times_gather_flops>::type
427 MPQCDataTimeVector_noflops_t;
428 const std::string timesresult =
429 writeTable<MPQCDataTimeMap_t, MPQCDataTimeVector_noflops_t >()(
430 results.Result_Time_fused, results.getMaxLevel());
431 LOG(0, "Times table is \n" << timesresult);
432 std::string filename;
433 filename += FRAGMENTPREFIX + std::string("_Times.dat");
434 writeToFile(filename, timesresult);
435 }
436}
437
438bool appendToHomologyFile(
439 const boost::filesystem::path &homology_file,
440 const FragmentationResults &results)
441{
442 /// read homology container (if present)
443 HomologyContainer homology_container;
444 if (boost::filesystem::exists(homology_file)) {
445 std::ifstream returnstream(homology_file.string().c_str());
446 if (returnstream.good()) {
447 boost::archive::text_iarchive ia(returnstream);
448 ia >> homology_container;
449 } else {
450 ELOG(2, "Failed to parse from " << homology_file.string() << ".");
451 }
452 returnstream.close();
453 } else {
454 LOG(2, "Could not open " << homology_file.string()
455 << ", creating empty container.");
456 }
457
458 /// append all fragments to a HomologyContainer
459 HomologyContainer::container_t values;
460 const size_t FragmentCounter = results.Result_perIndexSet_Energy.size();
461
462 // convert KeySetContainer to IndexSetContainer
463 IndexSetContainer::ptr ForceContainer(new IndexSetContainer(results.getForceKeySet()));
464 const IndexSetContainer::Container_t &Indices = results.getContainer();
465 const IndexSetContainer::Container_t &ForceIndices = ForceContainer->getContainer();
466 IndexSetContainer::Container_t::const_iterator iter = Indices.begin();
467 IndexSetContainer::Container_t::const_iterator forceiter = ForceIndices.begin();
468 /// go through all fragments
469 for (;iter != Indices.end(); ++iter, ++forceiter) // go through each IndexSet
470 {
471 /// create new graph entry in HomologyContainer which is (key,value) type
472 LOG(1, "INFO: Creating new graph with " << **forceiter << ".");
473 HomologyGraph graph(**forceiter);
474 LOG(2, "DEBUG: Created graph " << graph << ".");
475 const IndexSet::ptr &index = *iter;
476 HomologyContainer::value_t value;
477
478 // obtain fragment as key
479 std::map<IndexSet::ptr, std::pair< MPQCDataFragmentMap_t, MPQCDataFragmentMap_t> >::const_iterator fragmentiter
480 = results.Result_perIndexSet_Fragment.find(index);
481 ASSERT( fragmentiter != results.Result_perIndexSet_Fragment.end(),
482 "appendToHomologyFile() - cannot find index "+toString(*index)
483 +" in FragmentResults.");
484 value.first = boost::fusion::at_key<MPQCDataFused::fragment>(fragmentiter->second.first);
485
486 // obtain energy as value
487 std::map<IndexSet::ptr, std::pair< MPQCDataEnergyMap_t, MPQCDataEnergyMap_t> >::const_iterator energyiter
488 = results.Result_perIndexSet_Energy.find(index);
489 ASSERT( energyiter != results.Result_perIndexSet_Energy.end(),
490 "appendToHomologyFile() - cannot find index "+toString(*index)
491 +" in FragmentResults.");
492 // value.second = boost::fusion::at_key<MPQCDataFused::energy_total>(energyiter->second.first); // values
493 value.second = boost::fusion::at_key<MPQCDataFused::energy_total>(energyiter->second.second); // contributions
494 values.insert( std::make_pair( graph, value) );
495 }
496 homology_container.insert(values);
497
498 LOG(1, "INFO: Listing all present atomic ids ...");
499 std::stringstream output;
500 for (World::AtomIterator iter = World::getInstance().getAtomIter();
501 iter != World::getInstance().atomEnd(); ++iter)
502 output << (*iter)->getId() << " ";
503 LOG(1, "INFO: { " << output.str() << "} .");
504
505 // for debugging: print container
506 LOG(1, "INFO: Listing all present homologies ...");
507 for (HomologyContainer::container_t::const_iterator iter =
508 homology_container.begin(); iter != homology_container.end(); ++iter) {
509 LOG(1, "INFO: graph " << iter->first << " has Fragment "
510 << iter->second.first << " and associated energy " << iter->second.second << ".");
511 }
512
513 /// store homology container again
514 std::ofstream outputstream(homology_file.string().c_str());
515 if (outputstream.good()) { // check if opened
516 boost::archive::text_oarchive oa(outputstream);
517 oa << homology_container;
518 if (outputstream.fail()) { // check if correctly written
519 LOG(1, "Failed to write to file " << homology_file.string() << ".");
520 return false;
521 } else
522 outputstream.close();
523 } else {
524 LOG(1, "Failed to open file " << homology_file.string()
525 << " for writing.");
526 return false;
527 }
528 return true;
529}
530
531Action::state_ptr FragmentationFragmentationAutomationAction::performCall() {
532 boost::asio::io_service io_service;
533
534 // TODO: Have io_service run in second thread and merge with current again eventually
535
536 size_t Exitflag = 0;
537 std::map<JobId_t, MPQCData> fragmentData;
538 {
539 const size_t NumberJobs = FragmentJobQueue::getInstance().size();
540 MPQCFragmentController mpqccontroller(io_service);
541 mpqccontroller.setHost(params.host.get());
542 mpqccontroller.setPort(params.port.get());
543 // Phase One: obtain ids
544 mpqccontroller.requestIds(NumberJobs);
545
546 // Phase Two: add MPQCJobs and send
547 mpqccontroller.addJobsFromQueue();
548 mpqccontroller.run();
549
550 // Phase Three: calculate result
551 mpqccontroller.waitforResults(NumberJobs);
552 mpqccontroller.getResults(fragmentData);
553
554 Exitflag += mpqccontroller.getExitflag();
555 }
556
557#ifdef HAVE_VMG
558 if (params.DoLongrange.get()) {
559 if ( World::getInstance().getAllAtoms().size() == 0) {
560 ELOG(1, "Please load the full molecule into the world before starting this action.");
561 return Action::failure;
562 }
563
564 // obtain combined charge density
565 FragmentationChargeDensity summedChargeDensity(
566 fragmentData,
567 FragmentJobQueue::getInstance().getKeySets());
568 const std::vector<SamplingGrid> full_sample = summedChargeDensity.getFullSampledGrid();
569 LOG(1, "INFO: There are " << fragmentData.size() << " short-range and "
570 << full_sample.size() << " level-wise long-range jobs.");
571
572 // Phase Four: obtain more ids
573 std::map<JobId_t, VMGData> longrangeData;
574 {
575 VMGFragmentController vmgcontroller(io_service);
576 vmgcontroller.setHost(params.host.get());
577 vmgcontroller.setPort(params.port.get());
578 const size_t NoJobs = fragmentData.size()+full_sample.size();
579 vmgcontroller.requestIds(NoJobs);
580
581 // Phase Five: create VMGJobs
582 const size_t near_field_cells = params.near_field_cells.get();
583 const size_t interpolation_degree = params.interpolation_degree.get();
584 if (!vmgcontroller.createLongRangeJobs(
585 fragmentData,
586 full_sample,
587 summedChargeDensity.getFragment(),
588 near_field_cells,
589 interpolation_degree))
590 return Action::failure;
591
592 // Phase Six: calculate result
593 vmgcontroller.waitforResults(NoJobs);
594 vmgcontroller.getResults(longrangeData);
595 ASSERT( NoJobs == longrangeData.size(),
596 "FragmentationFragmentationAutomationAction::performCall() - number of MPQCresults+"
597 +toString(full_sample.size())+"="+toString(NoJobs)
598 +" and VMGresults "+toString(longrangeData.size())+" don't match.");
599 Exitflag += vmgcontroller.getExitflag();
600 }
601
602 // remove full solution corresponding to full_sample from map (must be highest ids), has to be treated extra
603 std::map<JobId_t, VMGData>::iterator iter = longrangeData.end();
604 for (size_t i=0;i<full_sample.size();++i)
605 --iter;
606 std::map<JobId_t, VMGData>::iterator remove_iter = iter;
607 std::vector<VMGData> fullsolutionData;
608 for (; iter != longrangeData.end(); ++iter)
609 fullsolutionData.push_back(iter->second);
610 longrangeData.erase(remove_iter, longrangeData.end());
611
612 // Final phase: sum up and print result
613 FragmentationResults results(
614 fragmentData,
615 longrangeData,
616 FragmentJobQueue::getInstance().getKeySets(),
617 FragmentJobQueue::getInstance().getFullKeySets());
618 results(
619 fragmentData,
620 longrangeData,
621 fullsolutionData,
622 full_sample);
623 printReceivedFullResults(results);
624
625 // append all keysets to homology file
626 if ((Exitflag == 0) && (!params.homology_file.get().empty())) {
627 const boost::filesystem::path &homology_file = params.homology_file.get();
628 if (homology_file.string() != "") {
629 LOG(1, "INFO: Appending HomologyGraphs to file " << homology_file.string() << ".");
630 if (!appendToHomologyFile(homology_file, results))
631 Exitflag = 1;
632 }
633 }
634
635 std::map<JobId_t, std::string> debugData;
636 {
637 if (!full_sample.empty()) {
638 // create debug jobs for each level to print the summed-up potential to vtk files
639 VMGDebugGridFragmentController debugcontroller(io_service);
640 debugcontroller.setHost(params.host.get());
641 debugcontroller.setPort(params.port.get());
642 debugcontroller.requestIds(full_sample.size());
643 if (!debugcontroller.createDebugJobs(full_sample))
644 return Action::failure;
645 debugcontroller.waitforResults(full_sample.size());
646 debugcontroller.getResults(debugData);
647 Exitflag += debugcontroller.getExitflag();
648 }
649 }
650 }
651#else
652 // Final phase: print result
653 {
654 LOG(1, "INFO: Parsing fragment files from " << params.path.get() << ".");
655 printReceivedMPQCResults(
656 fragmentData,
657 FragmentJobQueue::getInstance().getKeySets(),
658 FragmentJobQueue::getInstance().getFullKeySets(),
659 World::getInstance().getAllAtoms().size()));
660 }
661#endif
662
663 // now clear all present jobs as we are done
664 FragmentJobQueue::getInstance().clear();
665
666 return (Exitflag == 0) ? Action::success : Action::failure;
667}
668
669Action::state_ptr FragmentationFragmentationAutomationAction::performUndo(Action::state_ptr _state) {
670 return Action::success;
671}
672
673Action::state_ptr FragmentationFragmentationAutomationAction::performRedo(Action::state_ptr _state){
674 return Action::success;
675}
676
677bool FragmentationFragmentationAutomationAction::canUndo() {
678 return false;
679}
680
681bool FragmentationFragmentationAutomationAction::shouldUndo() {
682 return false;
683}
684/** =========== end of function ====================== */
Note: See TracBrowser for help on using the repository browser.